Local orthorhombic lattice distortions in the paramagnetic tetragonal phase of superconducting NaFe1−xNixAs

Understanding the interplay between nematicity, magnetism and superconductivity is pivotal for elucidating the physics of iron-based superconductors. Here we use neutron scattering to probe magnetic and nematic orders throughout the phase diagram of NaFe1−xNixAs, finding that while both static antiferromagnetic and nematic orders compete with superconductivity, the onset temperatures for these two orders remain well separated approaching the putative quantum critical points. We uncover local orthorhombic distortions that persist well above the tetragonal-to-orthorhombic structural transition temperature Ts in underdoped samples and extend well into the overdoped regime that exhibits neither magnetic nor structural phase transitions. These unexpected local orthorhombic distortions display Curie–Weiss temperature dependence and become suppressed below the superconducting transition temperature Tc, suggesting that they result from the large nematic susceptibility near optimal superconductivity. Our results account for observations of rotational symmetry breaking above Ts, and attest to the presence of significant nematic fluctuations near optimal superconductivity.

I ron pnictide superconductors are a large class of materials hosting unconventional superconductivity that emerges from antiferromagnetically ordered parent compounds [ Fig. 1a]. Unique to iron pnictides is the tetragonal-to-orthorhombic structural transition at T s , where the underlying lattice changes from exhibiting fourfold (C 4 ) above T s to twofold (C 2 ) rotational symmetry below T s , which occurs either simultaneously with or above the antiferromagnetic (AF) phase transition temperature T N [ Fig. 1b] 1,2 . The large electronic anisotropy present in the paramagnetic orthorhombic phase has been ascribed to an electronic nematic state 3-5 that couples with the shear strain of the lattice, the orthorhombicity δ [=(a − b)/(a + b), where a and b are in-plane orthorhombic lattice parameters], therefore acts as a proxy for the nematic order parameter. In the paramagnetic tetragonal state, the nematic susceptibility can be measured via determining the resistivity anisotropy induced by anisotropic inplane strain 6 or by measuring the elastic shear modulus 7,8 . By fitting temperature dependence of nematic susceptibility with a Curie-Weiss form, a nematic quantum critical point (QCP) with Weiss temperature T * → 0 has been identified near optimal superconductivity for different iron-based superconductors 6,8 . Theoretically, the proliferation of nematic fluctuations near the nematic QCP can act to enhance Cooper pairing [9][10][11][12] .
Although C 4 → C 2 symmetry breaking is typically associated with the structural transition at T s , there are numerous reports of its observation well above T s and in overdoped compounds [13][14][15][16][17][18][19] . These observations are either reflective of an intrinsic rotational symmetry-broken phase above T s , which can occur in bulk [13][14][15] or on the surface of the sample 16 , or simply a result from a large nematic susceptibility [17][18][19][20] . In the first case, there is a small, but nonzero nematic order parameter throughout the material above T s , although no additional symmetry breaking occurs below T s , despite the sharp increase of the nematic order parameter. For the latter scenario, only local orthorhombic distortions can be present and the system remains tetragonal on average. One way to differentiate the two scenarios is to directly and quantitatively probe the distribution of the interplanar atomic spacings (d-spacings) and its temperature dependence.
Ideally, when the system becomes orthorhombic, two different in-plane d-spacings, corresponding to different in-plane lattice parameters, can be resolved; on the other hand, when there are only local orthorhombic distortions, the d-spacing distribution only broadens, while the average structure remains tetragonal [ Fig. 1c]. However, experimentally, it can be very difficult to distinguish the two scenarios when δ is too small for a splitting to be resolved, then, a broadening is also seen even when the system goes through a tetragonal-to-orthorhombic phase transition. In such cases, it is more instructive to examine the temperature dependence of the experimentally obtained broadening, characterized either by δ or by the width of the d-spacing distribution, Δd/d [ Fig. 1c]. For a phase transition, the broadening should exhibit a clear order parameter-like onset; for local orthorhombic distortions in an average tetragonal structure, the broadening instead tracks the nematic susceptibility, which exhibits a Curie-Weiss temperature dependence 4 [ Fig. 1c]. An additional complication is that the AF order typically becomes spin-glasslike and sometimes incommensurate near the nematic QCP [21][22][23][24][25] , and given the strong magnetoelastic coupling in iron pnictides 5,8 , it is unclear how such changes in AF order affect the nematic order.
In this work, we use high-resolution neutron diffraction and neutron Larmor diffraction to map out the phase diagram of NaFe 1−x Ni x As 26 , focusing on the interplay between magnetic order, nematic order, and superconductivity near optimal superconductivity. Unlike most other iron pnictide systems, we find T N in NaFe 1−x Ni x As to be continuously suppressed toward T N ≈ T c near optimal doping, while the order remains long-range and commensurate. This allows us to demonstrate that T s and T N in NaFe 1−x Ni x As remain well separated near optimal superconductivity, indicating distinct QCPs associated with nematic and AF orders, similar to the quantum criticality in electrondoped Ba 2 Fe 2−x Ni x As 2 27 . Utilizing the high resolution provided by neutron Larmor diffraction 28,29 , we probed the nematic order parameter in underdoped NaFe 1−x Ni x As below T s and surprisingly, uncovered local orthorhombic distortions well above T s and in overdoped samples without a structural phase transition. Although the average structure is tetragonal in these regimes, broadening of the d-spacing distribution is unambiguously observed. Such local orthorhombic distortions were hinted in previous high-resolution neutron powder diffraction measurements on electron-overdoped NaFe 0.975 Co 0.025 As, where a small broadening of Bragg peaks at low temperature was observed 26 . Regardless of whether orthorhombic distortions are long-range due to a structural phase transition or local in nature, resulting from large nematic susceptibility, we find that they become suppressed inside the superconducting state, similar to AF order. Our results, therefore, elucidate the interplay between AF order, nematicity, and superconductivity in NaFe 1−x Ni x As; at the same time, our observation of local orthorhombic distortions with a Curie-Weiss temperature dependence across the phase diagram accounts for rotational symmetry breaking seen in nominally tetragonal iron pnictides. In addition, our measurements demonstrate that neutron Larmor diffraction can be used to determine the nematic susceptibility of free-standing iron pnictides without the need to apply external stress or strain. These results should stimulate future high-resolution neutron/X-ray diffraction work to study orthorhombic lattice distortion and its temperature dependence in the nominally tetragonal phase of iron-based superconductors.

