The initial flow fields of the supersonic underexpanded square jets with different nozzle pressure ratios were investigated numerically. The governing equations of large eddy simulation (LES) for compressible flow have been employed and solved numerically with the combination of high-order hybrid schemes. The numerical result provides the 3D shock wave structures embedded in the underexpanded square jet. Moreover, the effect of nozzle pressure ratio on the 3D characteristic behavior of near-nozzle shock wave and the jet boundary are discussed, it is found that the near-nozzle shock wave structure suppresses the jet boundary at the corner position expansion along the diagonal direction, and accelerates the jet boundary near the sidewall position to move in normal direction to the nozzle side, finally resulting in a remarkable cross-shaped jet cross section at the downstream. Besides, the 3D behavior of stream wise vortices during the self-induced deformation of primary vortex loop is identified using the iso-vorticity field, and their role in the axis-switching behavior of supersonic square vortex loop is analyzed in detail. (C) 2017 Elsevier Ltd. All rights reserved.