Spin-orbit-enhanced magnetic surface second-harmonic generation in Sr$_2$IrO$_4$

An anomalous optical second-harmonic generation (SHG) signal was previously reported in Sr$_2$IrO$_4$ and attributed to a hidden odd-parity bulk magnetic state. Here we investigate the origin of this SHG signal using a combination of bulk magnetic susceptibility, magnetic-field-dependent SHG rotational anisotropy, and overlapping wide-field SHG imaging and atomic force microscopy measurements. We find that the anomalous SHG signal exhibits a two-fold rotational symmetry as a function of in-plane magnetic field orientation that is associated with a crystallographic distortion. We also show a change in SHG signal across step edges that tracks the bulk antiferromagnetic stacking pattern. While we do not rule out the existence of hidden order in Sr$_2$IrO$_4$, our results altogether show that the anomalous SHG signal in parent Sr$_2$IrO$_4$ originates instead from a surface-magnetization-induced electric-dipole process that is enhanced by strong spin-orbit coupling.

The key feature of the SHG data in Sr 2 IrO 4 is the onset of a new radiation process below the Néel temperature (T N ∼ 230 K) that lowers the rotational symmetry about the c axis from C 4 to C 1 . This is incompatible with its reported Néel structure [1,[17][18][19], in which a canting-induced net ferromagnetic moment in each layer is stacked along the c axis in a − + +− order so as to preserve C 2 symmetry [ Fig. 1(a)]. Rather, it was found to be consistent with electric-dipole (ED) radiation from a noncentrosymmetric 2 /m or m1 magnetic point group [7], which bears the symmetries of a magnetoelectric loop-current phase [20]. However, alternative explanations were subsequently proposed [12], including laserinduced rearrangement of the magnetic stacking and enhanced sensitivity to surface rather than bulk magnetic order.
Studying the field dependence of the anomalous C 1 SHG signal can help reveal its origin. Previous works showed that applying an in-plane magnetic field H > H c (∼200 mT for temperatures well below T N ) induces a metamagnetic transition to a C 2 -breaking + + ++ state [ Fig. 1(a)], whose net ferromagnetic moment rotates with the in-plane field orientation [21]. In contrast, the proposed magnetoelectric order parameter can only couple linearly to H in the presence of an additional electric field [22]. Here we report field-dependent SHG rotational anisotropy (RA), magnetic susceptibility, as well as overlapping atomic force and SHG microscopy measurements on parent Sr 2 IrO 4 . While we do not rule out the presence of hidden order in Sr 2 IrO 4 [8][9][10], our data show that the C 1 SHG signal arises from an unexpectedly strong surface-magnetization-induced ED process in the − + +− state, rather than from laser-induced or intrinsic hidden orders.
Our experiments were performed on cleaved single crystals of Sr 2 IrO 4 [23]. The geometry of our magneticfield-tunable RA-SHG setup is shown in Fig. 1(b). In zero field, the low temperature RA (φ-dependence) pattern is produced from the interference between a C 4 crystallographic electric-quadrupole (EQ) process and the anomalous C 1 process [ Fig. 1(c)]. By tracking the intensity of the strongest lobe versus temperature, we confirm an order-parameter-like onset of the latter near T = 230 K in agreement with previous work [7]. In the presence of a weak in-plane field (H < H c ), no change is observed in the high-temperature EQ intensity, but there is slight enhancement of the low-temperature SHG intensity along with the emergence of a peak structure just below T N [ Fig. 1(d)]. This is reminiscent of bulk magnetization data acquired under similar field conditions, which exhibits a ferromagnetic upturn at T N followed by a sharp drop as the − + +− state is stabilized [inset Fig. 1(d)]. A closer comparison of their temperature dependence suggests that the low-field SHG curve is a superposition of the zero-field curve and a new contribution that tracks the bulk magnetization. This is further corroborated by the high-field (H > H c ) SHG data [ Fig. 1(e)] that shows a large intensity increase relative to the low-field case and a disappearance of the peak struc-  1. (a) Intralayer Néel order of the J eff = 1/2 pseudospins and the c-axis stacking pattern of the net in-plane moments (purple arrows) at H = 0 (− + +−) and H > Hc (+ + ++) below TN . (b) Schematic of RA-SHG setup. The angle of incidence (θ) is fixed to 10 • while the scattering plane angle (φ) varied. The fundamental input (ω) and SH output (2ω) beams can be selected as P -or S -polarized. An in-plane magnetic field H (green arrow) can be applied at varying angles β relative to [100]. Temperature dependence of P in-S out SHG intensities acquired in a field-cooled warming scheme for (c) H = 0 mT, (d) 75 mT, and (e) 370 mT, with φ fixed at the lobe of maximum intensity. The curves are normalized to the room-temperature value of the zero field data. Insets: (c) RA data for P in-S out polarization geometry measured at T = 80 K and H = 0 T. As reported in Ref. [7], this pattern is produced by the interference between a C4 crystallographic EQ term with a C1 magnetic ED term as illustrated on the right. Filled and white lobes denote opposite optical phase. The relative lobe intensities are highly sensitive to small changes in the EQ and ED tensor elements, which are found to vary slightly between this work and Ref. [7]. The data in Fig. 1 reveal three contributions to the field-dependent SHG response: a crystallographic EQ term and an anomalous C 1 term that both persist in zero field, plus a bulk magnetization-dependent term that is strongly enhanced in the + + ++ state. This eliminates the possibility of the sample already being in a laser-induced + + ++ state at zero field. To rule out the possibility of alternative laser-induced stacking patterns [12], we repeated our zero-field measurements on a fresh sample using two orders of magnitude weaker laser fluence and observed the same trend as in Fig. 1(c) [23]. Furthermore, we performed time-resolved SHG experiments and observed a strong suppression of the SHG intensity after a pump pulse, rather than an enhancement [23]. These tests show that the C 1 term arises from an intrinsic rather than laser-induced effect.
To examine how the anomalous C 1 term varies with magnetic field orientation, we collected RA-SHG patterns as a function of the in-plane field angle (β) with |H| fixed at 370 mT. Figure 2(a) shows that the RA patterns continue to exhibit C 1 symmetry in the + + ++ state and rotate by 180°upon reversal of the field direction (H → −H). This implies that both the bulk magnetization-dependent term and the anomalous C 1 term couple linearly to a magnetic field. By summing the RA-SHG data over all β to isolate the β-independent symmetries [ Fig. 2(b)], we confirm that no C 1 component remains, which is evidence against the magnetoelectric order parameter interpretation of the SHG response [23]. Surprisingly, despite the reported tetragonal symmetry of the lattice, the summed data exposes a twofold rotational symmetry, indicating an inequivalent response to fields along [100] and [010]. In Sr 2 IrO 4 , coupling of the pseudospins to an orthorhombically deformed lattice below T N naturally induces uniaxial anisotropy [21,29], which can lead to a twofold symmetry in β when H is below or near H c [30,31]. However, our SHG experiments were performed in a high-field regime where the bulk magnetization is constant and rotates rigidly with the field direction [21,23]. These observations were reproduced across multiple samples and up to temperatures as high as 220 K, where pseudospin-lattice coupling plays a smaller role [21,23]. This points to an a-b symmetry breaking distortion already existing above T N .
The nonlinear susceptibility tensors governing the bulk magnetization-dependent and anomalous SHG processes must both be c-type (time-reversal odd) by virtue of their observed linear coupling to H. Since the centrosymmetric 2 /m magnetic point group of the + + ++ state forbids ED SHG, it was proposed [12] that the former can be associated with a bulk magnetic-dipole (MD) process of the type where an incident electromagnetic field with frequency ω induces an oscillatory polarization at 2ω. To fit our high-field data [ Fig. 2(a)], we express the axial c tensor χ is a polar i tensor (time-reversal even) that respects the high temperature crystallographic point group and M = (M cos β, M sin β, 0) is the static magnetization of the sample [23]. This term is then coherently added to a constant crystallographic EQ term. A good simultaneous fit to the entire set of field-dependent RA patterns breaks tetragonal symmetry [23]. This reflects our observation that the magnetic field response along a and b differ.
Attributing the Néel-order-induced contribution to a c-type MD process naturally explains why there is strong bulk SHG from the + + ++ state but not the − + +− state. Since the MD susceptibility of each layer is time-reversal odd, layers with opposite magnetization radiate 180°out of phase, leading to overall destructive (constructive) interference for antiferromagnetic (ferromagnetic) stacking. It has been shown that − + +− stacking is stabilized by competing first and second nearest-neighbor interlayer exchange coupling [32]. Therefore, although we excluded bulk + + ++ stacking as the source of the zero-field SHG, one could posit that crystal termination forces a + + ++ stacking near the surface of a bulk − + +− ordered sample. In this scenario, we expect the zero-field RA patterns to map onto the high-field patterns by simply scaling χ MD(c) ijk (M). However, our fitting suggests this is not the case [23], leaving pure surface SHG as the remaining plausible candidate.
Since inversion symmetry is naturally broken at the surface of Sr 2 IrO 4 , a zero-field C 1 term can in principle arise from a surface-magnetization-induced c-type ED SHG process where χ ED(c) s,ijk (M s ) can be expressed in terms of the surface magnetization as χ ED(i) s,ijkl M s,l [23]. To explore this possibility, we performed SHG microscopy in conjunction with atomic force microscopy on ultraflat (001) surfaces of Sr 2 IrO 4 prepared by cleaving. These surfaces appear optically pristine in standard white light microscopy as well as wide-field SHG imaging at room temperature  [7], being careful to account for the shifted illumination between images, shows fixed boundaries similar to our data, suggestive of pinning to structural features.
To search for direct correlations between the spatial distribution of SHG patches and structural features, we show a large area topographic survey of the same sample in Fig. 3(f), which is a composite image of height maps assembled from several atomic force microscopy scans. The maps reveal an atomically smooth surface that is interrupted by long boundaries where the height changes by either a monolayer (∼0.65 nm) or bilayer (∼1.3 nm) [33]. We overlay these boundaries atop the SHG image in Fig. 3(e), with monolayer and bilayer steps marked by solid and dashed lines, respectively. We find that each topographical step correlates with the boundary of an SHG patch. But interestingly, the reverse is not true; there are clear patch boundaries that are not associated with topographical steps.
To understand these observations, we focus on regions where four distinct SHG patches intersect. Figure 4(a) shows a representative example. By performing scanning RA measurements over this area, we find that the four patches correspond to four different C 1 order parameter orientations denoted by arrows. Across both the 1 → 4 and 2 → 3 boundaries, the surface is atomically flat with no step edges, but the order parameter undergoes a 90°r otation [ Fig. 4(b)]. This is most likely caused by a crys- tallographic twin boundary [23] in which the a and b axes are interchanged. The presence of such twinning has been previously proposed to explain neutron diffraction results [17,18]. By measuring the magnetic-fielddependent RA response from each of the twin regions, we see evidence of a-b symmetry breaking, as displayed in Fig. 2(a), with axes switched between the regions [23]. Moreover, with repeated thermal cycling across T N , each region is observed to randomly undergo 0°or 180°, but never 90°or 270°, order parameter rotations [ Fig. 4(c)], confirming this interpretation.
On the other hand, upon crossing from region 1 → 2 or from region 3 → 4, one traverses a bilayer step edge that coincides with a 180°rotation of the order parameter [ Fig. 4(b)]. This is consistent with surface-magnetization-induced SHG from the established − + +− stacking since M s reverses every two layers. It also further rules out laser-induced + + ++ and − + −+ stacking that was previously proposed to explain the SHG data [12]. Upon successive thermal cycling across T N , the order parameter orientations in regions 1 and 2 are always anti-correlated, as is the case for regions 3 and 4 [ Fig. 4(c)], whereas no correlation exists across the twin boundary. This indicates that regions 1 and 2 form part of a single magnetic domain with − + +− stacking, and that regions 3 and 4 form part of another single magnetic domain that is 90°rotated. We note that a survey of all patches within the field of view shown in Figs. 3(e) and 3(f) reveal that most of the terrace steps correlate with twin boundaries, suggesting that Sr 2 IrO 4 crystals tend to fracture along them. We could not find any monolayer steps that were not associated with a twin boundary.
Surface-magnetization-induced ED SHG is typically weaker than its crystallographic counterpart. This is well documented for thin-film 3d transition metals such as Ni(110) [34] or Fe(110) [35], for which the ratio The canted Néel structure of Sr 2 IrO 4 , which is stabilized by a combination of conventional superexchange and Dzyaloshinskii-Moriya interactions [36], results in a small M s . Since the surface magnetization of Sr 2 IrO 4 (∼8 × 10 −2 µ B /Ir) [37] is approximately an order of magnitude weaker than that of Ni or Fe films [38], one naively expects a much smaller value of ρ. Previous SHG experiments on Sr 2 IrO 4 showed that the crystallographic contribution predominantly arises from a bulk χ EQ(i) process rather than a χ ED(i) s process [39], in keeping with weak inversion symmetry breaking at the surface due to its quasi-2D structure. Therefore, the existence of an appreciable χ ED(c) s (M s ) to χ EQ(i) ratio is surprising and explains why it was discounted in Ref. [7]. To understand this phenomenon, we appeal to a perturbative calculation of surface-magnetization-induced SHG [40], which showed that ρ ≈ M s Ab, where A is the spin-orbit coupling for the free atom [41] and b is a constant that depends on the band energies of the material. Since A increases significantly going from 3d to 5d transition metals, it makes sense that Sr 2 IrO 4 can still exhibit a large ρ despite M s being small.
Recent polarized neutron diffraction [8] and torque magnetometry [10] measurements have suggested the presence of a bulk magnetoelectric loop-current order in Sr 2 IrO 4 . There is also evidence that it survives and onsets above T N in Rh-doped samples [7][8][9][10]. However, these works point toward a magnetoelectric loop-current order that is ferroically stacked along the c axis, which is incompatible with our observation of an SHG signal that switches sign every two layers. While our work does not rule out the existence of such a hidden order, it firmly establishes that the anomalous SHG signal in parent Sr 2 IrO 4 is instead dominated by a process induced by surface magnetization. More generally, our results show that magnetization-induced SHG can be significantly enhanced in strongly spin-orbit coupled materials.

I. Experimental methods
Single crystals of Sr 2 IrO 4 were grown using an established flux method [1]. Dry starting powders of SrCO 3 , IrO 2 , and anhydrous SrCl 2 were placed in a 2:1:5.5 molar ratio within a 100 mL Pt crucible with a lid. The crucible was slowly heated to 1370°C in air, soaked for 5 hours, cooled to 850°C at 6°C per hour, and then furnace-cooled to room temperature. Single crystals were then obtained by dissolving the excess flux in deionized water. Magnetization data were collected on a cleaved and polished sample mounted on a quartz paddle with GE varnish, and measured with a Quantum Design MPMS3 system in the vibrating sample magnetometer mode. Samples for SHG were affixed to an oxygen-free high thermal conductivity copper mount using a small amount of Torr seal on the underside of the sample, then cleaved along the (001) face prior to measurement and pumped down below 10 −7 torr in a continuous flow optical cryostat. The crystal orientation was checked by x-ray diffraction. We found consistent results for the SHG magnetic field dependence across multiple samples. We note that the essential features of the zero-field SHG, namely, the presence of a C 1 term that onsets below ∼230 K, agree with prior reports [2]. However, the RA patterns differ in the relative magnitude of adjacent lobes in the S out geometries and the rotation away from the high symmetry directions for EQ SHG, which are markers of mirror symmetry breaking. These differences may stem from the high sensitivity of Sr 2 IrO 4 to slight variations in growth conditions [3]. The data presented here are representative of samples grown by the methods in Ref. [1].
RA-SHG measurements were performed using a rotating scattering plane technique [4] with 100 fs laser pulses at 800 nm from a Ti:sapphire amplifier (100 kHz repetition rate). A schematic of the experimental geometry is shown in Fig. S1. The fluence was fixed at 2 mJ/cm 2 with an angle of incidence of θ = 10°, except for the low-fluence Ti:sapphire oscillator measurements in Fig. S2(a), which used 10 µJ/cm 2 . In-plane magnetic fields were applied using SmCo permanent magnets placed adjacent to the sample in a home-built rotator apparatus [ Fig. S1], which allows for free rotation of the magnets about the cylindrical sample holder. The rotator can house several 3/8" × 3/8" × 1/8" SmCo magnets (T C ∼ 1000 K), with field strength adjustable by changing the number of magnets and their spacing. The magnetic field strengths were measured at room temperature with a gauss meter, and we assume constant magnetic field within our 80 K to 300 K range. The maximum field we achieved with SmCo magnets was 370 mT. The in-plane field direction is manually tunable using permanent magnets fixed to the outside rim of the cryostat. Wide-field SHG imaging was performed at θ ∼ 3°with a fluence of 4 mJ/cm 2 in P in -P out geometry.
H ω 2ω FIG. S1. Image of the home-built magnet rotator setup. A cylindrical copper sample mount holds an aluminum ring in place and allows for free rotation of the magnets about the cylinder axis.

II. Absence of laser-induced magnetic rearrangement
We examined the possibility that the SHG onset ∼230 K in Sr 2 IrO 4 is due to a metastable magnetic order induced by high-field laser pulses as suggested in Ref. [5]. To do so, we cleaved a fresh (not previously irradiated) Sr 2 IrO 4 sample and measured the temperature-dependent SHG intensity at a fluence of ∼10 µJ/cm 2 , 200 times lower than that used for the data shown in the main text, using a Ti:sapphire oscillator. To achieve sufficient signal-to-noise, we used 20 min integration times and pixel binning with an electron-multiplying charge-coupled device at a single angle of incidence where the intensity was at a maximum. Figure S2(a) shows that the SHG onset below ∼230 K is still clearly present. A laser-induced metastable stacking re-arrangement is therefore an unlikely explanation for the anomalous C 1 SHG term.
To rule out the possibility of a transient laser-induced effect, we performed an ultrafast pump-probe SHG experiment at 80 K to track the temporal change in SHG intensity following a pump pulse. Figure S2(b) shows a typical SHG transient in P in -S out geometry using a ∼100 fs pump at 1400 nm with a fluence of 360 µJ/cm 2 from an optical parametric amplifier. We observe an ultrafast suppression rather than enhancement of the SHG intensity after the pump pulse, which rules out any possibility that the enhanced SHG intensity below T N is due to a laser-induced effect on the timescale of the laser pulse.

Temperature (K)
Time (ps) The data is measured at the φ that gives maximum SH intensity for P in-S out. Solid line is a guide to the eye. Inset shows a schematic of the experimental geometry with a normally incident pump and oblique probe beam.
Using Eq. (1), we achieve good fits to the P in -S out RA patterns, as shown in Fig. 1(c) inset and Fig. 4(b). Fits to the other polarization geometries are shown in Fig. S4. In general, the components of the ED nonlinear susceptibility tensor in Eq. (2) will depend on the surface magnetization, M s , as shown explicitly in the following section.
IV. Fitting the high-field RA data In this section, we detail how the magnetization influences the second-harmonic response in Sr 2 IrO 4 , following the general framework laid out by Ref. [8] and Ref. [9]. Above the metamagnetic transition (∼200 mT), Sr 2 IrO 4 transitions from the − + +− ground state (magnetic point group 2/m1 ) into the ferromagnetic + + ++ state (magnetic point group 2 /m ) [10]. In the + + ++ state, there are three SHG contributions: bulk EQ and surface ED as described above in Section III, and a new magnetization-dependent MD process We note that the MD process is the second-harmonic nonlinear magnetization of the sample, but our data does not allow us to draw a distinction between this process and that of Eq. (7).
The total SHG intensity in the + + ++ state is given by where M (M s ) is the static bulk (surface) magnetization, A MD is a constant, and klm is the Levi-Civita symbol. The magnetization-dependent nonlinear susceptibility tensors may be expanded as χ ijk (M) = χ 0,ijk + χ ijkl M l , where we omit terms that are nonlinear in the magnetization. χ 0,ijk represents the magnetization-independent crystallographic contribution, which is a third-rank axial (polar) tensor for the MD (surface ED) process [6]. These processes give φ-independent terms in the SHG intensity, which are neglected because they cannot explain the high-temperature SHG patterns and they are not necessary for the fits [6]. In the MD (surface ED) case, we write χ ijkl as χ s,ijkl ), which is a fourth-rank polar (axial) tensor that respects the high-temperature crystallographic point group. Due to broken a-b symmetry, we must assume either orthorhombic or monoclinic symmetry for χ MD(i) ijkl and χ ED(i) s,ijkl . Under monoclinic symmetry [7], In our experiments, we apply an in-plane magnetic field that creates a static in-plane magnetization M = (M x , M y , 0), and leads to the third-rank axial c-tensor The MD contributions to the SH electric fields are then given by E MD ss (2ω) ∼ cos θ(χ yxy cos 3 φ + (χ xxy + χ yxx − χ yyy ) cos 2 φ sin φ + (χ xxx − χ xyy − χ yyx ) cos φ sin 2 φ − χ xyx sin 3 φ).
In the expressions above, χ ijk = M x χ ijkx + M y χ ijky . A similar exercise can be performed for the surface ED case, giving the third-rank polar c-tensor These tensor elements can be inserted into Eq. (3)-(6) to get the full surface-magnetization-dependent ED secondharmonic electric field expressions.
In Fig. 2(a) and Fig. S3, we display RA-SHG plots under an applied in-plane magnetic field of 0.37 T with varying field direction, β. Using Eq. (8), we achieve good simultaneous fits (as shown by the solid black lines) for each polarization geometry by only varying β and assuming constant bulk and surface magnetization magnitudes, M 0 and M s,0 , that both rotate with the field direction, M x = M 0 cos β, M y = M 0 sin β, M s,x = M s,0 cos β, and M s,y = M s,0 sin β. The clear pattern differences for β = 0°and β = 90°qualitatively show that tetragonal symmetry is broken. One must assume symmetry reductions to monoclinic or orthorhombic point groups to fit the data. V. β dependence of the anomalous C1 term In Fig. 2(b), we show the summation of the field-dependent RA patterns over all β, which isolates β-independent symmetries. We also show the results of two different fitting schemes, in which the anomalous C 1 term either rotates with β or is independent of β. In the first scheme, labeled "Free C 1 ", we use Eq. (8), which includes a fixed crystallographic EQ contribution, a β-dependent MD contribution, and an anomalous β-dependent C 1 ED contribution. As shown by Fig. 2(a) and Fig. S3, the fits are excellent, and as a result, the summation of the fielddependent fits matches the data well, as shown by the black curve in Fig. 2(b). In the second scheme, labeled "Fixed C 1 ", we first determine the β-independent EQ and anomalous C 1 ED tensor elements by fitting the zero-field RA pattern in Fig. 1(c) inset with Eq. (1). We then fit the 370 mT data by coherently adding the β-dependent MD term to the fixed EQ and ED terms. As shown by the red curve in Fig. 2(b), this fitting scheme imparts a noticeable C 1 dependence in the summation over β, which does not fit the data well. We can therefore conclude that the anomalous C 1 term rotates with β.
VI. Fitting RA data by scaling the MD contribution As discussed in the main text, another possible explanation we examined for the field-dependent SHG data is that the stacking order is altered to be + + ++ near the sample surface. In this scenario, we expect that all RA patterns can be described by only two terms: 4/m EQ and a magnetization-dependent 2 /m MD SHG. Furthermore, there should be a perfect mapping from the high-field to the zero-field case by simply changing the MD tensor by a complex scale factor. To test this, we first fit the zero-(high-) field data by itself, then attempted to fit the high-(zero-) field data only by scaling the MD contribution to account for the different magnetization. The results of this exercise are shown by the solid black lines in Fig. S4. While the zero-and high-field patterns are qualitatively similar, there are clear discrepancies between the black curves and the data, which is further evidence against the existence of an anomalous ferromagnetic stacking in part of the sample. On the other hand, the data fit well to the model in Section III and IV, as shown by the shaded areas. VII. SHG evidence for a-b symmetry breaking and twinning As we discussed in the main text and Section V, we can expose β-independent symmetries by summing the RA-SHG patterns over magnetic field angles β from 0 to 2π. A clear twofold symmetry is observed in Fig. 2(b), signifying a different magnetic field response along a and b. This symmetry persists at 220 K [ Fig. S5(a)], though with a slightly more subtle distortion due to the weakened magnetic term. We do not observe appreciable twofold distortions in the EQ SHG response at room temperature [ Fig. S5(b)]. There are two potential explanations for these data. The first is pseudospin-lattice coupling, which creates a uniaxial anisotropy in-plane along the direction of the layer magnetization. The effect of pseudospin-lattice coupling shows itself in the intermediate magnetic field range [11,12], where the applied field is comparable to the critical field of the metamagnetic transition, as seen, for example, in the twofold symmetry in the magnetic field direction dependence of out-of-plane anisotropic magnetoresistance [13,14]. This metamagnetic transition and anisotropy are strongly dependent on temperature [11]. At 80 K, the magnetization becomes fully saturated for all field directions above ∼0.3 T, while at 220 K, the saturation is achieved by ∼0.1 T. Above these fields, one expects a return to C 4 symmetry in the field direction dependence, which should be the case for our SHG experiments. As a result, the observed twofold symmetry points instead to an alternative explanation: the existence of a distortion below tetragonal symmetry that exists above T N . This is consistent with the observation in Fig. 4c that the order parameter only changes by 0°or 180°upon thermal cycling in zero field. However, the distortion is subtle enough that it is not detected by RA-SHG above T N . Figure S6 shows the magnetic field response from two different SHG domains with C 1 order parameters that are rotated by 90°at zero field. In the first region (domain A) with β = 0°(top left), the P in -S out pattern shows two lobes of equal intensity. When the field is rotated to β = 90°, however, we only observe a single strong lobe (bottom left). This is the same behavior that was found in Fig. 2(a). In contrast, for domain B, the field-direction-dependent response switches; one strong lobe is seen for β = 0°, while two are present for β = 90°. By inspecting the patterns, we see that region B is related to A by an interchange of the a and b axes. In other words, the SHG domains related by relative surface magnetizations of 90°are a-b twins, with locked magnetic and crystallographic ordering.