Results
Overall phase diagram of NaFe 1−x Ni x As. Our results are reported using the orthorhombic structural unit cell with lattice parameters a ≈ b ≈ 5.56 Å and c ≈ 7.05 Å for NaFeAs 30,31 . The momentum transfer Q = Ha * + Kb * + Lc * is denoted as Q = (H, K, L) in reciprocal lattice units (r.l.u.), with a * = b a2π/a, b * = b b2π/ b, and c * = b c2π/c. In this notation, magnetic Bragg peaks are at Q = (1, 0, L), with L = 0.5, 1.5, 2.5, …. Samples were mostly aligned in the [1, 0, 0] × [0, 0, 1] scattering plane, which allows scans of magnetic peaks along H and L; the x = 0.012 sample was also studied in the [1, 0, 1.5] × [0, 1, 0] plane. We have carried out neutron diffraction, neutron Larmor diffraction, and inelastic neutron scattering experiments on NaFe 1−x Ni x As (see Methods section for experimental details). Figure 1d shows the overall phase diagram determined from our experiments, with T s , T N , and T c marked. Although for optimal-doped and over-doped regimes, the samples on average exhibit a tetragonal structure at all temperatures, there are local orthorhombic distortions resulting in broadening of d-spacing distribution that can be characterized by δ or Δd/d. The orthorhombic distortion δ is plotted in a pseudo-color scheme as a function of temperature and doping near optimal-doping in Fig. 1d. Figure 1e shows the Ni-doping dependence of the ordered magnetic moment and δ at T = 5 K, and T = T c for superconducting samples. With increasing Ni-doping x, the AF ordered moment and T N decrease monotonically, and no magnetic order is detected in the x = 0.015 sample [ Fig. 1e]. While the magnetic order parameter for the x = 0.004 sample resembles that of NaFeAs [ Fig. 2e, f], the magnetic order becomes strongly suppressed upon entering the superconducting state for x = 0.010 [ Fig. 2g], similar to other iron pnictides 32,33 .
Reentry into the paramagnetic state in NaFe 1−x Ni x As with x = 0.012. For the x = 0.012 sample, magnetic order begins at T N ≈ 19 K and becomes strongly suppressed upon entering the superconducting state below T c and reenters into the paramagnetic state without any long-range order below T r ≈ 10 K [ Fig. 2h]. Given the sharp superconducting transition at T c (Supplementary Fig. 1 and Methods section), T r is well inside the superconducting state. This is similar to the behavior of nearly optimal-doped Ba(Fe 0.941 Co 0.059 ) 2 As 2 34 , although AF order in Ba (Fe 0.941 Co 0.059 ) 2 As 2 is short-range and incommensurate 21 . To confirm that the magnetic order in our x = 0.012 sample is longrange and commensurate, we carried out scans along [H, 0, 1  Fig. 1 The phase diagram of NaFe 1−x Ni x As determined from neutron scattering measurements. a Crystal structure and magnetic order of NaFeAs. b Schematic evolution of NaFe 1−x Ni x As in tetragonal, paramagnetic orthorhombic, and AF orthorhombic states. c Schematic of how d-spacing distribution changes from a tetragonal state at high temperatures to an orthorhombic state through phase transition [characterized by δ = (a − b)/(a + b)] or a state with local orthorhombic distortions (characterized by broadening of the d-spacing distribution Δd/d) that on an average remains tetragonal. For the orthorhombic state, when the splitting δ is too small to be resolved experimentally, a broadening is also observed (red dashed line). In such cases, the two situations can nonetheless be differentiated by examining the temperature dependence of δ or Δd/d. d The phase diagram of NaFe 1−x Ni x As. T s , T N , and T c are the transition temperatures for the tetragonal-to-orthorhombic structural phase transition, the AF phase transition, and the superconducting transition. The point for x = 0 is obtained from ref. 25 . e The Ni-doping dependence of the ordered magnetic moment and orthorhombic distortion δ at T = 5 K and T = T c for superconducting samples. The error bars in (d) are estimated errors from fits-to-order parameters and transition temperatures correlation lengths exceeding 100 Å) for the x = 0.012 sample near optimal superconductivity before disappearing near x = 0.015. These wave-vector scans also confirm the complete disappearance of long-range magnetic order below T r . For comparison, we note that magnetism in electron-doped Ba(Fe 1 −x Co x ) 2 As 2 (~6%) 21 , BaFe 2−x Ni x As 2 (~5%) 22 , and NaFe 1−x Co x As (~2.3%) 25 exhibits cluster spin glass and incommensurate magnetic order near optimal superconductivity likely related to impurity effects 23,35 . The absence of such behavior in NaFe 1 −x Ni x As is likely a result of significantly lower dopant concentration in NaFe 1−x Ni x As (~1.3%) near optimal doping. Our inelastic neutron scattering measurements on the x = 0.012 sample confirm that the presence of a neutron spin resonance, which can act as a proxy for the superconducting order parameter, is unaffected when cooled below T r (Supplementary Fig. 2 and Methods section).
Nematic order and local orthorhombic distortions in NaFe 1 −x Ni x As. Having established the evolution of AF order and its interplay with superconductivity in NaFe 1−x Ni x As, we examined the Ni-doping evolution of the nematic order in NaFe 1−x Ni x As. To precisely determine the evolution of orthorhombic distortion, we used high-resolution neutron diffraction and neutron Larmor diffraction to investigate the temperature evolution of the orthorhombic lattice distortion (Supplementary Figs. 3 and 4 and Methods section). For NaFe 1−x Ni x As with x ≤ 0.013, we can see clear orthorhombic lattice distortion below T s , also confirmed by the anomalies in temperature dependence of electrical resistivity measurements (Supplementary Fig. 5 and Methods section). Figure 3a-c shows temperature and Ni-doping dependence of the orthorhombic distortion δ. For NaFe 1−x Ni x As with x ≤ 0.013 at temperatures above T s , and for x ≥ 0.015 at all temperatures, the system is on an average tetragonal and should in principle have δ = 0. Surprisingly, we see clear temperature-dependent δ. Moreover, while δ below T s behaves as expected for an order parameter associated with phase transition, δ in temperature regimes with an average tetragonal structure exhibits a Curie-Weiss temperature dependence, suggesting that it arises from local orthorhombic distortions. In all cases, we find that δ decreases dramatically below T c , indicating that orthorhombic distortion, whether long-range or local, competes with superconductivity. The competition between superconductivity and long-range nematic order is similar to Ba(Fe 1−x Co x ) 2 As 2 36 and can be captured by a phenomenological Landau theory, based on an effective action in terms of the corresponding order parameters (see Methods section): where, the last term describes the competition between nematicity and superconductivity. As a result, the nematic order parameter is noticeably suppressed inside the superconducting phase, compared with its value (δ 0 ) in the normal phase, so that (see Methods section for the derivation) whereas the superconducting order parameter itself remains essentially unchanged due to tiny values of δ 0 (see Eq. (8) in Methods section). In the tetragonal phase (δ = 0), the competition between local orthorhombic distortions and superconductivity is reflective of the suppression of nematic susceptibility below T c 37 . We emphasize that the local orthorhombic distortions we uncovered in the tetragonal phase of NaFe 1−x Ni x As are distinct from the phase separation into superconducting tetragonal and AF orthorhombic regions found in Ca(Fe 1Àx Co x ) 2 As 2 under biaxial strain 38,39 . In the latter compound, the quantum phase transition between the superconducting tetragonal and AF orthorhombic phases is of first order, and the resulting phase separation into these two phases with different in-plane lattice parameters allows the material to respond to biaxial strain in a continuous fashion; this would occur even if there were no quenched disorder. In NaFe 1−x Ni x As, the quantum phase transition is of second order and, therefore an analogous phase separation does not occur. Instead, the local orthorhombic distortions we observe in NaFe 1−x Ni x As likely result from large nematic susceptibility near optimal superconductivity pinned by quenched disorder.  Given that the orthorhombic distortion with Curie-Weiss temperature dependence arises from local orthorhombic distortions, an alternative way to characterize such distortion is broadening of d-spacing distribution width, Δd/d (see Methods section). In Fig. 4a-d, we show Δd/d in NaFe 1-x Ni x As, obtained from our neutron Larmor diffraction measurements. Given that the local orthorhombic distortions arise from quenched disorder coupled with large nematic susceptibility near a nematic QCP, it should track the temperature dependence of nematic susceptibility, since the quenched disorder should depend weakly on temperature. Therefore, we have fitted Δd/d in Fig. 4a-d with the Curie-Weiss form Δd/d ∝ 1/(T − T * ) and extracted the Weiss temperature T * as a function of doping, as shown in Fig. 4e. Our Δd/d results are well described by the Curie-Weiss form, with T * changing from positive in underdoped to negative in overdoped regime [ Fig. 4e], suggesting a nematic QCP near optimal superconductivity. These results are reminiscent of temperature and doping dependence of nematic susceptibility from elastoresistance 6 and shear modulus measurements 8 , suggesting that temperature dependence of Δd/d is a direct measure of the nematic susceptibility without the need to apply external stress.

