The hyperfine coupling is given by: H h y p = ∑ i ∈ n . n . I ̂ ⋅ A i ⋅ S ̂ i , (1) where I ̂ α are the nuclear spin operators, S ̂ i are the Fe spins on the four nearest neighbors to the As nucleus, and the hyperfine coupling tensors, A i , have dipolar symmetry [18, 19]. In BaFe₂As₂, the eigenvalues of the tensor are A _aa = A _bb = 6.6 kOe/μ _B , and A _cc = 4.7 kOe/μ _B , but these values depend on details of the electronic structure and can vary between compounds. They can change as a function of pressure or doping, but the general symmetry of the A i does not. The magnitude of these values are greater than those for the direct dipole interaction, which reflects the electronic hybridization between the As 4p and Fe 3d orbitals. The tensor nature of this coupling gives rise to a hyperfine field H _hf ∥c when the Fe moments lie in the ab plane, as illustrated in top row of Figure 2. Alternatively, if the Fe moments ordered parallel to c, then H _hf lies in the plane, as shown in the middle row of Figure 2.

Other magnetic orderings can be present in the pnictides, beyond the stripe-type spin density wave illustrated in Figure 2. In particular, the spin-vortex crystal configuration, in which the Fe spins form either loop or hedgehog structures, can be stabilized when a glide symmetry across the Fe planes is broken [20]. This type of order was recently identified in Ni-doped and Co-doped CaKFe₄As₄ [21]. There are two crystallographically distinct As sites in this structure, and the hyperfine field either vanishes at both sites for loop spin-vortex crystal order, or is aligned along the c axis for one As site and vanishes for the other As site for hedgehog spin-vortex crystal order. In this system, doping induces the latter, which was identified by NMR.

Note that when the Fe is substituted by another transition metal, the symmetry among the hyperfine couplings of the four nearest neighbors is broken, and the hyperfine field vector can acquire a component perpendicular to the c axis [22]. This effect contributes to inhomogeneous broadening in doped samples. In the presence of nematic order, the C ₄ symmetry of the lattice is broken, and in principle A _aa ≠ A _bb . This asymmetry has not been observed directly. In FeSe, the Knight shift tensor becomes anisotropic, with K _aa ≠ K _bb , where K _αα = A _αα χ _αα , and χ _αα is the static magnetic susceptibility [5]. However, in this case the anisotropy is believed to reflect that of the susceptibility, i.e., χ _aa ≠ χ _bb in the nematic phase, rather than an asymmetry in the hyperfine coupling [23]. Curiously, χ _aa > χ _bb in FeSe, but χ _bb > χ _aa in BaFe₂As₂, LaFeAsO, and NaFeAs [24–26], which may reflect different natures of the nematic order in these materials [27]. In particular, the nematic order parameter in the latter can be understood in terms of differential occupations between the d _yz and d _xz orbitals, whereas in the former the nematic order may involve more complex superpositions between 3d orbitals, giving rise to bond-centered nematic order [23, 28, 29]. The sign of the resistivity anisotropies in BaFe₂As₂ and FeSe match those of the magnetic anisotropies [24, 30], suggesting a common origin for the two.

For nuclei with spin I > 1/2, the quadrupolar coupling is given by: H Q = h ν z z 6 3 I ̂ z 2 − I ̂ 2 + η I ̂ x 2 − I ̂ y 2 , (2) where the EFG tensor is given by: ν α β = e Q 12 h ∂ 2 V ∂ x α ∂ x β , (3) and V is the electrostatic potential at the As site. This quantity is dominated by the occupation of the As 4p orbitals, which in turn are hybridized with the d _xz,yz -orbitals of the neighboring Fe atoms [31]. The EFG asymmetry parameter is given by: η = ν y y − ν x x ν x x + ν y y . (4)

This quantity vanishes in the tetragonal phase because the As 4p _x and 4p _y orbitals are degenerate, hence ν _xx = ν _yy [32]. The EFG tensor is always traceless, thus in this case it can be characterized by a single quantity, ν _zz . The magnitude of ν _zz depends on the degree of hybridization of the As 4p _z orbital, and is sensitive to the c-axis length of the unit cell [33–35] and also varies with temperature [18, 36–38].

