Charge order textures induced by non-linear lattice coupling in a half-doped manganite

The self-organization of strongly interacting electrons into superlattice structures underlies the properties of many quantum materials. How these electrons arrange within the superlattice dictates what symmetries are broken and what ground states are stabilized. Here we show that cryogenic scanning transmission electron microscopy enables direct mapping of local symmetries and order at the intra-unit-cell level in the model charge-ordered system Nd$_{1/2}$Sr$_{1/2}$MnO$_{3}$. In addition to imaging the prototypical site-centered charge order, we discover the nanoscale coexistence of an exotic intermediate state which mixes site and bond order and breaks inversion symmetry. We further show that nonlinear coupling of distinct lattice modes controls the selection between competing ground states. The results demonstrate the importance of lattice coupling for understanding and manipulating the character of electronic self-organization and highlight a novel method for probing local order in a broad range of strongly correlated systems.

The self-organization of strongly interacting electrons into superlattice structures underlies the properties of many quantum materials. How these electrons arrange within the superlattice dictates what symmetries are broken and what ground states are stabilized. Here we show that cryogenic scanning transmission electron microscopy enables direct mapping of local symmetries and order at the intra-unit-cell level in the model charge-ordered system Nd 1/2 Sr 1/2 MnO 3 . In addition to imaging the prototypical site-centered charge order, we discover the nanoscale coexistence of an exotic intermediate state which mixes site and bond order and breaks inversion symmetry. We further show that nonlinear coupling of distinct lattice modes controls the selection between competing ground states.
The results demonstrate the importance of lattice coupling for understanding and manipulating the character of electronic self-organization and highlight a novel method for probing local order in a broad range of strongly correlated systems.
Strong interactions between electrons and the atomic lattice often lead to their selforganization into ordered spatial patterns [1][2][3][4]. One well-known example is charge ordering, the spatial modulation of the electronic charge density which forms superlattices and governs the properties of many exotic materials, from oxides to transition-metal chalcogenides to charge-transfer salts [5][6][7]. In general, charge ordering is studied at scales larger than the superlattice, with a focus on the average periodicity of the charge modulation or the degree of long range order [8,9]. On the other hand, the microscopic arrangement at sub-unit-cell length scales, such as whether the electrons reside on the atomic sites or bonds ( Fig. 1a and   b), is significantly more challenging to measure, but dictates the symmetry of the system in addition to the mechanism underlying electronic order [10][11][12].
Perovskite 3d transition-metal oxides are a class of materials in which charge order plays an especially important role, influencing antiferromagnetic order, high-temperature superconductivity, and colossal magnetoresistance [5,13,14]. The strong hybridization between the transition metal (site) and oxygen (bond) in these systems contributes to their rich, intricate electronic structure landscape [15]. In cuprates, for instance, recent bulk measurements suggest that the charge-ordered phase may contain a site-centered modulation in addition to the putative bond-centered d-wave modulation [10,16]. A similarly complex and fundamental debate concerns the half-doped manganites, an ideal playground for exploring coupled charge and orbital orders and phase competition [17]. The widely accepted model, proposed over 60 years ago, invokes a site-centered modulation in which electrons localize on manganese sites and form a zig-zag orbital order pattern which doubles the unit cell ( Fig. 1c) [18][19][20]. An alternative bond-centered state has also been proposed [21]; in this case a charge modulation does not occur so the superlattice is generated from orbital ordering (Fig. 1d). Subsequent experimental studies, however, have not conclusively confirmed this scenario [20,22,23]. Theory further predicts an even more exotic intermediate state which mixes both site and bond characters, leading to the breaking of inversion symmetry through the formation of uncompensated electric dipoles [24,25]. Determining if such a state exists or whether charge order corresponds to a pure site-or bond-centered state remains a fundamental challenge, with broad and urgent implications for other classes of materials.
Manganite compounds exhibit large interactions between the lattice and the electronic degrees of freedom, and so most experimental proposals for site-or bond-centered charge order models have relied on obtaining the average lattice distortions and crystal symmetry [19,21]. Here we employ a novel methodology for probing charge order, based on atomicresolution cryogenic scanning transmission electron microscopy (cryo-STEM) and the direct measurement of picoscale lattice displacements in real space. By mapping the displacement patterns within the charge order superlattice, this method enables us to visualize distinct, coexisting arrangements of charge order and, importantly, establish the existence of the intermediate phase. Absent the real space visualization, this phase would remain hidden, partly due to the challenge of disentangling nanoscale domains of pure bond/site order from genuine mixed states in bulk-averaged measurements. Going further, cryogenic STEM imaging and theoretical calculations reveal coupled secondary order parameters that govern the nature of the charge order phase. These atomic-scale visualizations not only advance our understanding of charge ordering through a novel lens but also reveal a tuning knob for manipulating electronic self-organization via lattice coupling.
The material we focus on is a Nd 1/2 Sr 1/2 MnO 3 thin film which was grown using pulsed laser deposition on a (110)-oriented SrTiO 3 substrate. The particular substrate orientation, which imparts in-plane anisotropy on the film, reproduces the electronic and structural transitions found in the bulk counterpart including the charge order transition [26][27][28]. The temperature dependence of the magnetization in this film matches that of previous reported epitaxial films. Figure 1e shows a projection HAADF-STEM image of the crystal along the [010] orientation (in the P nma space group setting) below the charge order transition temperature. The atomically-resolved Nd/Sr columns appear brighter than the Mn columns because the contrast in HAADF scales strongly with the atomic number. Near the crystalline peaks in the Fourier transform amplitude (Fig. 1f), superlattice peaks (blue arrows) appear at low temperature, indicating the formation of a modulated structure. These are located at Q CO = (1/2, 0, 0) as expected for charge ordering at half doping. The modulation is unidirectional, however, some regions of the sample exhibit bi-directional modulations.
These might arise from the coexistence of small charge order domains (either in-plane or out-of-plane) within a single crystal twin or from the presence of crystalline twins that establish the direction of the charge order wavevector. By mapping both the charge order and the crystalline order parameters, we find that the orthogonal charge order domains are coupled to crystalline twins in the sample ( Supplementary Fig. 3).
The exact intra-unit-cell arrangement within charge order superlattices is key to understanding their microscopic origins and interactions with other electronic phases. In the The size of the circle indicates the electronic density on the site. b, Charge order modulation centered on the bonds. Ellipses represent stronger, more electron-rich bonds. The sites remain equivalent in density. c, Prototypical site-centered charge and orbital order model in the halfdoped manganites, including orbital order. d, Alternative bond-centered model in which all Mn sites have the same charge configuration (Mn 3.5+ ). In this case, the supercell is generated by the orbital order. e, High-angle annular dark-field (HAADF) scanning transmission electron microscopy (STEM) projection image at 93 K, where we expect the system to be in the charge-ordered phase.
Bright atomic columns represent Nd/Sr sites and dark atomic columns represent Mn sites. The scale bar corresponds to 2 nm. f, Fourier transform amplitude of the STEM image. At 93 K, superlattice peaks (arrows) with Q CO = (1/2, 0, 0) are evident, indicating the formation of a twofold superlattice in real space (box in E). majority of theoretical treatments of charge-ordered phases, the ground state is discussed in terms of pure electronic degrees of freedom. Site-centered charge order in half-doped manganites, for instance, is described as the alternation of Mn 4+ and Mn 3+ species [18,19]; however, the degree of charge disproportionation is much smaller and is better described by Mn valences of 3.5 + δ (Mn 4+ -like) and 3.5 − δ (Mn 3+ -like) with δ 0.5 [20,29]. In reality, the crystal structure also undergoes a variety of complex atomic displacements, such as Jahn-Teller and breathing distortions, which alter the bonding network and hence the electronic configuration. Therefore, emergent charge and orbital textures are closely linked to the pattern and symmetry of said displacements [30,31]. For site-centered order, bulk X-ray or neutron structural data suggests that the crystal adopts P 2 1 /m space group symmetry with bond distortions consistent with charge localization on the Mn sites [19,20,22].
The experimental report for the bond-centered model, on the other hand, found a different space group symmetry (P nm2 1 ) and a distinct displacement pattern associated with it [21].
The ability to probe intra-unit-cell lattice distortions can therefore determine or even reveal novel ground states.
To firmly connect patterns of atomic displacements to the reported crystal symmetries and hence to the models of electronic order, we first explore the possible distorted structures that emerge from the high-symmetry phase (space group P nma) using group theory. Any distorted structure must double the unit cell, therefore we require atomic displacements whose wavevector is q = (1/2, 0, 0). The relevant distortion consistent with this requirement is the X1 displacement mode, or irreducible representation (irrep), which affects two inequivalent Mn sublattices in the high-symmetry unit cell (See Methods section 4). In other words, the X1 irrep is two-dimensional with the first dimension corresponding to the first Mn sublattice (Mn-1 in blue) and the second dimension to the second Mn sublattice (Mn-2 in red), as shown in Fig. 2a. If X1 displacements, which consist of a complex set of atomic distortions including the transverse Mn displacements shown in Fig. 2, affect only the first sublattice, the resulting crystal structure has P 2 1 /m symmetry, which is consistent with site-centered order (Fig. 2b,c). If they occur in both sublattices and with equal magnitude (Fig. 2b,d), the resulting crystal structure has P nm2 1 symmetry, which matches that of bond-centered order. A third structure with P m symmetry can be obtained by having displacements on both sublattices but with different magnitude (Fig. 2b,  We now pinpoint the underlying ordering model in Nd 1/2 Sr 1/2 MnO 3 , by mapping the picoscale lattice degrees of freedom using HAADF-STEM at low temperature. A reference image lacking the periodic modulation is generated by removing the contribution of the Q CO superlattice peaks in the Fourier transform [32,33]. The displacements are then extracted by mapping the atomic positions in the original image and the reference image. Figure 3a shows a STEM image overlaid with arrows indicating Mn column displacements. The area of the triangle denotes the magnitude of displacement and the color represents the angle of the displacement relative to the wavevector. The dominant displacements have transverse polarization and generate a twofold superlattice. The displacement pattern shown in Fig. 3a matches site-centered order. The first Mn sublattice has a large displacement amplitude (8.2(9) pm) while the second has a comparatively negligible amplitude (1.3(6) pm), similar to the group theory prediction of site-centered order. The manganese displacements as a function of position, r, are sinusoidal ( Fig. 3c) and may be described by ∆ Mn (r) = A sin(Q CO .r + φ) where ∆ Mn (r), A, Q CO , and φ are the displacement, amplitude and wavevectors, and the phase, respectively. The phase determines the centering of the modulation relative to the Mn sites, with site-and bond-centered order corresponding to φ = nπ (n is an integer) and φ = nπ/2 (n is an odd integer), respectively.
A deviation of the phase from these two limits corresponds to a state intermediate between site and bond order. Fitting the displacements in the region shown in Fig. 3a further supports that the modulation is predominantly site-centered with φ = 0.04(1) π. From this measurement, we thus confirm the presence of the prototypical site-centered state.
Remarkably, we discover that within the same sample another region has a distinct displacement pattern: prominent Mn displacements occur in both Mn sublattices (Fig. 3b).
Importantly, the displacements have different amplitudes on each sublattice, unlike the predicted pattern for purely bond-centered order. By comparing to the group theory analysis of displacement patterns, the observed structure is consistent with a state which is intermediate between pure site and bond order (P m structure). Figure   An intriguing consequence of the overlap of site-and bond-centered order is that additional crystal symmetries may be broken [24]. In oxides or charge-transfer salts, intermediate order is predicted to break inversion symmetry due to the formation of uncom- Another key insight from these local visualizations is that ground states with distinct symmetries may coexist within the same system. Both Figs. 3A and B are taken from the same field of view image (Fig. 3g), in which the site-centered phase (left side) transforms into the intermediate phase (right side) over a few unit cells. This suggests that the energies of site and intermediate orders are comparable and likely linked to subtle spatial fluctuations inherent to strongly interacting oxide materials [27]. Such nanoscale coexistence further highlights the long-standing difficulty in determining the correct charge order model from bulk-averaged measurements, a situation exacerbated by the presence of crystalline twins (see for example Supplementary Fig. 3) and other inhomogeneity.
Informed by the atomic-scale evidence for site and intermediate orders, we next examine the possible origin of these two states. While we have focused on the X1 displacements so far, the low temperature phase may also contain additional structural responses that are allowed by symmetry (Fig. 4a). In particular, a Γ + 4 mode, which appears as an oxygen breathing distortion, and a Γ − 2 mode, which appears as antipolar displacements of Mn atoms, are such symmetry-allowed responses (see Methods section 4). To understand their role we expand the Landau free energy about the high-symmetry phase up to fourth order in powers of the Γ + 4 , Γ − 2 , and X1 displacement mode amplitudes (see equation (3) in Methods section 5 for the expansion). Significantly, this expansion contains two relevant third-order coupling terms: Our study shows that cryogenic STEM is an emerging methodology for probing complex electronic ordering phenomena through the important but often neglected lens of the lattice degrees of freedom. Using this approach, we discovered the co-existence of distinct charge order configurations including a pure site-centered phase and a more exotic intermediate phase which breaks inversion symmetry. We also revealed that the stability of these different configurations depends on unique lattice couplings. Such a rich electronic and structural landscape should be relevant to other charge-ordered systems, including cuprates and nickelates where the exact microscopic arrangement associated with various electronic instabilities remains under intense scrutiny. Similar to the current study, novel insights may be achieved through direct, real space visualizations that can spatially disentangle different ordering models and characterize the intra-unit-cell structure and symmetry in detail.

COMPETING INTERESTS
The authors declare that they have no competing financial interests.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author on reasonable request.

CORRESPONDENCE
Correspondence and requests for materials should be addressed to L.F.K.

Magnetization
The magnetization data were measured using a superconducting quantum interference device (SQUID) magnetometer with an in-plane magnetic field. The sample was cooled down from room temperature to 10 K in zero external field, and its magnetization was measured in a 1000 Oe field while warming.

Scanning transmission electron microscopy & data analysis
Experimental details. Thin electron transparent samples were fabricated on a FEI Strata 400 focused ion beam using the standard liftout method with the stage tilted by 45 • .
The stage tilt enabled preparation of TEM specimens in the correct crystal orientation, as the films are grown on (110)-oriented SrTiO 3 substrates. High-angle annular dark-field (HAADF) STEM imaging was performed in an aberration-corrected FEI Titan Themis microscope operating at 300 kV and 50 pA probe current. The convergence semi-angle was 21.4 mrad and the inner collection angle was 68 mrad. Cryogenic STEM measurements were performed using a Gatan double-tilt side-entry liquid nitrogen TEM holder (Model 636) with a base temperature of 93 K. We acquired multiple fast-acquisition images (0.5 or 1 µs per pixel and 1024×1024 pixels per image) to overcome sample drift which occurs at low temperature, and aligned the image series using a rigid registration method optimized for noisy images [2].
Characterization of epitaxial film quality and fluctuations. Room temperature data.
The 80 nm epitaxial Nd 1/2 Sr 1/2 MnO 3 /SrTiO 3 thin film is characterized using HAADF-STEM. Supplementary Fig. 1a shows a low magnification image that highlights the epitaxial growth of the film, with the substrate appearing dimmer than the film due to the atomic number contrast in HAADF-STEM. At room temperature, no superlattice peaks associated with charge ordering are observed in the Fourier transform of STEM data ( Supplementary   Fig. 2), in agreement with the magnetization data which shows that charge order occurs below 150 K. We note the presence of broad superlattice peaks ( Supplementary Fig.2b, green circles) which occur along the pseudocubic direction at Q cat = ( 1 2 , 0, 0) pc . These peaks exist above and below the charge ordering transition temperature and are associated with local short-range A-site (Nd/Sr) cation ordering.
Crystal twins. Structural inhomogeneities in bulk crystals, such as crystal twins and chemical fluctuations, often impede, or at the very least complicate, diffraction experiments since they require complex models or assumptions for structural refinements [3]. Real space imaging, on the other hand, allows us to isolate and characterize local structural or chemical fluctuations in the manganite, from the occurrence of crystal and charge order domains to short-range cation ordering. Supplementary Figure 3a shows a 40 nm region of the film in which we observe superlattice peaks along two orthogonal directions, Q CO 1 and Q CO 2 , in the Fourier transform amplitude (Supplementary Fig. 3e). In the case of bi-directional order, it is important to determine whether the modulations are overlapping or spatially separated [4,5]. By Fourier filtering charge order superlattice peaks near the (200) and (002) Bragg reflections, we obtain the local amplitude of the modulation along each direction, A α (r) with α = 1, 2. We define a local charge-order unidirectionality parameter which determines whether the ordering pattern in locally unidirectional (positive for the Q CO 2 direction, negative for the Q CO 1 direction and close to zero for overlapped). The two charge order domains are separated by a well-defined boundary ( Supplementary Fig. 3d), suggestive of a crystal twin effect rather than anomalous competition between unidirectional order parameters as seen in a previous study on a related system [5]. Indeed, we observe a splitting of crystalline Bragg peaks which indicates the presence of orthorhombic twins ( Supplementary Fig. 3c). A map of the local crystal unidirectionality parameter, Σ Crystal (r), determined using the split (402) Bragg peaks shows that it correlates with the charge order domains (Supplementary Fig. 3e). The crystal twinning, and hence the charge order twinning, is likely due to the relaxation of the epitaxial film since it occurs along the growth direction.
Mapping primary displacements. The method for extracting periodic lattice displacements and rigorous numerical tests are described in great detail in previous work [5,6]. To briefly summarize the method, we canvas the low temperature Fourier transform amplitude for all Q CO = (1/2, 0, 0) superlattice peaks associated with charge ordering and damp their amplitude to the background level. We then apply an inverse Fourier transform, generating a reference image in which the modulation is removed. The positions of the atomic columns are determined by fitting 2D Gaussian functions to STEM images [5][6][7] (Supplementary Fig. 7). Incoherent stacking of the disordered displacements cannot a yield well-defined, periodic displacement pattern in projection. Therefore, the X1 displacements in the intermediate phase are intrinsic and coherent along the imaging direction and do not reflect a stacking of phase-shifted domains of site-centered displacements.
Mapping secondary displacements. The Γ − 2 mode is a zone-center mode of the orthorhombic P nma cell and can be extracted using the same method if applied to the relevant peaks [5,6]. To determine these peaks, we note that the Γ − 2 modes (antipolar displacements of Mn atoms) correspond to zone-boundary modes of the cubic cell ( Supplementary Fig. 6, dashed red square), therefore they have well defined and finite q peaks in the Fourier transform amplitude (Supplementary Fig. 6a, pink circles). Indeed, these modes' periodicity is encoded in the q = (1, 0, 0)-type Fourier peaks ( Supplementary Fig. 6, pink circles). To extract these displacements, we generate a reference image by removing the contribution of these particular peaks.
The obtained displacements are shown in Supplementary Fig. 6c (full field of view) and e (zoom-in). The projection of the displacements along the Γ − 2 distortion mode (see Methods section 4 for details about this mode and Supplementary Fig. 6d for the displacement pattern) is shown in Supplementary Fig. 6f.
Over the full field of view, the amplitude and coherence of this mode fluctuates in lockstep with the fluctuation in X1 displacements (Supplementary Fig. 7). In the site-centered region, Γ − 2 displacements are disordered and weak ( Supplementary Fig. 7c and d). However, they are strong and coherent in the intermediate region ( Supplementary Fig. 7e and f)  A purely odd signal will have a purely imaginary and odd F(q) while a purely even signal will have a purely real and even F(q). An experimental image of a crystal structure, however, is much more complicated and has contributions from imperfect centering of the signal and phase offsets, which means that the Fourier transforms will have both real and imaginary components. Since both charge-ordered structures occur within the same image, we can focus on the prominent differences in Im{F(q)} between the two models. After taking the local Fourier transforms, we extract integrated line cuts across the Bragg peaks (pseudo-cubic directions). The width of the integration window corresponds to the width of the Bragg peaks.

Landau Theory
Subgroups. Group theory analysis and symmetry decomposition were performed using the ISOTROPY software suite [12,13]. First, we enumerate the possible distorted structures, or subgroups, that emerge from the high-symmetry structure in the P nma phase. To double the unit cell along the [100] direction, we require a zone-boundary displacement mode with q = ( 1 2 , 0, 0). There are only two such order parameters, with one transforming like the irreducible representation (irrep) X1 and the other like X2 (notation of Miller and Love [12,14]). The subgroups that are generated by distinct directions of the X1 irrep are P 2 1 /m, P nm2 1 and P m while those that can be generated by X2 are P 2 1 /c, P na2 1 and P c.
Displacement patterns. We do not observe any X2-type displacements in the STEM data and no previous structural refinements have adequately described the charge-ordered structure in terms of the X2 displacement modes or their associated subgroups [15,16].
Thus, the X1 mode is the relevant order parameter. While this mode involves displacements of both cations and oxygen atoms, we will focus on the Mn atoms (sites). For the Mn atoms, the X1 mode appears as periodic lattice displacements which can either be longitudinal, transverse and in-plane, or transverse and out-of-the plane. Our data, as well as previous refinements, show that the transverse in-plane component dominates; therefore we represent the X1 mode in terms of the transverse in-plane component, as done in Figs. 1 and 2 in the main text.
As discussed in the main text, the X1 order parameter is two-dimensional with each dimension corresponding to a Mn sublattice in the P nma cell. The P 2 1 /m subgroup is generated by an X1 distortion along a single dimension (s 1 ,s 2 =0); P nm2 1 by X1 distortions along both dimensions and with equal magnitude (s 1 ,s 2 =s 1 ); and P m by X1 distortions along both dimensions but with different magnitude (s 1 ,s 2 = s 1 ).
Energy expansion and secondary modes. We performed the free energy expansion about the high-symmetry phase (space group P nma) in terms of polynomials of the X1 order parameter along the general direction, (s 1 ,s 2 ), up to fourth order [17]. Importantly, we also include the secondary modes that are allowed by symmetry [18,19]. The possible secondary order parameters are limited to Γ modes with q = (0, 0, 0) with respect to the P nma cell, and include Γ + 1 , Γ + 4 , Γ − 2 , and Γ − 3 . The Γ + 1 mode is a strain or volume change which appears, by symmetry, in any phase transition. The Γ + 4 mode is a displacement mode appearing as breathing distortion of the oxygen octahedra, which, from a chemical bond perspective, corresponds to charge localization [20,21]. Note that this "charge localization" mode alone does not double the unit cell so it must require coupling to the X1 mode. Further, according to DFT calculations, this mode by itself raises the energy of the system and only acquires a finite but small amplitude by coupling to other modes (see Methods section 5 and ref. [21]). The Γ − 2 mode involves predominantly oxygen displacements as well as antipolar displacements of the metal cations. Since only the cations are visible in HAADF-STEM we will focus on the Mn displacement pattern associated with the Γ − 2 irrep. Finally, Γ − 3 corresponds to a polar displacement mode.
The free energy will contain the usual quadratic and quartic terms for each of these listed order parameters (X1 and all the listed Γ modes) and, more importantly, coupling terms whose form depends on the symmetry of the order parameters. The full coupling energy is given by where Q V , Q B , Q AP , Q P , and (s 1 , s 2 ) are the amplitude of the Γ 1 +, Γ + 4 , Γ − 2 , Γ − 3 , and X1 distortions, respectively; and δ vss , δ bss , δ tss δ btp are the coupling coefficients. The coefficients can be determined using DFT calculations (see Methods section 5 and Supplementary Table III). We find that δ bss , δ tss and δ btp are negative, which means that the corresponding coupling terms lower the energy.
Interpretation of the coupling terms. The first coupling term, Q V (s 2 1 + s 2 2 ), reflects the coupling of strain and volume change to the phase transition. Moreover, the Γ + 1 mode maintains the symmetry of the parent P nma structure and does not favor a particular direction in the X1 order parameter. Its main effect is quantitative, in that it slightly alters the coefficients of the usual s 1 /s 2 quadratic terms in the Landau energy expansion.
The interplay between the second and third coupling terms, on the other hand, is important for determining the nature of the ground state, as discussed in the main text. The second term induces the Γ + 4 (Q B = 0) distortion and a (s 1 , 0) amplitude in the X1 mode, which yields site-centered order with P 2 1 /m symmetry. The third term induces the Γ − 2 (Q AP = 0) distortion and a (s 1 , s 1 ) X1 amplitude, which yields bond-centered order with P nm2 1 symmetry. Only when both terms occur do we obtain the intermediate phase with Finally, we note that in the main text we mainly focused on the second and third terms in F 3 because they are directly related to determining the stability of site and intermediate orders. That said, the other coupling terms should be considered if a more quantitative analysis is required, especially in future studies of electric field or strain effects on charge order.

Density Functional Theory
Computational details. We perform density functional theory (DFT) + U calculations as implemented in the Vienna ab initio simulation package (VASP) [25,26] with the PBEsol exchange-correlation functional [27] and a 500 eV energy cutoff for the plane-wave basis. We  Table I). Note that in the DFT-calculated structures, the computational requirement for A-site (Nd/Sr) ordering changes the symmetry.
Displacement modes. To assess the symmetries and amplitudes of the structural distortions that contribute to the P 2 1 ma structure, we report in Supplementary Table II a struc-tural decomposition of the DFT-relaxed P 2 1 ma structure into symmetry adapted modes of the P 4/mmm reference structure [12,13]. Note that since we are using a P 4/mmm (rather than P nma) reference structure here, the labelling of the irreps is different, although they represent the same physical atomic displacements. The M − 5 mode, which is an octahedral tilting mode (contained in P nma) together with antipolar Mn displacements, and the Σ 2 mode, which corresponds to the X1-type periodic lattice displacements, have the largest amplitudes, followed by the M + 4 oxygen breathing distortion. The polar (Γ − 5 ) distortion is also present. For each distortion, we use the displacement of each distinct atomic type to construct the relevant symmetry-adapted mode vector, which we use for our energy surface calculations below [21].
Landau expansion about the P 4/mmm reference structure. With the P 4/mmm irreps defined above, we perform the energy expansion of the free energy, F 1 , about the high-symmetry P 4/mmm reference structure [19]. To fourth order, the energy is where F 3 is the third-order coupling energy Here, Q AP , Q B , Q P , and (s 1 ,s 2 ) are the amplitudes of the M − 5 , M + 4 , Γ − 5 and Σ 2 displacement modes, respectively (see Table II in Methods section for the amplitudes). Note that the coupling terms in F 3 have the same functional form as in equation (2) in Methods section 4.
Determining the coupling coefficients. We calculate energy surfaces by "freezing in" increasing amplitudes of each mode and calculating the total energy with DFT, as shown in Supplementary Fig. 8. By fitting these energy surfaces we extract the Landau expansion coefficients α i and β i in equation (3). Note that we freeze in the two-dimensional Σ 2 order parameter along the direction (0.907, 0.421), which is the direction this order parameter takes in the DFT-relaxed P 2 1 ma structure with A-AFM order. That is, we freeze in amplitude Q s = s 2 1 + s 2 2 where s 1 = 0.907Q s and s 2 = 0.421Q s . For this reason, we omit the γ s coefficient from equation (3).
To calculate the coupling coefficients δ tss and δ bss in equation (4), we fix the Σ 2 distortion amplitude to the value obtained in the DFT-relaxed structure (Q s = 0.68Å), freeze in the M − 5 distortion (M + 4 breathing distortion), and then fit the resulting energy surface ( Supplementary Fig. 9a,b). To calculate the trilinear δ tbp coefficient, we fix the M − 5 and M + 4 amplitudes to the values in the DFT-relaxed structure and freeze in the polar mode ( Supplementary Fig. 9c). The resulting coupling coefficients are reported in Supplementary   Table III. The α s , α b , α p coefficients are all positive which means that the X1-like displacements, the breathing distortion, and the polar distortion by themselves do not lower the energy. The breathing distortion and the polar distortion significantly raise the energy so we limited the energy expansion to the second powers of Q B and Q P . Meanwhile, the coupling coefficients δ tss , δ bss and δ btp are negative, which means that the coupling terms lower the energy of the system by non-linearly inducing distortions that by themselves are energetically unfavorable.
The energy surfaces and line plots in Supplementary Fig. 10 and Fig. 4 of the main text are     Supplementary Fig. 9. Energy surfaces used to calculate the coupling coefficients in equation (4).
a,b, The Σ 2 (X1-like) distortion is set to a fixed amplitude and then Q AP , Q B is frozen in. c, Both Q AP and Q B are fixed, and the polar mode Q P is frozen in.  c, By including only the δ tss Q AP s 1 s 2 coupling term, the energy is minimized at the (s 1 ,s 2 = s 1 ) coordinate of the X1-like order parameter (bond-centered order). d, When both coupling terms are present, the energy is minimized at the (s 1 , s 2 = s 1 ) coordinate of the X1-like order parameter (intermediate order). We set the amplitudes of Q B , Q AP , and Q P to the values in Supplementary  Fig. 4 in the main text were taken from panels b and d, along the (0,s 2 ) directions and passing though the minimum energy positions (crosses).