Discussion
In NaFe 1−x Ni x As, the orthorhombic distortion and the structural phase transition temperature are δ ≈ 1.7 × 10 −3 and T s ≈ 58 K for x = 0 25,31 ; for x = = 0.012, they become δ ≈ 7 × 10 −4 and T s ≈ 33 K. We find no evidence of structural phase transitions for samples with x ≥ 0.015, suggesting the presence of a putative nematic QCP at x = x c , where x c ≳ 0.015. These results are consistent with recent Muon spin rotation and relaxation study of the magnetic phase diagram of NaFe 1−x Ni x As 40 . The doping-dependence of T s and δ are also consistent with the Ni-doping dependence of T * determined from Curie-Weiss fits to temperature dependence of Δd/d, which changes from positive to negative near x ≈ 0.015 [ Fig. 4e]. Since our neutron Larmor diffraction measurements were carried out using polarized neutron beam produced by an Heusler monochromator, which has an energy resolution of about ΔE ≈ 1.0 meV 28,29 , the local orthorhombic distortions captured in our measurements are either static or fluctuating slower than a time scale of τ $ h=2ΔE $ 0:3 ps, where ħ is the reduced Planck's constant 41,42 . One possible origin of such slow fluctuations may be the in-plane transverse acoustic phonons that exhibit significant softening in the paramagnetic tetragonal phase when approaching a nematic instability 43 . Future neutron scattering experiments with energy resolutions much better than ΔE ≈ 1 meV are desirable to separate the static and slowly fluctuating contributions. Our results also indicate that the nematic QCP would occur at a x value that is distinctively larger than that of the magnetic QCP in the absence of superconductivity. In the phase diagram of iron pnictides with decoupled T s and T N , due to the competition between superconductivity with both nematic and magnetic orders, magnetic order forms a hump peaked at T c near optimal doping [ Fig. 1d], and the structural phase transition disappears in a similar fashion at a larger x.
Theoretically, a determinantal quantum Monte Carlo study of a two-dimensional sign-problem-free lattice model reveals an Ising nematic QCP in a metal at finite fermion density 44 . In the nematic phase, the discrete lattice rotational symmetry is spontaneously broken from fourfold to twofold, and there are also nematic quantum critical fluctuations above the nematic ordering temperature. Within the numerical accuracy of the determinantal quantum Monte Carlo study, the uniform nematic susceptibility above the nematic ordering temperature has Curie-Weiss temperature dependence, signaling an asymptotic quantum critical scaling regime consistent with our observation 44 . Alternatively, the observed Curie-Weiss temperature-dependent behavior of nematic susceptibility can be understood from spin-driven nematic order theory, where magnetic fluctuations associated with static AF order induce formation of the nematic state 45 . In this picture, the effect of lattice strain coupled to the nematic order parameter produces a mean-field Curie-Weiss-like behavior, arising from the nemato-elastic coupling which has direction-dependent terms in the propagator for nematic fluctuations. The Curie-Weiss temperature-dependent nematic susceptibility should occur in the entire phase diagram, where there is a significant softening of the elastic modulus 45 . This means that Curie-Weiss temperature dependence of local orthorhombic distortions that we observe is a signature of nemato-elastic coupling, which does not suppress the magnetic fluctuations that cause the nematic order, but transforms the Ising nematic transition into a mean-field transition 45 . T c ≈ 17K δ is obtained by assuming that it is 0 at T = 50 K, and broadening at lower temperatures are fit with two split peaks, with widths of the single peak at T = 50 K. Open symbols correspond to measurements where a splitting is definitively observed, and solid symbols represent measurements that only resolve a broadening due to experimental limitations (Methods section and Supplementary Fig. 4). All vertical error bars in the figure represent least-square fits to the raw data with errors of 1 s.d Our discovery of local orthorhombic distortions exhibiting Curie-Weiss temperature dependence across the phase diagram of NaFe 1−x Ni x As results from the proliferation of nematic fluctuations and large nematic susceptibility near the nematic QCP. Quenched disorder that are always present in such doped materials act to pin the otherwise fluctuating local nematic domains, resulting in static (or quasi-static) local orthorhombic distortions that can lead to observations of rotational symmetry breaking seen with multiple probes [13][14][15][16][17][18][19] . We have definitively observed local nematic distortions in NaFe 1−x Ni x As that are static or quasistatic, in contrast to local distortions seen in Sr 1−x Na x Fe 2 As 2 , using pair distribution function analysis that contain significantly more dynamic contributions 46 and which would not cause rotational symmetry breaking seen by static probes. Our observation of local nematic distortions highlights the presence of nematic fluctuations near the nematic QCP, which can play an important role in enhancing superconductivity of iron pnictides 9-12 , while the intense Ising nematic spin correlations near the nematic QCP may be the dominant pairing interaction [47][48][49] .