The EFG is strongly affected by doping. Replacing the Fe by another transition metal, or replacing the As by P, will give rise to spatially-varying strain fields that will distort the electrostatic potential and alter the EFG tensor. This effect can lower the symmetry such that ν _xx ≠ ν _yy , and create non-zero off-diagonal terms of the tensor. In some cases the effect is sufficiently large that separate resonances can be detected for sites that are adjacent to the substitutional site [31, 39]. This sensitivity to disorder tends to significantly broaden the NMR resonance frequencies in doped samples relative to the parent compound. In doped Ba(Fe,Co)₂As₂ the width of the quadrupolar satellites is dominated by a variation of η that exhibits a Curie-Weiss temperature dependence, which has been associated with the growth of nematic correlations surrounding the dopant atoms [40]. A similar effect has been observed in FeSe where the magnetic linewidth is broadened by the presence of crystal defects (there is no quadrupolar interaction for Se) [41]. In this case, the Knight shift anisotropy reflects the growth of an Edwards-Anderson order parameter of the nematicity [42].

In the presence of long range nematic order, the C ₄ symmetry of the EFG tensor is broken and ν _xx ≠ ν _yy [43]. In CaFe₂As₂ η reaches 0.4 at low temperature, but in BaFe₂As₂ and SrFe₂As₂, η reaches 1.2 and 1.4, respectively, reflecting the fact that ν _xx and ν _yy have opposite signs [33]. Note that it is sometimes customary to define the principal axes of the EFG such that the eigenvalues |V _zz |≥|V _yy |≥|V _xx |, so that 0 ≤ η ≤ 1. Here we do not follow this convention, maintaining the V _zz associated with the crystal c axis as in the tetragonal phase. This choice is the reason that η > 1, which would correspond to a confusing rotation of orientation of the principal axes under the customary definition. These changes reflect a dramatic rearrangement of the charge distribution around the As nucleus below the structural transition, despite an orthorhombicity of only approximately 1%. The EFG tensor is dominated by the occupations of the As 4p orbitals, which in turn are hybridized with the Fe 3d _yz and 3d _xz orbitals [31, 32]. These orbitals are degenerate in the tetragonal phase, but develop a splitting on the order of 40 meV in the nematic phase [44]. As a result, the relative occupations of the As p _x and p _y orbitals change and alter the EFG.

In order to resolve the in-plane asymmetry that emerges in the Knight shift or EFG tensors, it is important to measure the NMR spectra of carefully oriented single crystals. In an external magnetic field, the As nuclear spins experience the sum of three interactions: H Z + H h y p + H Q , where H Z = − γ ℏ I ̂ ⋅ H 0 is the Zeeman interaction. The spectrum will consist of three resonances. For sufficiently large Zeeman interaction (γH ₀ ≫ ν _αα ), the frequencies are given to first order by: f s a t 1,2 = γ H 0 1 + K θ , ϕ ± ν z z 3 ⁡ cos 2 ⁡ θ − 1 + η ⁡ sin 2 ⁡ θ ⁡ cos 2 ϕ f c e n = γ H 0 1 + K θ , ϕ + Δ f 2 θ , ϕ (5) where f _sat1,2 and f _cen are the frequencies of the two satellites and the central transition, θ is the angle between H ₀ and c, and ϕ is the azimuthal angle of the projection of H ₀ in the ab plane. The Knight shift is given by: K θ , ϕ = K z z ⁡ cos 2 ⁡ θ + K x x ⁡ sin 2 ⁡ θ ⁡ cos 2 ⁡ ϕ + K y y ⁡ sin 2 ⁡ θ ⁡ sin 2 ⁡ ϕ . (6)

