Dynamic spin-lattice coupling and nematic fluctuations in NaFeAs

We use inelastic neutron scattering to study acoustic phonons and spin excitations in single crystals of NaFeAs, a parent compound of iron pnictide superconductors. NaFeAs exhibits a tetragonal-to-orthorhombic structural transition at $T_s\approx 58$ K and a collinear antiferromagnetic (AF) order at $T_N\approx 45$ K. While longitudinal and out-of-plane transverse acoustic phonons behave as expected, the in-plane transverse acoustic phonons reveal considerable softening on cooling to $T_s$, and then harden on approaching $T_N$ before saturating below $T_N$. In addition, we find that spin-spin correlation lengths of low-energy magnetic excitations within the FeAs layer and along the $c$-axis increase dramatically below $T_s$, and show weak anomaly across $T_N$. These results suggest that the electronic nematic phase present in the paramagnetic tetragonal phase is closely associated with dynamic spin-lattice coupling, possibly arising from the one-phonon-two-magnon mechanism.


I. INTRODUCTION
Spin waves (excitations) and phonons are two fundamental quasiparticles in a solid describing propagating disturbance of the ordered magnetic moment and lattice vibrations, respectively [1]. They are the by-products of linearized theories that ignore all the terms that are higher order than the quadratic term and that neglect quasiparticlequasiparticle interactions [2], unlikely interacting with each other [3]. However, in systems with multiple coupled degrees of freedom, the interaction affects the self-energy of the quasiparticles and changes their intrinsic dispersion relationship. Therefore, discovering and understanding how the interaction and coupling change the spectra of spin excitations and phonons and further influence the electronic properties of solids is one of the most important themes in modern condensed matter physics.
In general, spin-lattice (magnon-phonon) coupling can modify the spin excitations in two different ways. First, the static lattice distortion coupled with magnetic order, as seen in some parent compounds of iron-based superconductors [4,5], can dramatically change the effective magnetic exchange anisotropy in a near square lattice of the antiferromagnetic (AF) ordered phase [6]. Second, the dynamic lattice vibrations (phonons) interacting with time-dependent spin excitations may create energy gaps in the magnon dispersion at the nominal intersections of the magnon and phonon modes [1], as seen in the insulating noncollinear antiferromagnet ðY; LuÞMnO 3 [7]. Similarly, spin waves can hybridize with optical phonon modes near the zone boundary and form broadened magnon-phonon excitations in metallic ferromagnetic manganites [8]. In both cases, magnon-phonon interactions only form hybridized excitations at the magnon-phonon intersection points due to the repulsive magnon-phonon dispersion curves.
In this paper, we report a different form of magnon-phonon interaction in NaFeAs, which is a parent compound of an iron-pnictide superconductor and exhibits a tetragonal-toorthorhombic structural transition at T s ≈ 58 K followed by a collinear AF order at T N ≈ 45 K [ Fig. 1(a)] [9][10][11][12]. For temperatures between T N and T s , NaFeAs has an electronic nematic phase, a form of electronic order that breaks the rotational symmetries without changing the translational symmetry of the underlying lattice [13][14][15][16][17][18][19][20][21][22][23][24][25][26], which may play an important role in high-temperature superconductivity [27,28]. Instead of the usual magnon-phonon interaction where the magnon dispersion would intersect with the phonon dispersion in reciprocal space, we find evidence for a magnon-phonon interaction in NaFeAs, where the low-energy acoustic phonons near the nuclear-zone-center Γ point in NaFeAs can interact with the low-energy spin excitations near the AF zone-center M point at Q AF in reciprocal space without direct phonon-spin excitation crossing [Figs. 1(b) and 1(c), Figs. [2][3][4][5]. Since both the in-plane transverse acoustic (IPTA) phonons [ Fig. 1(b)] and spin excitations show dramatic anomaly across T s , our results indicate that the magnon-phonon interaction may be responsible for the low-temperature electronic nematic phase in NaFeAs [13][14][15][16][17][18][19][20][21][22][23][24][25][26]. Although not conclusively established, one possible microscopic interpretation of these results is the one-phonon-two-magnon scattering mechanism [29,30]. By considering dynamical spin exchange coupling tuned by local lattice vibration, we calculate the effect of magnon-phonon coupling in a one-dimensional (1D) spin chain. Further phonon dispersion calculations using density  Red spots represent Fe atoms, and the tetragonal (orthorhombic) unit cell is shown by the green (blue) shaded square. The middle and right panels are the two orthorhombic domains. The yellow and red areas are the tetragonal unit cells with opposite shear distortions. The arrows in the middle panel indicate the spin direction in the AF ordered phase. The orthorhombic lattice parameters are a o ¼ 5.589 Å and b o ¼ 5.569 Å. (b) A particular IPTA phonon mode with a momentum propagating horizontally to the right but with Fe atoms oscillating vertically. Here, λ is the size of the probing domain, and it is inversely related to q. (c) The three different equivalent positions in reciprocal space used in our neutron scattering measurements in order to determine all three acoustic phonon branches. The purple arrows with q are the momentum of the measured phonons along the ½H; 0; 0 direction, and the double-headed arrows represent the phonon polarizations along the a (red, LA), b (blue, IPTA), and c (green, OPTA) axes. (d) The corresponding dispersions of the three phonon modes at 100 K estimated from data in Ref. [37]. The inset shows an expanded view of the low-energy part of the dispersion. functional theory (DFT) on NaFeAs demonstrate that softening of the IPTA mode occurs at q → 0 in the presence of stripe-type AF spin-spin correlations, consistent with our observation. Therefore, instead of a pure spin [14][15][16][17] or lattice (orbital) [31][32][33][34][35][36] degree of freedom, dynamic spinlattice coupling is important for the electronic nematic phase and anisotropic electronic properties in the paramagnetic state of NaFeAs.

A. Neutron scattering
We chose to study NaFeAs because transport [18,19], angle-resolved photoemission spectroscopy (ARPES) [20,21], scanning tunneling microscopy [22,23], and neutron scattering experiments [24][25][26] provide ample evidence for the presence of an electronic nematic phase in T N ≤ T ≤ T s . By assuming weak spin-orbit and spinlattice coupling, most theoretical understandings of the nematic phase involve pure "spin" or "orbital" degrees of freedom. In the first case, the electronic nematic phase originates from spin fluctuations that couple quadratically to the nematicity but have no direct coupling to the lattice [14][15][16]. Here, key predictions of the theory, such as the enhanced spin excitation intensity and correlation length at the AF ordering wave vector Q AF in the paramagnetic orthorhombic phase of iron pnictides [38,39], have been confirmed by inelastic neutron scattering experiments [24][25][26]. Alternatively, ferro-orbital order or orbital fluctuations [31][32][33][34][35][36], seen in ARPES measurements as scans along the cut-3 direction as shown in Fig. 3(a) at 45 K, 50 K, 55 K, and 60 K. Background scans are obtained by measuring at q ¼ ð0.4; 0.4; LÞ and ð0.6; 0.6; LÞ and are subtracted from the raw data. The scan at 60 K is featureless. On cooling to 55 K, two weak peaks show up around L ¼ AE0.5 r.l.u. Upon further cooling, these peaks become prominent at T ¼ 45 K. They disappear below T N due to the opening of a spin gap. The total magnetic scattering intensity increases on cooling. The solid lines are fits with a periodic Lorentz function [37]. (e)-(g) Temperature dependence of the linewidths and correlation lengths of the low-energy spin fluctuations along the c-axis direction [cut 3 in Fig. 3(a)] at 2 meV, 3 meV, and 5 meV. (h) Temperature dependence of the integrated intensity along the cut-3 direction at different energies [37]. a pronounced energy splitting of bands with d xz and d yz orbital characters above T N [20,21], may induce the electronic nematicity through the direct ferro-orbital order (fluctuations) and orthorhombic lattice distortion (or vibration) coupling [40,41].
The difficulty in sorting out the microscopic origin of the nematic phase in iron pnictides lies in the fact that the relevant degrees of freedom-such as spin, lattice, and orbital-are closely entangled [40,41]. However, DFT calculations of iron pnictides [42][43][44] show a remarkable sensitivity of superconductivity and magnetism to the lattice, particularly the iron pnictogen height, thus suggesting a strong dynamic spin-lattice coupling in this system. Measurements and calculations on the phonon dispersions in CaFe 2 As 2 and BaFe 2 As 2 suggest that magnetism must be taken into account when calculating the phonon frequencies, and phonons are strongly coupled not to the static AF order but to spin fluctuations [45][46][47]. In addition, inelastic x-ray and neutron scattering measurements on AF ordered BaFe 2 As 2 and SrFe 2 As 2 , where T N ≈ T s , show a softening of the IPTA phonon on cooling to T s [Figs. 1(a) and 1(b)], followed by a dramatic hardening below T N =T s [48,49]. For electron-doped superconducting BaðFe 1−x Co x Þ 2 As 2 with T c < T N < T s , the softening of the IPTA phonon stops at T s on cooling, and the hardening of the mode passes through T N unabated before softening again below T c [50]. On the other hand, a sharp enhancement in the in-plane spin-spin correlation length was found below T s at Q AF in LaFeAsO and BaðFe 1−x Co x Þ 2 As 2 [25], suggesting a strong feedback effect of the structural transition and nematic order on the low-energy spin fluctuations. Although these results suggest a correlation between the evolution of the transverse acoustic phonon and low-energy spin fluctuations, the measurements were carried out on different samples by different groups, thus making it difficult to establish a direct connection between lattice vibrations and spin fluctuations.
To resolve this problem, we use inelastic neutron scattering to study phonons and spin excitations in single crystals of NaFeAs, which is the parent compound of the NaFe 1−x Co x As family of iron-based superconductors [51]. Our measurements were carried out on the C5 triple-axis spectrometer at the Canadian Neutron Beam Center, Chalk River Laboratories, with fixed final neutron energies of E f ¼ 8.04 meV and 14.56 meV for phonon measurements and E f ¼ 13.7 meV for spin fluctuation measurements. Some scans of the IPTA phonon were carried on the HB-3 tripleaxis spectrometer at the High Flux Isotope Reactor (HFIR), Oak Ridge National Laboratory, with E f ¼ 14.7 meV. The wave-vector transfer Q in three-dimensional (3D) reciprocal where H, K, and L are Miller indices. In the paramagnetic tetragonal state, we have a Å, and AF spin fluctuations occur at Q AF ≈ ð0.5; 0.5Þ in the tetragonal notation [4]. The phonon measurements were made on a single piece of NaFeAs single crystal weighing about 2.5 g with a mosaic less than 1°aligned in the ½H; K; 0 and ½H; 0; L scattering planes [ Fig. 1(c)]. The spin fluctuation measurements were carried out using a well-aligned pack of single crystals (∼10 g) with a mosaic less than 3°a ligned in the ½H; H; L scattering plane [ Fig. 3(a)].
To measure IPTA, out-of-plane transverse acoustic (OPTA) and longitudinal acoustic (LA) phonon modes in NaFeAs, we probed wave vectors Q ¼ ðq; 2; 0Þ, ðq; 0; 5Þ, and ð2 þ q; 0; 0Þ, respectively, where q is the reduced momentum transfer away from the Γ point in reciprocal space [ Fig. 1(c)] [48,49]. In the long wavelength (small jqj) limit, one would expect linear dispersions for acoustic phonons E ¼ vq, where E is the energy of the mode and v is the sound velocity. Figure 1(d) shows our measured dispersions of LA, OPTA, and IPTA acoustic phonon modes at 100 K [37]. While LA and OPTA phonons behave as expected and pass through the zero energy at q ¼ 0, the IPTA mode, which corresponds to the q ¼ 0 shear modulus mode C 66 seen in the ultra sound spectroscopy measurements [52][53][54], has a lower sound velocity as q → 0 [dashed lines in Fig. 2(a)]. The solid lines are linear fits to the large wavevector data. Figure 2(a) shows dispersions of the IPTA phonons determined from a series of constant-Q scans at different temperatures [37]. While the high-q (>0.2 r.l.u.) part of the dispersion hardens slightly with decreasing temperature because of a stiffer lattice at low temperature (but showed no anomaly across T s and T N ), the low-q (<0.2 r.l.u.) part of the dispersion is strongly temperature dependent [ Fig. 2(b)] [37]. This is seen in the temperature dependence of the constant-Q scans at q ¼ 0.05 r.l.u. [Fig. 2(c)], where the peak center shifts from 1.25 meV at 100 K to 0.9 meV at 60 K near T s . On further cooling to below T s , the phonon hardens with decreasing temperature. The temperature dependence of the phonon energy at different q is summarized in Figs. 2(d) and 2(e).
In general, acoustic phonons with different polarizations should exhibit similar linear dispersion relationship at similar temperatures. When phonons couple with electrons or other elementary excitations, their peak positions, widths, and intensities might change [55]. While LA and OPTA phonon dispersions in Fig. 1 Fig. 1(d)]. Since the IPTA phonon at q ¼ 0.2 r.l.u. is weakly temperature dependent [Figs. 2(a) and 2(b)], its temperature-dependent effect for q < 0.2 r.l.u. will be absorbed into the parameter b [ Fig. 2(a)] [37]. Figure 2(f) shows temperature dependence of b for the LA, OPTA, and IPTA phonons obtained by fitting with E ¼ vq þ b. While the LA and OPTA phonons behave as expected with b ¼ 0, the values of b for the IPTA mode are negative in both the AF ordered and paramagnetic state, thus suggesting the presence of the IPTA phonon-elementary excitation (either spin or orbital) coupling that is independent of the magnetic ordering. Since such coupling is absent for the LA and OPTA phonon modes, we conclude that the nonlinear dispersion curves seen in the IPTA phonon mode at small q and their strong temperature dependence around T s are due to the IPTA phonon and elementary (spin or orbital) excitation coupling, with the absolute value b reflecting the coupling strength.
In addition to phonons, spin fluctuations are also important elementary excitations in iron pnictides [4,13]. Distinct from acoustic phonons which are observed near nuclear Bragg peak positions such as (2, 0, 0), low-energy spin fluctuations in NaFeAs occur near the Q AF ¼ ð0.5; 0.5; 0.5Þ AF zone center in reciprocal space [ Fig. 3(a)] [4]. Figure 3 summarizes temperature dependence of the in-plane spin-spin correlations for spin fluctuations with energies of E ¼ 3 meV, 5 meV, and 8 meV along the cut-1 and cut-2 directions as shown in Fig. 3(a). Figures 3(b)-3(e) show temperature dependence of the linewidth and the corresponding spin correlation length ξ, respectively. A sharp reduction in the peak linewidth and a dramatic increase in the spinspin correlation lengths are seen below T s [37]. They are accompanied by the increase of the intensities at the peak centers, reflecting a redistribution of the magnetic spectral weight in the nematic phase. While these results are somewhat different from a neutron scattering work on NaFeAs [56], which finds no dramatic anomaly across T s , they are consistent with previous work on LaFeAsO and BaðFe 0.953 Co 0.047 Þ 2 As 2 , which was interpreted as resulting from the long-range nematic order below T s [25]. We note that a phonon at q ¼ 0.05 probes coherent lattice vibrations of about 20 unit cells, giving a dynamic length scale of about 40-80 Å. These values are similar to the observed spin-spin correlation lengths above T s , suggesting that correlated spin domains might be limited by sizes of the dynamic structural domains.
To further understand spin excitations across T s , we carefully studied the temperature dependence of spinspin correlations along the c axis. Figures 4(a)-4(d) are the background-subtracted wave-vector scans along the ½0.5; 0.5; L direction at a spin excitation energy of E ¼ 2 meV and different temperatures. At T ¼ 60 K (>T s ), the scattering is essentially featureless and decreases with increasing L due to the iron magnetic form factor, suggesting that magnetic excitations are two-dimensional (2D) rods along the c axis in reciprocal space [ Fig. 4(d)]. On cooling to T ¼ 55 K (<T s ), Two broad peaks emerge around L ¼ AE0.5 r.l.u., indicating an increased c-axis spinspin correlation length [ Fig. 4(c)]. Upon further cooling to T ¼ 50 K and 45 K (>T N ), these peaks become progressively better established [ Figs. 4(b) and 4(a)]. The solid lines of Figs. 4(a)-4(d) are fits to these peaks with a periodic Lorentzian function [37]. Figures 4(e)-4(g) plot the temperature dependence of the fitted width and the corresponding spin correlation length at E ¼ 2 meV, 3 meV, and 5 meV. Inspection of the figures reveals a clear change of the c-axis spin correlation length below T s , suggesting that a spin fluctuation crossover from 2D to 3D occurs below T s and above T N . Figure 4(h) shows the temperature dependence of the integrated intensities over the range from L ¼ −1 to L ¼ 1, again revealing a clear increase in magnetic scattering below T s . From the results described in Figs. 3 and 4, we conclude that the low-energy spin fluctuations in NaFeAs respond dramatically to the structural transition and change its spectral distribution below T s .

