In　this　paper　we　propose　and　analyze　a　second　order　accurate　(in　time)　numerical　scheme　for　the　square　phase　field　crystal　equation,　a　gradient　flow　modeling　crystal　dynamics　at　the　atomic　scale　in　space　but　on　diffusive　scales　in　time.　Its　primary　difference　with　the　standard　phase　field　crystal　model　is　an　introduction　of　the　4-Laplacian　term　in　the　free　energy　potential,　which　in　turn　leads　to　a　much　higher　degree　of　nonlinearity.　To　make　the　numerical　scheme　linear　while　preserving　the　nonlinear　energy　stability,　we　make　use　of　the　scalar　auxiliary　variable　(SAV)　approach,　in　which　a　second　order　backward　differentiation　formula　is　applied　in　the　temporal　stencil.　Meanwhile,　a　direct　application　of　the　SAV　method　faces　certain　difficulties,　due　to　the　involvement　of　the　4-Laplacian　term,　combined　with　a　derivation　of　the　lower　bound　of　the　nonlinear　energy　functional.　In　the　proposed　numerical　method,　an　appropriate　decomposition　for　the　physical　energy　functional　is　formulated,　so　that　the　nonlinear　energy　part　has　a　well-established　global　lower　bound,　and　the　rest　terms　lead　to　constant-coefficient　diffusion　terms　with　positive　eigenvalues.　In　turn,　the　numerical　scheme　could　be　very　efficiently　implemented　by　constant-coefficient　Poisson-like　type　solvers　(via　FFT),　and　energy　stability　is　established　by　introducing　an　auxiliary　variable,　and　an　optimal　rate　convergence　analysis　is　provided　for　the　proposed　SAV　method.　A　few　numerical　experiments　are　also　presented,　which　confirm　the　efficiency　and　accuracy　of　the　proposed　scheme.