The second order correction to the central transition is given by: Δ f 2 θ , ϕ = ν z z 2 96 γ H 0 7 η 2 − 18 − 6 η 2 − 30 cos 2 ⁡ θ − 9 18 + η 2 cos 4 ⁡ θ + 12 η ⁡ cos 2 ϕ 9 ⁡ cos 4 ⁡ θ − 8 ⁡ cos 2 ⁡ θ − 1 − 9 η 2 ⁡ cos 4 ϕ sin 4 ⁡ θ . (7)

Measurements of the two satellite transitions as a function of field orientation can be used to extract the Knight shift and EFG tensors.

In the tetragonal phase, there are three distinct resonances for each As in the unit cell (e.g., one As site in AFe₂As₂). However, in the nematic phase there are multiple domains in an untwinned crystal, giving rise to extra resonances. If the magnetic field is oriented in the plane, then there may be two sets of resonances associated with these two domains because the azimuthal angle, ϕ, differs for each domain, as observed in Figure 3. This behavior is illustrated in Figure 4. In general, there are five distinct transitions (six if the Knight shift is also anisotropic), but in certain directions the spectra from the two domains overlap and three resonances are recovered. Large values of |η| cause the satellites to switch positions from below to above the central transition as the in-plane angle ϕ is rotated.

When η is non-zero the nuclear spin Hamiltonian does not commute with I ̂ z , and the eigenstates, |i⟩, are superpositions of the |I _z ⟩ states. In such cases it may be possible to observe hidden transitions when the matrix element 〈 ψ i | I ̂ ± | ψ j 〉 does not vanish. This effect can arise when ν _zz ≳ γH, but in almost all cases there is an internal field, H _int , that develops in the antiferromagnetic state that accompanies the orthorhombic distortion. Typically γH _int ≫ ν _zz is sufficiently large that such hidden transitions should not be evident. For CaFe₂As₂, however, ν _zz /γH ∼ 0.7 thus it may be possible to observe these transitions for small external fields H ₀ applied in the plane, as illustrated in Figure 5. To date no such transitions have been reported.

The hyperfine field, H _hyp , and the EFG not only give rise to static spectral signatures, but their fluctuations can also drive relaxation of the nuclear spins. The nuclear spin energies are generally several order of magnitude smaller than the relevant electron energies, and thus the nuclear spin ensemble can easily be manipulated to have a non-equilibrium distribution among the energy levels. The nuclear spins will return to thermal equilibrium with the “lattice,” or the electronic spin system, over a time scale, T ₁, known as the spin-lattice relaxation time. This rate, T 1 − 1 , can be measured with high precision and is determined by the couplings described above in Section 2.1, Section 2.2. In essence, the nuclear spins act as quantum sensors that probe the noise spectra of the electronic spin and charge degrees of freedom. There are three important spectral densities that may play a role: J ω = γ 2 2 ∫ ∞ ∞ h + , h − e i ω τ d τ (8) J 1,2 ω = e Q ℏ 2 ∫ ∞ ∞ V + 1,2 , V − 1,2 e i ω τ d τ (9) where A , B = 〈 A ( τ ) B ( 0 ) + B ( τ ) A ( 0 ) 〉 / 2 is an ensemble average. Here h _± = H _hyp,x ± iH _hyp,y , V _{± 1} = V _xz ± iV _yz , and V _{± 2} = (V _xx − V _yy )/2 ± iV _xy [45]. J(ω) and J ^(1,2)(ω) are the spectral densities of the hyperfine and EFG fluctuations, respectively. The T 1 − 1 measurements probe these spectral functions at the nuclear resonance frequency, ω _N , typically on the order of 100 MHz ( ∼ 0.4 μ eV). This frequency is nearly vanishing in comparison to the relevant electronic energies, typically on the order of meV. If there is only one type of fluctuation, then T 1 − 1 is directly related to the appropriate spectral density. For example, for pure magnetic fluctuations, (e.g., V _{± 1,2} = 0), then T 1 − 1 = J ( ω N ) , or if h _± = 0 and V _{± 1} = 0, then T 1 − 1 = J ( 2 ) ( ω N ) , or if h _± = 0 and V _{± 2} = 0 then T 1 − 1 = J ( 1 ) ( ω N ) . On the other hand, if multiple types of fluctuations are present, then the nuclear spin relaxation is a complex function of multiple spectral functions and it is not straightforward to disentangle the different contributions to T 1 − 1 [46–48]. It is common for the magnetic fluctuations to dominate, in which case the contribution from quadrupolar contributions can be ignored. However, as we discuss below in Section 2.4.1, there is a significant quadrupolar contribution to the T 1 − 1 of the As in the pnictides.

