  We utilize Density Matrix Renormalization Group (DMRG) and Fully Augmented
Matrix Product States (FAMPS) methods to investigate the Valence Bond Solid
(VBS) phase in the $J_1$-$J_2$ square lattice Heisenberg model. To
differentiate between the Columnar Valence Bond Solid (CVBS) and Plaquette
Valence Bond Solid (PVBS) phases, we introduce an anisotropy $\Delta_y$ in the
nearest neighboring coupling in the $y$-direction, aiming at detecting the
possible spontaneous rotational symmetry breaking in the VBS phase. In the
calculations, we push the bond dimension to as large as $D = 25000$ in FAMPS,
simulating systems at a maximum size of $14 \times 14$. With a careful
extrapolation of the truncation errors and appropriate finite-size scaling,
followed by finite $\Delta_y$ scaling analysis of the VBS dimer order
parameters, we identify the VBS phase as a PVBS type, meaning there is no
spontaneous rotational symmetry breaking in the VBS phase. This study not only
resolves the long-standing issue of the characterization of the VBS order in
the $J_1$-$J_2$ square lattice Heisenberg model but also highlights the
capabilities of FAMPS in the study of two-dimensional quantum many-body
systems.