B. Theoretical considerations and density functional theory calculation
Given the simultaneous occurrence of the dramatic increasing spin correlation length at Q AF and a hardening of the IPTA phonon at the Γ point below T s (Figs. 2-4), it is interesting to ask if the low-energy spin fluctuations are coupled to the IPTA phonons. To address this question, we consider the differences between paramagnetic tetragonal and orthorhombic phases. In the paramagnetic tetragonal phase, iron spins are completely random, spending an equal amount of time along all directions. On cooling to the paramagnetic orthorhombic phase below T s , the direction of the spin alignment gets approximately fixed to the lowtemperature AF ordered collinear phase due to dynamic spin-lattice coupling; thus, within a structural domain, iron spins in the nearest neighbors are predominantly antiparallel along the orthorhombic a o axis and parallel along the b o axis. As a consequence, the magnetic correlation length changes abruptly across T s but does not become infinite due to the lack of long-range magnetic order. Regardless of a detailed microscopic spin-lattice coupling mechanism, this picture is valid and suggests a magnetic origin of the nematic phase [16].
To provide a possible microscopic description of the spin-lattice coupling, we consider the magnitude of the effective magnetic exchange couplings in NaFeAs. By fitting the spin-wave spectra throughout the Brillouin zone at 5 K using a Heisenberg Hamiltonian, we find that the effective magnetic exchange couplings along the orthorhombic a o and b o axes are SJ 1a ¼ 40 AE 0.8 meV and SJ 1b ¼ 16 AE 0.6 meV, respectively, where S ≈ 1 is the spin of the system [57]. Since the distances between the nearest Fe-Fe neighbors in the AF orthorhombic phase are a o =2 ¼ 2.795 Å and b o =2 ¼ 2.785 Å [10], the large changes in the magnetic exchange in the AF ordered phase may arise from tiny differences in the in-plane Fe-Fe distances within the Heisenberg Hamiltonian [ Fig. 5(a)] [58]. In the paramagnetic tetragonal state, the shear modulations of the IPTA phonon mode [Figs. 1(a) and 1(b)] may induce large changes in magnetic exchange through dynamic spinlattice coupling [ Fig. 5(b)]. In previous ultrafast optical spectroscopy measurements on BaðFe 1−x Co x Þ 2 As 2 , it was found that coherent optical phonon motions of As atoms on ultrafast timescales induced by a femtosecond optical pulse can produce transient dynamic AF and/or nematic order at temperatures above T s through possible modifications of the Fermi surfaces and dynamic spin-lattice coupling [60,61]. For conventional AF insulators such as ðY; LuÞMnO 3 , the effect of dynamic spin-lattice coupling is to create energy gaps in the magnon dispersion at the nominal intersections of the magnon and phonon modes [7]. However, this cannot explain our observation in NaFeAs since the low-energy acoustic phonons and spin fluctuations are located at different wave vectors in reciprocal space [ Fig. 3(a)].
To solve this problem, we consider the one-phonon-twomagnon scattering mechanism [29,30]. Although such a mechanism is no different from the usual magnetoelastic coupling or magnon-phonon coupling, it provides a simple pictorial way to visualize such a coupling. Theoretically, it has been found that spin-lattice coupling due to the onephonon-two-magnon scattering process in an antiferromagnet can affect the spin-wave broadening or damping, dispersion, and intensity [62]. Since spin fluctuations in iron pnictides exhibit spin-wave-like features above T N up to very high temperatures [59], such magnon-phonon coupling must also exist in the paramagnetic tetragonal phase. Figure 5(b) schematically illustrates how such a dynamic spin-lattice scattering mechanism might work. In the paramagnetic tetragonal phase with strong AF spin fluctuations, the neighboring spins are aligned antiparallel on a timescale determined by the energy of the fluctuations, and the average moment of the system at elastic position (E ¼ 0 or infinite time) is still zero. When an IPTA phonon with momentum q and energy E 0 is emitted near the Γ point [ Fig. 1(b)], the Fe-Fe distance stretches, and thus the exchange coupling increases from J 0 to J 0 þ ΔJ [ Fig. 5(b)] [37]. This represents an energy transfer from the lattice to the spin system, annihilating one phonon and producing spin fluctuations, and vice versa.
To test the possible spin-lattice coupling, we built a toy model of a 1D spin chain in which the spin exchange coupling is tuned by the local lattice vibration induced by a passing phonon (see Ref. [37], Sec. IV). Similar to the recent calculations on the effect of one-phonon-twomagnon scattering mechanism to the indirect K-edge resonant inelastic x-ray scattering spectrum of a 2D Heisenberg antiferromagnet [62], we find that spin waves in a 1D chain are also affected by the lattice distortion. In particular, the periodic lattice distortion affects spin waves, which in turn leads to the phonon softening at small q, similar to our observation in NaFeAs (Fig. S21). Since analytical analysis in the 2D case is difficult, we perform DFT calculations to compare the phonon dispersions of NaFeAs with and without a magnetic moment on the Fe site. The open triangles, squares, and circles in Fig. 5(f) are DFT-calculated LA, OPTA, and IPTA phonon dispersions with a Fe magnetic moment, respectively. The solid lines in the same figure show identical calculations without magnetism. While the main effect of magnetism is to soften the OPTA and IPTA phonons near the zone boundary (q ¼ 0.5), consistent with earlier work [45][46][47], its effect on the lowenergy OPTA and IPTA phonons can be seen by comparing phonon energies with magnetism (E M ) and without magnetism (E NM ). Figure 5(g) shows wave-vector dependence of E M − E NM , revealing a clear difference between OPTA and IPTA phonons below q ¼ 0.075. These calculations suggest that the presence of a magnetic moment on Fe has a much larger effect on the low-energy IPTA phonons, consistent with our observation of the phonon softening determined by the intercept b [ Fig. 2(b)]. To test if the small orthorhombic lattice distortion in the nematic phase of NaFaAs can affect acoustic phonons without magnetism, we show in Fig. 5(h) the calculated phonon dispersions in orthorhombic and tetragonal states in the nonmagnetic state. Within the errors of our calculation, we find no evidence of phonon softening due to lattice orthorhombicity. From these results, we conclude that the presence of magnetism and spin-spin correlations selectively influences the dispersion of the IPTA phonon at q → 0.