When the electronic system is close to an thermodynamic instability, critical fluctuations will slow down and may have a dramatic effect on T 1 − 1 . As an example, consider the case where the magnetic fluctuations exhibit Lorentzian autocorrelation function: h + , h − = h 0 2 e − | τ | / τ c , where τ _c is the correlation time. In this case: 1 T 1 = γ 2 h 0 2 τ c 1 + ω N 2 τ c 2 . (10)

Usually ω _N τ _c ≪ 1, but close to a phase transition where there is critical slowing down, τ _c can approach infinity as a power law. As a result, T 1 − 1 exhibits a peak at the phase transition temperature, as illustrated in Figure 6. This property has been exploited to investigate the critical dynamics of many correlated electron systems, especially the pnictides [39, 49].

As discussed above in Section 2.3, the nuclear spin eigenstates become superpositions of the I _z states for non-zero values of η. In this case, it is important to consider spectral densities other than those described in Eqs 6, 7 to account for fluctuations in other directions. For example, hyperfine field fluctuations, H _hyp,z , may play a role in relaxing the nuclear spins. In such cases the form of the relaxation function is more complicated, and Eq. 8 will contain more terms.

In BaFe₂As₂, uniaxial strain dramatically alters the EFG tensor and η ∝ɛ _xx , as illustrated in Figure 7 [32]. In the absence of strain, the crystal remains tetragonal down to the nematic transition temperature T _s = 135 K, in which case the EFG tensor spontaneously develops an asymmetry such that ν _xx ≠ ν _yy . However, in the presence of a finite strain field, η becomes finite above T _s . In fact, η ∝ɛ _xx and the constant of proportionality is the nematic susceptibility. As seen in Figure 7, the response is largest just above T _s , where the nematic susceptibility diverges. This behavior is illustrated in Figure 8, which reveals that dη/dɛ _xx ≈ 300 at T _s .

Uniaxial stress has also been utilized to study the NMR spectra in LaFeAsO and FeSe [23]. In this case, compressive stress on the order of 10–20 MPa was applied by tightening a screw. A similar approach in which a crystal is suspended between the prongs of a horseshoe-shaped device was also utilized for NMR studies under tensile strain [52]. Neither approach directly measures the strain, precluding the possibility for feedback control. Moreover, differential thermal contraction may also give rise to finite strain fields even without applied stresses via the tightening screw. Nevertheless, such approaches can be used to mechanically detwin the crystal in the nematic phase [53]. If the crystal is cooled in a finite strain field, then a single nematic domain should nucleate below the transition. For a twinned crystal, there should be two sets of resonances as discussed above in Section 2.3. In detwinned or partially de-twinned crystals, the relative intensity changes. In FeSe, this approach enabled the identification of which resonance corresponds to which domain [23].

The spin-lattice-relaxation rate varies strongly with uniaxial strain in BaFe₂As₂, as shown in Figure 9 [15]. Strain enhances the relaxation rate and the effect is anisotropic: the enhancement is stronger for in-plane fields than for out-of-plane fields. The enhancement can be understood qualitatively by realizing that strain enhances the antiferromagnetic ordering temperature, T _N , thus for T > T _N , the correlation length should grow with increasing strain. The correlation time τ _c grows with increasing correlation length, thus enhancing T 1 − 1 as illustrated in Figure 6.