Methods
Elastic neutron scattering experimental details. Elastic neutron experiments were carried out on the Spin Polarized Inelastic Neutron Spectrometer (SPINS) at the NIST Center for Neutron Research, United States and the HB-1A triple-axis spectrometer at the High-Flux-Isotope Reactor (HFIR), Oak Ridge National Laboratory (ORNL), United States. We used pyrolytic graphite [PG(002)] monochromators and analyzers in these measurements. At HB-1A, the monochromator is vertically focused with fixed-incident neutron energy E i = 14.6 meV and the analyzer is flat. At SPINS, the monochromator is vertically focused and the analyzer is flat with fixed-scattered neutron energy E f = 5 meV. A PG filter was used at HB-1A and a Be filter was used at SPINS to avoid contamination from higherorder neutrons. Collimations of 40′-40′-sample-40′-80′ and guide-40′-sample-40′open were used on HB-1A and SPINS, respectively.
To measure the structural distortion in NaFe 1−x Ni x As (x = 0.012) at SPINS, we changed the collimation to guide-20′-sample-20′-open to improve the resolution and removed the Be filter. Our measurement was carried out nominally around a weak nuclear Bragg peak Q = (2, 0, 0), but the measured intensity at this position mostly come from higher-order neutrons [Q = (4, 0, 0) for λ/2 neutrons and Q = (6, 0, 0) for λ/3 neutrons]. While we do not resolve two split peaks in the orthorhombic state, clear broadening can be observed. Typical scans along the [H, 0, 0] direction centered at Q = (2, 0, 0) are shown in Supplementary Fig. 3. δ in Fig. 3b is obtained by assuming δ = 0 at T = 50 K and fitting broadening at lower temperatures as two split peaks with fixed widths of the peak at T = 50 K.
Inelastic neutron scattering experimental details. Our inelastic neutron scattering experiment was carried out on the HB-3 triple-axis spectrometer at HFIR, ORNL, United States. Vertically focused pyrolytic graphite [PG(002)] monochromator and analyzer with fixed-scattered neutron energy E f = 14.7 meV were used. A PG filter was used to avoid higher-order neutron contaminations. The collimation used was 48′-40′-sample-40′-120′.
Using inelastic neutron scattering, we studied the neutron spin resonance mode 2,50 in NaFe 1−x Ni x As, with x = 0.012. Energy scans at Q = (1, 0, 0.5) above (T = 35 K) and below T c (T = 1.5 and 9 K) are compared in Supplementary Fig. 2a. The scans below T c after subtracting the T = 35 K scan are compared in Supplementary Fig. 2b. A clear resonance mode at E r = 7 meV similar to optimal- doped NaFe 1−x Co x As 51 is observed, with almost identical intensities at T = 1.5 and 9 K. Constant energy scans along [H,0,0.5] at different temperatures are compared in Supplementary Fig. 2c, confirming the results in Supplementary Figs. 2a, b. Temperature dependence of the resonance mode is shown in Supplementary  Fig. 2d, over-plotted with temperature dependence of orthorhombicity and AF order parameter. Intensity of the resonance mode increases smoothly below T c and T r , displaying no response when AF order is completely suppressed below T r . These results demonstrate the coexistence of robust superconductivity and nematic order without AF order in NaFe 1−x Ni x As (x = 0.012) below T r .
Larmor diffraction experimental details. Our neutron Larmor diffraction measurements were carried out at the three axes spin-echo spectrometer at Forschungs-Neutronenquelle Heinz Maier-Leibnitz (MLZ), Garching, Germany. Neutrons are polarized by a super-mirror bender, and higher-order neutrons are eliminated using a velocity selector. We used double-focused PG(002) monochromator and horizontal-focused Heusler (Cu 2 MnAl) analyzer in these measurements. Incident and scattered neutron energies are fixed at The detailed principles of neutron Larmor diffraction has been described in detail elsewhere 29,52 . In such experiments, polarization of the scattered neutrons P is measured as a function of the total Larmor precession phase ϕ tot . By analyzing the measured P(ϕ tot ), information about the sample's d-spacing distribution can be obtained.
For an ideal crystal with d-spacing described by a δ function, P is independent of ϕ tot , with P(ϕ tot ) = P 0 . P 0 accounts for the non-ideal polarization of neutrons and can be corrected for by Ge crystal calibration measurements. In real materials, due to internal strain and sample inhomogeneity, or in the case of iron pnictides, a twinned orthorhombic phase, the d-spacing should instead be described by a distribution f(ϵ), with ϵ = δd/d. δd is the deviation from the average d-spacing d. P (ϕ tot ) is then described by Thus, P(ϕ tot ) can be regarded as the Fourier transform of the lattice d-spacing distribution f(ϵ). By measuring P(ϕ tot ), it is possible to resolve features with a resolution better than 10 −5 in terms of ϵ, limited by the range of accessible ϕ tot .
The distribution of d-spacing f(ϵ) is commonly described as a Gaussian function with full-width-at-half-maximum (FWHM) ϵ FWHM , also denoted as Δd/d in the rest of the paper. Eq. (3) then becomes In iron pnictides with a nonzero nematic order parameter, due to twinning, f(ϵ) becomes the sum of two Gaussian functions. Assuming that the two Gaussian peaks have identical FWHM ϵ FWHM , Eq. (4) becomes where, r and (1 − r) denotes the relative populations of the two lattice d-spacings a and b, and Δϵ = 2(a − b)/(a + b) = 2δ 53 . Therefore, the nematic order parameter can be extracted by fitting P(ϕ tot ) using Eq. (5). When δ is too small to be directly resolved by Larmor diffraction, P(ϕ tot ) can be well described by either Eq. (4) or (5). In such cases, we either extract Δd/d from Eq. (4) (Fig. 4) or extract δ by assuming at T = 50 K, δ = 0 and extract ϵ FWHM , then fit to Eq. (5) by fixing ϵ FWHM to this value (Figs. 1e, 3). Measurements of P(ϕ tot ) at several different temperatures for NaFe 1−x Ni x As (x = 0.013) are shown in Supplementary Fig. 4, and fit to Eq. (5) as described.
A key feature of Eq. (5) is an oscillation in P(ϕ tot ), which can be seen in raw data in Supplementary Fig. 4d-i (open symbols in Fig. 3c); in these cases, the measurement provides definitive evidence of an orthorhombic state. For other panels in Supplementary Fig. 4, due to limited range of ϕ tot , P(ϕ tot ) can be equally well described by Eq. (4) (solid symbols in Fig. 3c); for such data, we cannot differentiate between a true splitting and a broadening from measurements done at a single temperature.
Magnetic susceptibility and electrical resistivity measurements. To ensure that T r for NaFe 1−x Ni x As (x = 0.012) is well inside the superconducting state, we show in Supplementary Fig. 1 its magnetic susceptibility as a function of temperature. As can be seen, the sample displays a sharp superconducting transition at T c ≈ 17 K, with a width ΔT c ≈ 2 K. T r is well inside the superconducting state, unaffected by the width of the superconducting transition.
The temperature and doping dependence of the in-plane electrical resistivity ρ (T) were measured using the standard four-probe method, the results are normalized to ρ(200 K) and summarized in Supplementary Fig. 5. The superconducting transitions for all measured samples are sharp. The kinks associated with the structural transition at T s can be clearly identified in underdoped samples ( Supplementary Fig. 5a-d), similar to NaFe 1−x Cu x As 54 . These kinks are progressively suppressed with increasing Ni concentration and disappear in overdoped samples. T s determined from electrical resistivity measurements are in good agreement with those obtained from Larmor diffraction.
Coexistence of superconductivity with lattice nematicity. Here, we first consider the case without any long-range magnetic order, as is realized in NaFe 1 −x Ni x As for x > 0.012. In that case, the effective Landau free energy can be written in terms of only the superconducting order parameter Δ and the orthorhombicity δ≡(a − b)/(a + b): Here, we assume that the superconducting order parameter transforms under the tetragonal point symmetry, i.e., it does not break the C 4 rotational symmetry of the lattice. Since the lattice-nematic order parameter breaks this symmetry, the coupling to superconductivity is quadratic in δ. Above, the coefficient C is in fact the elastic shear modulus C 66 , which is the inverse of the nematic susceptibility. The latter has a Curie-Weiss behavior (see Fig. 4 in the main text): Here, C ð0Þ 66 is the "bare" value of shear modulus in the absence of nematic transition. Note that, the above formula can been derived rigorously from an effective model of lattice orthorombicity δ coupled with an electronic nematic order parameter 29 .
Here, we simply take T * to be the phenomenological Curie-Weiss temperature extracted from fitting the d-spacing in Fig. 4e. Note that, if T * is positive (for x < 0.016), we identify it with the nematic transition temperature T s such that 0 > C = − C j j is below T s . Minimizing this effective action with respect to the two-order parameters ∂F/ ∂Δ = 0 = ∂F/∂δ we obtain in the mixed state with T < {T s , T c } nonzero values of both parameters: where, Δ 0 = ffiffiffiffiffiffiffi ffi α=β p and δ 0 = ffiffiffiffiffiffiffiffiffiffiffiffi C j j=D p are the values of the order parameters in the absence of coupling between them. In the coexistence phase, the free energy becomes: where, F ð0Þ SC ¼ Àα Δ 0 j j 2 =4. Note that, for the coexistence phase to be stable, the last term in the above expression must be positive, which is only possible if 4γ 2 βD <1, or equivalently, βD > 4γ 2 .
There is no perceptible change in the superconducting transition temperature below T s , implying Δ j j ' Δ 0 j j. Substituting this into Eq. (8), we obtain: By contrast, the suppression of the orthorhombicity below T c is substantial, δ ≈ 0.5δ 0 (see Fig. 3b, c), meaning that 2γ D À Á Δ 0 j j 2 % δ 2 0 from Eq. (9). Substituting this into Eq. (11), we obtain: in other words, we can approximate the denominator in Eqs. (8) and (9) to be 1. This is also consistent with the requirement from Eq. (10) for the coexistence phase to be stable. In summary, the phenomenological Landau free energy explains qualitatively the experimental data in the coexistence phase of superconductivity and nematicity. Furthermore, comparison with the experiment allows us to impose strong condition on the smallness of the coupling constant γ in terms of inequality (12).
Coexistence of three phases. Below x > 0.012, NaFe 1−x Ni x As has a long-range AF order, and the free energy in Eq. (10) has to be modified to include the magnetic order parameter M: where, we have included phenomenological coupling constants μ and w. The sign of w is positive, in accord with our experimental observation that AF order and superconductivity compete with each other (see Fig. 2g, h in the main text). The sign in front of μ on the other hand is negative, indicating magnetoelastic coupling that favors the coexistence of magnetism and orthorhombic distortion. Because of this coupling, it is clear that δ will acquire an additional component proportional to M 2 inside the AF phase: because M 2 and Δ j j 2 repel each other via the last term in Eq. (13), this implies, in view of Eq. (14), that a new term proportional to ΔF / δ j j Δ j j 2 will be generated in the action, coupling the square of the superconducting order parameter linearly to the lattice orthorhombicity.
Working with full free energy in Eq. (13) is impractical because of the large number of phenomenological parameters that are difficult to determine experimentally. Nevertheless, it offers a qualitative insight into the coexistence between AF, lattice nematicity, and superconductivity, as the above discussion shows.
As a parenthetical remark, we note that the term Àμ δ j j Á M 2 in free energy may appear surprising at first sight, as one might have expected that lattice distortion and magnetization should couple biquadratically. The reason for linear coupling is because the stripe AF order in iron pnictides breaks the lattice C 4 symmetry, as does the shear strain δ 29,55-57 . Note that this conclusion holds independently of whether the microscopic origin of nematicity is purely magnetic 55,56 or due to orbital ordering of Fe d xz /d yz orbitals 57-60 .
Data availability. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.