III. DISCUSSION AND CONCLUSION
Assuming that the spin-lattice coupling in NaFeAs is due to the one-phonon-two-magnon mechanism, the scattering process must satisfy the energy and momentum conservation [30,63]. Two magnons with momenta k and k 0 and energies E and E 0 with a relationship q ¼ k þ k 0 and E 0 ¼ E þ E 0 (positive E means generating, and negative E means annihilating) must be introduced, as shown in Fig. 5(c). The scattering amplitude is constrained by the momentum triangle, and it depends on the detail of the spin fluctuation dispersion [30]. In the paramagnetic tetragonal state, the magnetic scattering is a rod centered around Q AF ¼ ð0.5; 0.5Þ in the FeAs plane [Figs. 3(a), 4(d), and 5(d)]. This implies a large phase space for the dynamic spinlattice coupling to occur. However, when a 2D-to-3D crossover in spin fluctuations occurs below T s , the spinlattice coupling phase space is greatly reduced, as shown in Fig. 5(e). Since the IPTA phonon mode dynamically changes the nearest-neighbor Fe-Fe distance, which is directly coupled to the nearest-neighbor magnetic exchange coupling, the reduced scattering phase space below T s explains the recovery of the phonon softening (phonon hardening) at low temperatures in NaFeAs. Therefore, the observed IPTA phonon hardening below T s is a sign of reduced phonon-magnon coupling induced by the 2D-to-3D spin fluctuation crossover. In this picture, the narrowing of the low-energy in-plane spin fluctuations below T s is associated with the orthorhombic lattice distortion and nematic order, which increases the effective exchange coupling along a and therefore the spinwave velocity. Since finite magnetic exchange coupling along the c axis is required to establish the eventual 3D AF long-range order, the presence of the quasi-2D spin fluctuations in the paramagnetic tetragonal phase and the resulting large dynamic spin-lattice coupling are directly associated with the observed nematic fluctuations in a variety of experiments [19][20][21][22][23][24][25][26]. The energy scale of the phonon softening of about 1 meV [64] is consistent with the slow recovery of the nematic order on a timescale of several or decades of ps observed by time-resolved polarimetry [61], further confirming its close relationship with the nematic fluctuations.
Although the one-phonon-two-magnon mechanism described above may provide a pictorial way to explain the experimentally observed phonon and spin excitation anomaly across T s in NaFeAs, we recognize that any effect of the dynamic spin-lattice coupling can be described as such a process, and we have not established a quantitative perturbative theory based on such a process. Regardless of the microscopic origin of the observed phonon and spin excitation anomaly across T s , our experimental results reveal a clear coupling between low-energy acoustic phonon and spin excitations at T s . Therefore, instead of the pure spin or lattice (with associated orbital ordering) as the driving force for the electronic nematic phase in iron pnictides, a strong dynamic spin-lattice-orbital coupling in the paramagnetic tetragonal phase is responsible for the observed electronic nematic order and anisotropic behavior.