A more complete picture of the nematicity and the antiferromagnetism can be obtained by analyzing the reciprocal-space structure of the spin susceptibility tensor, χ _αβ (q). In the disordered tetragonal phase, χ _xx (q) = χ _yy (q), but in the presence of nematicity these components are no longer equal. In the ordered phase, the susceptibility diverges at the antiferromagnetic ordering wavevectors, Q ₁ = (π, 0) in one domain, and Q ₂ = (0, π) in the other domain. A nematic order parameter can defined as ϕ α β = χ α α − 1 ( Q 2 ) − χ β β − 1 ( Q 1 ) [15]. By symmetry, the only non-vanishing components are ϕ _xy , ϕ _yx and ϕ _zz . In the ordered state, only ϕ _xy condenses, which corresponds to the spins oriented parallel to the antiferromagnetic ordering vector, as shown in the upper row of Figure 2. ϕ _zz order is illustrated in the second row of Figure 2. In this case, the spins are oriented along the c direction, and the hyperfine field is in-plane. Although this channel does not condense, ϕ _zz fluctuations are present and contribute to T 1 − 1 . ϕ _yx is illustrated in the bottom row, where the spins are oriented in-plane and are perpendicular to the antiferromagnetic ordering vector. For ϕ _yx order the hyperfine field vanishes, so T 1 − 1 is insensitive to these types of fluctuations.

The anisotropy of T 1 − 1 contains key information about these two nematic channels, ϕ _xy and ϕ _zz . The fluctuations h _α (τ) of the hyperfine field that contribute to the spin-lattice-relaxation rate as described in Eq. 10 are perpendicular to the applied field, H ₀. If H ₀ ∥c, then the in-plane fluctuations driven by ϕ _zz fluctuations will dominate T 1 − 1 . If H ₀ ⊥c, then both ϕ _zz and ϕ _xy fluctuations will contribute. As seen in Figure 9, T 1 − 1 measured for H ₀ ⊥c has a much stronger variation with applied strain. This suggests that ϕ _zz fluctuations are enhanced under strain more than ϕ _xy fluctuations. In fact, a detailed analysis of the nematic susceptibilities, χ x y n e m = ∂ ϕ x y / ∂ ε B 2 g and χ z z n e m = ∂ ϕ z z / ∂ ε B 2 g indicates that the latter is largest, despite the fact that ϕ _xy condenses at zero strain [15]. These results imply that for sufficient strain, the system would order with the spins oriented along the c-axis. This interpretation was recently confirmed by polarized inelastic neutron scattering results [54].

Measurements of the doped pnictides under strain are challenging because the doping introduces both static and dynamic inhomogeneities. The quadrupolar linewidths of the satellites broaden significantly with doping, thus it is more difficult to discern static changes in the resonances due to strain. Measurements of the spin-lattice-relaxation rate are more straightforward and do not require sharp resonances, but doping introduces a distribution of local hyperfine fields and autocorrelation times, τ _c , so that there is no longer a homogeneous T 1 − 1 value throughout the bulk of the materials [16].

In principle, an external strain field of sufficient magnitude might overcome the local strain fields and give rise to more homogeneous relaxation. This hypothesis was tested in Ba(Fe_1−x Co _x )₂As₂ with x = 0.048 under tensile strain [52]. Although the magnitude of T 1 − 1 changed in response to the strain field, no changes were observed to the stretching exponent. This observation suggests that the intrinsic strains around the dopants exceed the homogeneous external strain, so that P ( W 1 ) remains unchanged under strain. In fact, the lattice strain around Co dopants in Ba(Fe_1−x Co _x )₂As₂ can reach several percent [79]. The applied strain was only a fraction of a percent, thus it is not surprising that no change was observed in the width of the distribution.

Although the width of the P ( W 1 ) does not change with externally applied strain, the median of the distribution does, and in a manner that depends on the strain history. Figure 14 illustrates a protocol to investigate such hysteretic effects. A crystal is strain-cooled either under compressive (or tensile) strain from room temperature down to a base temperature. The spin lattice relaxation is then measured for a series of different levels of increasing (or decreasing) strain at constant temperature. The sample is then returned to room temperature (where the relaxation is homogeneous and there is no glassy behavior), and then the process is repeated under the opposite direction. For a fully linear response, the spin lattice relaxation as a function of strain should be independent of whether it was tensile-cooled or compressive-cooled.

Figure 15 shows the results of such a study on Ba(Fe_1−x Co _x )₂As₂ for four different base temperatures for x = 0.059. As expected, β decreases as temperature is reduced, but remains independent of the applied strain, and is the same for both tensile- and compressive-cooling. On the other hand, the median of the distribution, T 1 − 1 , is strain-dependent and is different under the two strain conditions at low temperature.

Note that the horizontal axes in Figure 15 are displacement, not strain. Strain is defined as ɛ = (d − d ₀)/L, where d is the displacement of the capacitative dilatometer, d ₀ is the equilibrium unstrained displacement, and L is the sample length. d ₀ can be determined at room temperature, but thermal contractions of both the sample and the strain device make it difficult to determine d ₀ under cryogenic conditions. It is best to utilize an independent measure of strain, such as the EFG tensor discussed in Section 3.1, to find the displacement corresponding to ɛ = 0. Such a calibration was not possible for the T 1 − 1 measurements in Figure 15. It is also important to note that the level of strain during the cooling part of the protocol in Figure 14 is likely temperature-dependent, rather than constant as illustrated. Again, this is due to differential thermal contractions between the sample and the strain device. The displacement measured by the dilatometer is the quantity that is held constant.

Despite these complications, these measurements reveal that P ( W 1 ) , and presumably the distribution of local nematic domains, depends on the past history of applied strain. Hysteresis is commonly observed in ordered ferromagnets, in which the static bulk magnetization is hysteretic with the applied field. If the applied field exceeds the coercive field, the magnetization will saturate, and when the applied field is reduced to zero, a remnant magnetization remains. Similar behavior should be expected in an ordered nematic without quenched disorder, in which nematic domains can reorient in an external strain field [80]. The NMR measurements described here, on the other hand, are in a nematic glass phase with quenched disorder, and probe the distribution of local correlation times of different fluctuating domains. Hysteresis has been studied in magnetic materials with quenched disorder in both the weak-disorder limit in the paramagnetic phase, as well in the strong-disorder limit in the spin glass phase [81]. The particular domain pattern is likely determined by the combination of local strains and the external strain field applied during the cool down. The fluctuation time of the domains is a complicated function of these strain fields as well as the behavior of their neighboring nematic domains. As the external strain field changes, the local correlation times change, but in such a manner that the median of the distribution of correlation times changes, but not the standard deviation. If the system were homogenous and not glassy, then the nematic domains would be mobile and there would be a single average correlation time. The presence of a distribution of correlation times that is altered by external strain provides evidence for the formation of a nematic glass, in which some domains exhibit very large correlation times, whereas others have short correlation times and may be reoriented by small external strain fields.

A quantitative theory to connect the measured distribution of the spin-lattice-relaxation rate, P ( W 1 ) , to the nematic domains and their dynamics in the presence of external strain, particularly in the context of the iron pnictides, is highly desirable. Although no such theory yet exists, there are important models that capture elements of the glassy behavior in these materials. Nematic ordering in the presence of random strain fields can be mapped onto the random field Ising model (RFIM) [80, 82, 83]. Detailed numerical simulations have uncovered a broad distribution of local correlation times, which agrees qualitatively with the observations of stretched relaxation in NMR experiments [16].

On the other hand, the iron pnictides exhibit both antiferromgnetism and nematicity, and there are two types of spin order for each of two possible nematic domains, as illustrated in Figure 2. Recently, a random Baxter-field model was proposed as an alternative to capture both the nematic and antiferromagnetic fluctuations in the presence of quenched random strain fields [84]. In this case, there are two Ising variables at each lattice site associated with the magnetic and the nematic degree of freedom. These two variables are coupled due to the next nearest neighbor coupling, J ₂, in the J ₁ − J ₂ model for the iron pnictides [4]. Random strain fields couple only to the nematic degree of freedom, but the magnetic variable is also affected and consequently the system breaks up into domains of all four configurations illustrated in the top row of Figure 2.