Strain tuning and anisotropic spin correlations in iron-based systems

TL;DR

This paper investigates how strain influences nematic and magnetic phase transitions in iron-based superconductors, revealing that strain sensitivity arises from spin fluctuations and anisotropic magnetic excitations, with implications for understanding related vestigial phases.

Contribution

It provides a theoretical framework explaining strain effects on nematicity and magnetic transitions, emphasizing the role of spin fluctuations and anisotropic magnetic spectra in iron-based systems.

Findings

Strain sensitivity is primarily a fluctuation effect.

Antisymmetric strain strongly affects magnetic excitations.

The formalism can be extended to other magnetic vestigial phases.

Abstract

Nematic order in the iron-based superconductors is closely tied to a lattice distortion and a structural transition from tetragonal to orthorhombic symmetry. External stress of the appropriate symmetry acts as a conjugate field of the nematic order parameter and can be utilized to detwin nematic domains but also smears an otherwise sharp nematic transition. On the other hand, applying stress in proper symmetry channels allows one to tune the nematic phase transition. Recent experiments analyzed the stress-induced changes of the nematic and magnetic phase transition temperature. Here we show that the observed trends can be understood in terms of spin-induced nematicity. The strain sensitivity is shown to be a fluctuation effect. The strong sensitivity to antisymmetric strain is a consequence of the anisotropic nature of the magnetic excitation spectrum. The formalism presented here can be naturally generalized to determine the strain-sensitivity of vestigial phases related to other magnetic states that have been observed in the iron-based systems, such as e.g.\ the spin-charge density wave and the spin-vortex crystals.

Equations

$$T_{\mathrm{nem}}=T_{\mathrm{nem}}^0(1-D_{B_{2g}}{\epsilon }_{B_{2g}}^2),$$

$$T_{\mathrm{nem}}=T_{\mathrm{nem}}^0(1+C_{A_{1g}}{\epsilon }_{A_{1g}}-D_{A_{1g}}{\epsilon }_{A_{1g}}^2).$$

$$H=J_1\sum _{⟨ij⟩}\mathbf{s}({\mathbf{r}}_i)\cdot \mathbf{s}({\mathbf{r}}_j)+J_2\sum _{⟨⟨ij⟩⟩}\mathbf{s}({\mathbf{r}}_i)\cdot \mathbf{s}({\mathbf{r}}_j).$$

$$J_2\to J_2^{\pm }\equiv J_2(1\pm \kappa {\epsilon }_{B_{2g}}),$$

$$H=\sum _{\mathbf{k}}J(\mathbf{k}){\mathbf{s}}_{\mathbf{k}}\cdot {\mathbf{s}}_{-\mathbf{k}},$$

$$J(\mathbf{k})$$

$$+2J_2^{+}\cos (k_x+k_y)+2J_2^{-}\cos (k_x-k_y).$$

$$E_{cl.}=-(J_2^{+}+J_2^{-})s^2=-2J_2{s}^2.$$

$$\mathbf{s}(\mathbf{r})$$

$$S$$

$$\displaystyle\qquad\qquad+2\varepsilon_{B_{\mathrm{2g}}}(\partial_{x}\boldsymbol{n}_{1}\!\cdot\!\partial_{y}\boldsymbol{n}_{1}+\partial_{x}\boldsymbol{n}_{2}\!\cdot\!\partial_{y}\boldsymbol{n}_{2})\Big{\}}.$$

$$\displaystyle\!S\!=\!\frac{1}{2\tilde{g}}\!\int\!\!d\tilde{\boldsymbol{r}}\,\Big{\{}\!\Big{[}\!\!\sum_{j=1,2}\!\!(\tilde{\nabla}\boldsymbol{n}_{a})^{2}\Big{]}\!+2\tilde{\alpha}(\partial_{\tilde{x}}\boldsymbol{n}_{1}\!\cdot\!\partial_{\tilde{x}}\boldsymbol{n}_{2}\!-\!\partial_{\tilde{y}}\boldsymbol{n}_{1}\!\cdot\!\partial_{\tilde{y}}\boldsymbol{n}_{2}),$$

$$J_2\to J_2(1-{\epsilon }_{B_{2g}}^2{)}^{1/2}<J_2.$$

$$\displaystyle\!\!\xi\sim a_{0}e^{2\pi/\tilde{z}\tilde{g}}\approx a_{0}\exp\Big{\{}\frac{2\pi}{zg}[1\!-\!(1+\alpha/4)\varepsilon_{B_{\mathrm{2g}}}^{2}/2)]\Big{\}},$$

$$\omega (\mathbf{k})=4s{J}_2\sqrt{A(\mathbf{k}{)}^2-B(\mathbf{k}{)}^2},$$

$$A(\mathbf{k})$$

$$B(\mathbf{k})$$

$$b(\mathbf{k})=\frac{\alpha }{2}\sum _{\sigma =\pm }(1+\sigma {\epsilon }_{B_{2g}})\cos (k_x+\sigma k_y).$$

$$\displaystyle\Delta F=T\sum_{\boldsymbol{k}}\ln\Big{[}\sinh\Big{(}\frac{\omega(\boldsymbol{k})}{2T}\Big{)}\Big{]},$$

$$\Delta E={\gamma }_Q\frac{J_{2}s}{2}{\alpha }^2[1+({\mathbf{n}}_1\cdot {\mathbf{n}}_2{)}^2],$$

$${\gamma }_Q={\gamma }_Q^{(0)}(1+\zeta {\epsilon }_{B_{1g}}^2).$$

$$\varphi =⟨{\mathbf{n}}_1\cdot {\mathbf{n}}_2⟩,$$

$$\displaystyle\Big{[}\frac{\xi(T_{\mathrm{nem}})}{a_{0}}\Big{]}^{2}\Delta E$$

$$\mathbf{s}(\mathbf{R})={\mathbf{m}}_x(\mathbf{R})e^{i{\mathbf{Q}}_x\cdot \mathbf{R}}+{\mathbf{m}}_y(\mathbf{R})e^{i{\mathbf{Q}}_y\cdot \mathbf{R}},$$

$$S=\sum _l\int d\mathbf{r}\mathcal{L}(\mathbf{m},{\partial }_{\beta }\mathbf{m})$$

$${\mathcal{L}}_{\parallel }$$

$$-\frac{g}{4}({\mathbf{m}}_l{\tau }_3{\mathbf{m}}_l{)}^2+\frac{v}{4}({\mathbf{m}}_l{\tau }_1{\mathbf{m}}_l{)}^2+\frac{1}{2}{\partial }_{\beta }{\mathbf{m}}_l{\tau }_0{\partial }_{\beta }{\mathbf{m}}_l$$

$$+\frac{\alpha }{2}({\partial }_x{\mathbf{m}}_l{\tau }_3{\partial }_x{\mathbf{m}}_l-{\partial }_y{\mathbf{m}}_l{\tau }_3{\partial }_y{\mathbf{m}}_l)$$

$${\mathcal{L}}_{\perp }=q_0^2{\mathbf{m}}_l{\tau }_0{\mathbf{m}}_{l+1}.$$

$$S_{\epsilon }=S-\sum _l\int d\mathbf{r}\sum _{\alpha \beta }{\epsilon }_{\alpha \beta }{T}_{\alpha \beta }.$$

$$T_{\alpha \beta }=\frac{\partial \mathcal{L}}{\partial ({\partial }_{\alpha }\mathbf{m})}\cdot {\partial }_{\beta }\mathbf{m}-{\delta }_{\alpha \beta }\mathcal{L}.$$

$${\epsilon }_{A_{1g}}$$

$${\epsilon }_{B_{1g}}$$

Full text

I Introduction

Nematicity is a well-established state of electronic order in Fe-based superconductors [Fang et al., 2008; Xu et al., 2008; Si and Abrahams, 2008; Fernandes et al., 2010; Cano et al., 2010; Chuang et al., 2010; Chu et al., 2010; Lv and Phillips, 2011; Fernandes et al., 2012; Fernandes and Schmalian, 2012; Chu et al., 2012; Liang et al., 2013; Stanev and Littlewood, 2013; Fernandes R. M. et al., 2014; Böhmer et al., 2014]. Nematic fluctuations were identified [Fernandes et al., 2010] via a significant softening of the elastic shear modulus $C_{66}$, a behavior that was interpreted in terms of spin-induced nematicity [Fang et al., 2008; Xu et al., 2008; Si and Abrahams, 2008; Fernandes et al., 2010, 2012; Fernandes and Schmalian, 2012; Fernandes R. M. et al., 2014]. The latter can be understood as partially melted striped spin density-waves that break rotational order, without breaking time-reversal symmetry. Thus, within the spin-induced scenario, nematic order is caused by fluctuations of $C_{2}$-symmetric single-$\boldsymbol{Q}$ spin order. The strong coupling between elastic and nematic degrees of freedom has been exploited in elastoresistivity measurements [Chu et al., 2012] and three-point bending measurements [Böhmer et al., 2014]. It is rooted in the fact that the Ising-nematic order parameter $\phi$ couples bi-linearly to the corresponding strain $\varepsilon_{B_{\mathrm{1g}}}$. Notice, by ’Ising’-order we mean a single-component scalar order-parameter that falls into the universality class of Ising systems.

Whereas this work shall focus on the nematicity associated with fluctuations of a spin-density wave order, the nematic order may also arise due to orbital or charge order Yamase and Zeyher (2013); Onari and Kontani (2012). One may further speculate that nematicity is at work in the hidden order observed in iridates Zhao et al. (2015), or even at the origin of the pseudogap phase Fechner et al. (2016); Fradkin et al. (2010); Orth et al. (2019). A transition of a parent disordered phase into a phase of electronic or spin nematicity that breaks a crystal symmetry is associated with a softening of elastic moduli and results in anisotropic transport responses. On the other hand a standard phonon-driven structural transition has almost identical anisotropic responses and may also facilitate the subsequent appearance of a magnetic or charge order. Resolving the problem of which mechanism is driving the transition has conclusively been answered in the pnictide materials, primarily thanks to elastoresistivity measurements Chu et al. (2010, 2012).

Recently, Ikeda et al. [Ikeda et al., 2018] investigated the impact of strain on the nematic and magnetic phase transitions of $\mathrm{Co}$-doped $\mathrm{Ba}\mathrm{Fe}_{2}\mathrm{As}_{2}$ in different symmetry channels such as111In this paper we use a unit cell with iron-iron bonds along the coordinate axes. Thus, the irreducible representations $B_{\mathrm{1g}}$ and $B_{\mathrm{2g}}$ are interchanged if compared to Ref. [Ikeda et al., 2018], see also Fig. 1. $B_{\mathrm{2g}}$ and $A_{\mathrm{1g}}$, see Fig. 1. In distinction to $\varepsilon_{B_{\mathrm{1g}}}$ that couples directly to $\phi$, strain in other symmetry channels will not wash out the nematic transition, but shift its value. Thus, the nematic transition temperature $T_{\mathrm{nem}}$ remains sharply-defined. Among the key observations of Ref. [Ikeda et al., 2018] are a quadratic suppression of $T_{\mathrm{nem}}$ with antisymmetric strain, i.e.

[TABLE]

with $D_{B_{\mathrm{2g}}}>0$, and a dominant linear variation for symmetric strain

[TABLE]

With regards to the quadratic response of the type $\propto\varepsilon^{2}$ it was further found that $D_{A_{\mathrm{1g}}}\!\ll\!D_{B_{\mathrm{2g}}}$, i.e. the change in the nematic transition temperature due to antisymmetric strain is stronger than due to symmetric strain.

In this paper we analyze the tuning of nematic order through to critical symmetric and antisymmetric strain within the theory of spin-driven nematicity. We demonstrate that the effect of strain is a fluctuation effect, a behavior that is caused by the frustrated nature of the striped magnetic order. We then show that the behavior observed in Ref. [Ikeda et al., 2018] follows naturally within the approach of spin-driven nematicity. The suppression of $T_{\mathrm{nem}}$ by strain is shown to be a consequence of strong classical magnetic fluctuations. As such, strain-enhanced magnetic fluctuations provide a second route to suppress $T_{\mathrm{nem}}$, in addition to increasing quantum fluctuations, see Ref. Maharaj et al. (2017). The finding that $D_{A_{\mathrm{1g}}}\!\ll\!D_{B_{\mathrm{2g}}}$, see Fig. 2, is shown to be a consequence of the anisotropic, i.e. quasi two-dimensional nature of these fluctuations. The analysis is performed for a model of local spins and within a long-wavelength model of collective magnetic fluctuations with composite order. The results of both approaches are fully consistent with each other. The key findings of our analysis of the strain dependence of the nematic ordering temperature are summarized in Fig. 2. Anticipating the similarity in the results for the localized-spin and the long-wavelength approach, we note the generality of the result, i.e., applying equally to systems with localized or itinerant magnetism.

II Strain Tuning nematic order of localized spins

We start our analysis with a lattice model of localized spins on a two-dimensional $J_{1}$-$J_{2}$ Heisenberg model. This model is a specific microscopic realization of the more general field theory for anisotropic three-dimensional systems that we discuss in section III and offers a more microscopic insight into the strain tuning of nematic order. To be specific, we only consider the modification of the nematic transition temperature due to $\varepsilon_{B_{\mathrm{2g}}}$-strain. Without strain, the Hamiltonian is given as

[TABLE]

Here $\boldsymbol{s}(\boldsymbol{r}_{i})$ denotes the spin at the lattice site $\boldsymbol{r}_{i}$, $\langle ij\rangle$ refers to a pair of nearest neighbor spins and $\langle\!\langle ij\rangle\!\rangle$ to next-nearest neighbor spins. The important regime for us is where $J_{2}$ dominates over $J_{1}$. The $J_{1}$-$J_{2}$ Heisenberg model was shown in Ref. [Chandra et al., 1990] to exhibit Ising-nematic order. We closely follow the analysis of Ref. [Chandra et al., 1990], yet add a external strain $\varepsilon_{B_{\mathrm{2g}}}$ that does not couple to the Ising nematic order parameter. Such strain will change the next-nearest neighbor exchange interaction according to

[TABLE]

where the two signs refer to the two diagonal couplings of a square lattice, that corresponds to the single-iron unit cell description of the iron-based materials. The coupling constant $\kappa$ can be determined from a microscopic first-principles calculation. As it always as a prefactor to $\varepsilon_{B_{\mathrm{2g}}}$, we drop $\kappa$ in what follows. Fourier transformation of the Hamiltonian yields

[TABLE]

with

[TABLE]

The classical ground state energy per lattice site

[TABLE]

of this model is unchanged by strain. The energy gain due to the enhancement of one diagonal coupling $J_{2}^{+}$ is off-set by the loss in exchange coupling along the orthogonal diagonal $J_{2}^{-}$, see Fig. 3. This reveals that for this frustrated lattice, the effect of strain only arises due to fluctuation effects.

In what follows we will evaluate how fluctuations affect the correlation length $\xi$ and lead to a the reduction $\Delta E$ of the ground state’s energy density. This program will require a separate discussion of long-wavelength and short-distance fluctuations. When combined, the energy scale $(\Delta E)\xi^{2}$ reveals the characteristic temperature where nematic order appears.

II.1 Long-wave length fluctuations

The long-wavelength fluctuations are described by considering $\boldsymbol{s}(\boldsymbol{k})$ for $\boldsymbol{k}$ near $\boldsymbol{Q}_{x}^{\text{2d}}\!=\!(\pi,0)$ or $\boldsymbol{Q}_{y}^{\text{2d}}\!=\!(0,\pi)$. The associated real-space modulation then reads

[TABLE]

where $\boldsymbol{n}_{a}(\boldsymbol{r})^{2}=1$ are unit vectors and the coordinate $\boldsymbol{r}=(x,y)$ becomes a continuous variable. The two coupled, interpenetrating Néel sub-lattices magnetizations $\boldsymbol{n}_{1,2}$ can be modeled in a nonlinear sigma model [Chandra et al., 1990],

[TABLE]

Here $g=T/2J_{2}s^{2}$ is the stiffness and $\alpha=J_{1}/J_{2}$. One can eliminate the strain term proportional to $\varepsilon_{B_{\mathrm{2g}}}$ through two successive coordinate transformations. First, we introduce $z_{\mp}=(x\mp y)/[2(1\mp\varepsilon_{B_{\mathrm{2g}}})]^{1/2}$, i.e. a rotation by $\pi/4$ with a simultaneous stretching/compression along the two directions. The second transformation $(\tilde{x},\tilde{y})=(z_{+}+z_{-},z_{+}-z_{-})/\sqrt{2}$ rotates the coordinates back (by $-\pi/4$) (without compressing or stretching). The action then takes the form

[TABLE]

with $\tilde{g}=g/(1-\varepsilon_{B_{\mathrm{2g}}}^{2})^{1/2}$ and $\tilde{\alpha}=\alpha/(1-\varepsilon_{B_{\mathrm{2g}}}^{2})^{1/2}$. This is precisely the model without strain, yet with a reduced (effective) next nearest neighbor interaction

[TABLE]

It is important to note that fluctuations are affected by strain while the classical ground-state energy is not. This gives rise to a shortening of the magnetic correlation length $\xi$. For the two-dimensional spin model follows from the usual renormalization group procedure [Chandra et al., 1990] that

[TABLE]

with $a_{0}$ of the order of the lattice constant, $\tilde{z}\!=\!z(\tilde{\alpha})$ and $z=z(\alpha)\equiv 2\alpha/(\arcsin(\alpha)+\alpha\sqrt{1-\alpha})$. The strain-induced reduction of the spin-wave stiffness is the is the predominant effect of strain on long wavelength magnetic fluctuations.

In performing the above coordinate transformations one has to be careful as they may change the boundary conditions and thus the symmetry of the system. This turns out to be a problem if one considers $\varepsilon_{B_{\mathrm{1g}}}$ strain that couples to the Ising-nematic order parameter. For $\varepsilon_{B_{\mathrm{2g}}}$ discussed here this problem does not exist.

II.2 Short-distance fluctuations

Whereas $B_{\mathrm{2g}}$ strain modifies the stiffness of long-wave fluctuations, is will affect the system’s energy through short distance fluctuations. This energy depends on the angle $\theta\equiv\arccos(\boldsymbol{n}_{1\!}\cdot\boldsymbol{n}_{2})$ between the two sub-lattices, see Fig. 3. To quantify the effect of strain we perform a $1/s$ spin-wave analysis of the Heisenberg model (similar to Ref. [Chandra et al., 1990]) with distinct strain-induced exchange interactions along the two diagonals. The spin-wave spectrum of this problems is given as

[TABLE]

where

[TABLE]

and

[TABLE]

This allows us to analyze the free energy corrections

[TABLE]

due to spin-wave excitations, which—in the limit $T\!\to\!0$—corresponds to the correction to the ground state energy $\Delta E=\Delta F(T=0)$. Performing the momentum integration, we obtain the additional biquadratic exchange energy

[TABLE]

where

[TABLE]

The coefficient $\gamma_{Q}^{(0)}\approx 0.26025$ describes the situation without strain and was already given in Ref. [Chandra et al., 1990]. For our considerations it is more important to determine the change in the bi-quadratic interaction due to strain which is characterized by the coefficient $\zeta\approx 0.17273>0$. Thus, the biquadratic spin interaction increases due to strain.

II.3 The nematic transition temperature

Combining long and short wave-length excitations finally allows us to determine the nematic transition temperature $T_{\mathrm{nem}}.$ The order parameter of the nematic state is the Ising variable

[TABLE]

which relates to the angle $\theta$ between the two sublattices via $\cos(\theta)=\phi$. To determine $T_{\mathrm{nem}}$ we consider the typical interaction energy of a region of size $\xi^{2}$. If this energy is comparable to the temperature, one expects the nematic phase transition to take place [Chandra et al., 1990]. This gives rise to the criterion

[TABLE]

We obtain two opposite trends due to external strain: On the one hand, the correlation length gets smaller which reduces the transition temperature. On the other hand, $\Delta E$ gets larger, which tends to enhance $T_{\mathrm{nem}}$. For any system that is near a second-order magnetic phase transition the correlation length is about to diverge. Even for weak first-order transitions, the correlation length above the magnetic ordering temperature is exponentially large, see Eq. (12). As a result, the strain modification of the correlation length is always the dominant one. Thus, we find that within the spin-induced nematic theory, $B_{\mathrm{2g}}$ strain clearly decreases the nematic transition temperature, in agreement with the experimental observation [Ikeda et al., 2018]. This suppression of $T_{\mathrm{nem}}$ is quadratic in $\varepsilon_{B_{\mathrm{2g}}}$ as suggested by the expressions (12) and (19).

Since the analysis of this section was performed for a two-dimensional system with exponentially growing magnetic correlation length, it is important to analyze the role of three-dimensional, albeit anisotropic spin correlations, relevant for many iron-based materials. This analysis will be performed in the next section.

III Strain tuning for spin-induced vestigial order

In this section, we consider the long-wavelength theory of collective magnetic degrees of freedom along the lines of Ref. [Fernandes et al., 2012]. As mentioned above, this approach may equally apply to the collective response of localized spins, or to magnetism of itinerant moments. We consider a generic spin configuration

[TABLE]

where the $\boldsymbol{m}_{x,y}(\boldsymbol{R})$ vary smoothly in space. In distinction to Eq. (8) the fields $\boldsymbol{m}_{a}$ are not unit vectors. As we are analyzing anisotropic three-dimensional systems, the ordering vectors are now given as $\boldsymbol{Q}_{x\!}\!=\!(\pi,0,0)$ and $\boldsymbol{Q}_{y\!}\!=\!(0,\pi,0)$. We perform a continuum’s description for the coordinates within the planes, but keep the discrete lattice structure for the third dimension, with layer index $l$. Thus, we express the three-dimensional coordinates $\boldsymbol{R}=(\boldsymbol{r},la_{z})$ in terms of the two-dimensional vector $\boldsymbol{r}=(x,y)$ and the discrete layer index $l$. Following Ref. [Fernandes et al., 2018] we combine the two vectors into $\boldsymbol{m}=(\boldsymbol{m}_{x},\boldsymbol{m}_{y}).$ and obtain the effective action of the problem

[TABLE]

where the Lagrangian $\mathcal{L}\!=\!\mathcal{L}_{\parallel}+\mathcal{L}_{\perp}$ consists of an intra-layer term

[TABLE]

and a coupling between nearest neighboring layers

[TABLE]

The Pauli matrices $\tau_{\alpha}$ act in the space of two ordering vectors, e.g. $\boldsymbol{m}\tau_{3}\boldsymbol{m}\!=\!\boldsymbol{m}_{x}\!\cdot\!\boldsymbol{m}_{x}\!-\!\boldsymbol{m}_{y}\!\cdot\!\boldsymbol{m}_{y}$. Depending on the sign and the magnitude of the coupling constants $u$, $g$ and $v$, several magnetic phases and their associated vestigial orders have been discussed [Fernandes et al., 2016]: These are the stripe antiferromagnetic ($\boldsymbol{m}\tau_{3}\boldsymbol{m}$), the charge-spin density-wave ($\boldsymbol{m}\tau_{1}\boldsymbol{m}$). A third nematic order $\!\propto\!\boldsymbol{m}_{1}\!\times\!\boldsymbol{m}_{2}$ (not considered here), associated with a spin-vortex crystal has been identified. Microscopic expressions for the phenomenological parameters are given elsewhere, see Ref. [Fernandes et al., 2012]. Note that our analysis assumes non-critical responses of these phenomenological parameters to strain. More specifically their strain-dependence is neglected in the following. Tuning of nematicity is then caused by the modification of critical fluctuations at finite strain.

In the presence of strain $\varepsilon_{\alpha\beta}$ the response of collective spin modes is governed by a modified action

[TABLE]

with the stress tensor

[TABLE]

In our subsequent analysis we will work at constant strain. This is clearly adequate to describe measurements such as the elasto-resistivity [Chu et al., 2012]. On the other hand, experiments on unstrained samples should rather be performed at fixed stress. For the description of Ising nematic order this requires to include harmonic elastic terms, characterized by bare elastic constants. As shown in Ref. [Fernandes et al., 2010] this gives rise to an enhancement of the nematic coupling constant, hence enhancing the nematic transition temperature compared to the value of a purely electronic system. For the experiments of Ref. [Ikeda et al., 2018] one has to keep in mind that fixed strain $\varepsilon_{B_{\mathrm{2g}}}\neq 0$ still corresponds to fixed stress in the other symmetry channels. Thus, the above renormalizations of the nematic coupling due to fluctuations of $\varepsilon_{B_{\mathrm{1g}}}$ should nevertheless be included. We will expand on these issues below, when we make contact to experiment.

In the following, we focus on in-plane strain; in particular we introduce the strain combinations

[TABLE]

and similar for the stress tensor. In this notation the strain-stress term in Eq. (26) takes the form

[TABLE]

For the trivially transforming $A_{\mathrm{1g}}$ combination follows

[TABLE]

while the stress associated to the two nontrivial irreducible representations $B_{\mathrm{1g}}$ and $B_{\mathrm{2g}}$ read

[TABLE]

The combination $(T_{xy}-T_{yx})/2=\alpha\partial_{y}\boldsymbol{m}\tau_{3}\partial_{x}\boldsymbol{m}$ describes the response to a rotation about the $z$ axis and will not be discussed further.

In our treatment of anisotropic three-dimensional systems, the dispersion in the third direction is characterized by the coupling between neighboring layers of Eq. (25). To simplify our notation we use $\boldsymbol{q}^{2}=q_{x}^{2}+q_{y}^{2}+q_{z}^{2}$ where $q_{z}^{2}$ stands in fact for $2q_{0}^{2}[1-\cos(q_{z}/q_{0})]$. Here $q_{0}$ is a measure of the in-plane versus out-of-plane anisotropy. The 2d case is obtained in the limit $q_{0}\!\to\!0$, while for $q_{0}\!\to\!\infty$ the isotropic 3d case is recovered. By parametrizing the in-plane momentum vector as $(q_{x},q_{y})=q(\cos\varphi,\sin\varphi)]$ (consequently $\boldsymbol{q}^{2}=q^{2}+q_{z}^{2}$), and after introducing an ulta-violett cut-off $q<\Lambda$, the momentum integral of a function $f(\boldsymbol{q})$ takes the form

[TABLE]

III.1 Analysis without strain

The subsequent analysis extends the approach used in Ref. [Fernandes et al., 2012] to the case of finite strain. For completeness, we briefly summarize the zero-strain case. After introducing the Hubbard-Stratonovich fields $\eta$ and $\phi$ that are conjugate to $\boldsymbol{m}\tau_{0}\boldsymbol{m}$ and $\boldsymbol{m}\tau_{3}\boldsymbol{m}$ respectively, we can integrate out the quadratic action for the fields $\boldsymbol{m}$ and obtain an effective action for the conjugate fields

[TABLE]

where the Green’s functions associated with $\boldsymbol{m}_{1,2}$ read

[TABLE]

In the large-$N$ limit, the partition function may be evaluated at the saddle-point solution. Then

[TABLE]

emerges as the Ising nematic order parameter.

Expressing the saddle-point equations in Fourier space, we assume here that $\eta$ and $\phi$ are not coordinate-dependent, we obtain the self-consistency equations

[TABLE]

where we introduced $r=r_{0}-i3u\eta$ and used the convention of Eq. 35.

In order to highlight solution of this set of equations, let us consider the case $\alpha=0$ first. The more general case $\alpha\neq 0$ will be treated below when we include finite strain. The renormalization of $r_{0}\!\to\!r$ gives rise to a change in the in-plane magnetic correlation length $\xi$, via

[TABLE]

In absence of nematic order $\phi\!=\!0$ we use $\bar{r}\!\equiv\!r(\phi=0)$ which obeys

[TABLE]

where $I_{1}(\bar{r})$ belongs to a series of reappearing integrals $I_{n}(\bar{r})$, for which we introduce a unified notation, see Appendix A. Treating the onset of nematicity perturbatively we write $r\!\approx\!\bar{r}+d\phi^{2}$ and find $d\!=\!(3g)^{2}6uI_{3}(\bar{r})/[1+6uI_{2}(\bar{r})]$ from Eq. (39). Similarly, Eq. (40) takes the form

[TABLE]

This is the equation of state for $\phi$ that can be interpreted as being due to an effective $\phi^{4}$-theory,

[TABLE]

for the nematic order parameter: At high temperature, when the coefficient $1-6gI_{2}(\bar{r})$ is positive, no nematic order exists. A vanishing of this coefficient defines the nematic transition temperature $T_{\mathrm{nem}}$ via the temperature dependence of $\bar{r}$. Substituting $I_{2}(\bar{r})$ by its expression in an anisotropic 3d system, as derived in Appendix A, the nematic phase sets in when $\sqrt{\bar{r}(\bar{r}+4q_{0}^{2})}=3gq_{0}/2\pi$. The sign of the square bracket in Eq. (43) at the onset of nematicity [proportional to $3g(u-2g)-16\pi^{2}q_{0}^{2}$ in the anisotropic 3d system] discriminates between a first (negative) and second order (positive) phase transition. One readily sees that an isotropic three dimensional system always undergoes a first order nematic transition, while split second order transitions occur in sufficiently anisotropic systems, see also Ref. [Fernandes et al., 2012].

III.2 Strain tuning Ising nematic order

III.2.1 $A_{\mathrm{1g}}$-strain

When applying an external strain in the $A_{\mathrm{1g}}$-channel, $\varepsilon_{xx}\!=\!\varepsilon_{yy}\!=\!\varepsilon_{A_{\mathrm{1g}}}$ the inverse Green’s functions (37) are modified to

[TABLE]

with $\gamma=(1+2\varepsilon_{A_{\mathrm{1g}}})$ renormalizing the bare parameters $r_{0}$, $u$, and $g$ as well as the gradient term. From this renormalization follows immediately, that $\varepsilon_{A_{\mathrm{1g}}}$ affects the transition temperatures to linear order. The self-consistency equations for $r$ and $\phi$ now read

[TABLE]

with the integrand’s denominator

[TABLE]

Considering the effects of $\alpha$ and $\varepsilon_{A_{\mathrm{1g}}}$ as perturbations to the system, motivates the Ansatz

[TABLE]

Specifically, the above self-consistency equations respond linearly to an external strain while the combination $q_{x}^{2}-q_{y}^{2}$ only allows for even powers in $\alpha$. By expanding the first equation and equating coefficients we find $d$ as in the case without strain as well as

[TABLE]

with the functions $J_{n}^{{\ell,m}}(\bar{r})$ defined in Appendix A. Notice, if the index $n$ subscript exceeds the sum of the superscripts $\ell+m$ by more than 1, the integral is convergent for the anisotropic 3d case. For $n\!=\!\ell\!+\!m\!+\!1$ the integral is logarithmically divergent. Expanding the second self-consistency-equation (47) in the small parameters $\psi$, $\alpha$ and $\varepsilon_{A_{\mathrm{1g}}}$ yields the equation of state $0\!=\!a_{A_{\mathrm{1g}}}\!(\varepsilon_{A_{\mathrm{1g}}})\phi+b\phi^{3}$, see Eq. (43), with the linear coefficient

[TABLE]

The coefficient $a_{A_{\mathrm{1g}}}\!(\varepsilon_{A_{\mathrm{1g}}})$ now determines the ground state at finite strain and, in particular, the onset of the nematic order. As anticipated the $\mathrm{A}_{\mathrm{1g}}$-strain changes linearly the onset of nematic order, as featured by a term $\propto\!\varepsilon_{A_{\mathrm{1g}}}$ in the linear coefficient $a$ of the effective $\phi^{4}$ theory.

Anticipating a discussion with the results obtained in other strain channels, the expansion in $\varepsilon_{A_{\mathrm{1g}}}$ has been carried to quadratic order. In this context it is important to note that the coefficient $f_{A_{\mathrm{1g}}}$ and the square bracket determining the quadratic correction $\propto\!\varepsilon_{A_{\mathrm{1g}}}^{2}\!$ remain bounded for strongly anisotropic systems. As we shall see, this is in contrast to the quadratic terms in the $\mathrm{B}_{\mathrm{1g}}$ and $\mathrm{B}_{\mathrm{2g}}$ channels where the coefficient $\propto\!\varepsilon_{B_{\mathrm{2g}}}^{2}\!$ logarithmically diverges upon approaching the two-dimensional limit.

III.2.2 $B_{\mathrm{1g}}$-strain

Next, we consider an external strain with the components $\varepsilon_{xx}\!=\!-\varepsilon_{yy}\!=\!\varepsilon_{B_{\mathrm{1g}}}$. Following the derivation scheme discussed above, the layer-diagonal part of inverse Green’s functions (37) will be augmented by additional terms while the inter-layer terms remain unchanged

[TABLE]

The self-consistency equations for $r$ and $\psi$ now read

[TABLE]

with

[TABLE]

In the equation for $\phi$, the term $\propto\!\alpha\varepsilon_{B_{\mathrm{1g}}}$ of the integrand’s numerator implies that the nematic order exists at any temperature. In fact, this term generically produces a finite value for $\phi$. As shall become clear below, the symmetry-breaking term acts as an Ising field $h_{\mathrm{eff}}$ lifting the phase transition. An expansion of Eq. (53) using the Ansatz $r\approx\bar{r}+d\phi^{2}+d^{\prime}\alpha^{2}+f_{B_{\mathrm{1g}}}\varepsilon_{B_{\mathrm{1g}}}^{2}$ yields the known values for $d$ and $d^{\prime}$ as well as $f_{B_{\mathrm{1g}}}=4d^{\prime}$. Expanding the equation of state (54) for $\phi$ provides

[TABLE]

with effective source field $h_{\mathrm{eff}}\!=\!4\alpha\varepsilon_{B_{\mathrm{1g}}}[J_{2}^{0,1}(\bar{r})-2J_{3}^{2,0}(\bar{r})]$ induced by strain. Since $\phi$ and $\varepsilon_{B_{\mathrm{1g}}}$ transform according to the same irreducible representation, it is natural to expect the nemato-elastic term $\lambda_{\mathrm{n.e.}}\phi\varepsilon_{B_{\mathrm{1g}}}$ in the free energy expansion. Our analysis shows that the nemato-elastic coupling constant is

[TABLE]

and simplifies to $\lambda_{\mathrm{n.e.}}=\alpha q_{0}/2\pi$ in the 2d limit. Subleading corrections in $\varepsilon_{B_{\mathrm{1g}}}$ or $\alpha$ to the Landau parameters $a$ and $b$ have been omitted here222Quadratic corrections in $\alpha$ and $\varepsilon_{B_{\mathrm{1g}}}$ on the left-hand side of Eq. (56) turn out to be identical to the corrections in Eq. (66). Only in the accidental case where either $\varepsilon_{B_{\mathrm{1g}}}$ or $\alpha$ identically vanishes, the nematic transition survives and is merely shifted by a term proportional to $\alpha^{2}$ or $\varepsilon_{B_{\mathrm{1g}}}^{2}$ respectively.. The parameters $\alpha$ and $\varepsilon_{B_{\mathrm{1g}}}$ now appear to linear order on the left-hand side of the above equation and generically induce a finite nematic order at any temperature. The nematic order is

[TABLE]

The denominator indicates the softening of $\phi$ when approaching the underlying (for $\alpha\varepsilon_{B_{\mathrm{1g}}}\!\!=\!0$) nematic transition. This trend is also reflected in the susceptibility

[TABLE]

which diverges at the transition temperature.

Let us note here that our analysis is performed at fixed strain. The nematic transition hence corresponds to that observed in measurements of the elastoresistivity, with resistivity anisotropy $\Delta\rho$, as well as that obtained from the Raman response $R_{B_{\mathrm{1g}}}(\omega)$ in the $B_{\mathrm{1g}}$ channel [Gallais et al., 2013; Gallais and Paul, 2016; Karahasanovic et al., 2015]

[TABLE]

This susceptibility $\tilde{\chi}_{\mathrm{nem}}$ differs from the ’true’ thermodynamic nematic quantity

[TABLE]

where $h_{\mathrm{nem}}$ is a field conjugate to the nematic order parameter and $g_{\mathrm{ren}\!}\!=\!g\!+\!\lambda_{\mathrm{n.e.}}/C^{(0)}_{66}$ is the nematic coupling constant renormalized by the coupling to elastic degrees of freedom [Fernandes et al., 2010]. Here $C^{(0)}_{66}$ is the bare (high-temperature) value of the elastic constant. The renormalized elastic constant is then $C_{66}\!=\!C_{66}^{(0)}/(1+\lambda_{\mathrm{n.e.}}^{2}\chi_{\mathrm{nem}})$. The true susceptibility $\chi_{\mathrm{nem}}$ diverges at the nematic transition temperature for constant stress. The fact that $g\!<\!g_{\mathrm{ren}}$ is the reason why the elastoresistivity or $\tilde{\chi}_{\mathrm{nem}}$ obtained from Raman measurements remain finite, and displays a temperature dependence $\tilde{\chi}_{\mathrm{nem}}\!\propto\!(T-\Theta_{\mathrm{C}})^{-1}$ with the Curie temperature $\Theta_{\mathrm{C}}$ below $T_{\mathrm{nem}}$ [Gallais and Paul, 2016; Karahasanovic et al., 2015]. The crucial distinction between these two susceptibilities was exploited in Ref. [Chu et al., 2012] to conclude that the origin of the nematic transition in the iron-based systems is electronic.

III.2.3 $B_{\mathrm{2g}}$-strain

The nematic symmetry is preserved for strain with the components $\varepsilon_{xy}\!=\!\varepsilon_{yx}\!=\!\varepsilon_{B_{\mathrm{2g}}}$. Here, the intra-plane inverse Green’s functions (37) are simply augmented by the term $4\varepsilon_{B_{\mathrm{2g}}}\delta(\boldsymbol{r}\!-\!\boldsymbol{r}^{\prime})\partial_{x^{\prime}}\partial_{y^{\prime}}$. Now the self-consistency equations for $r$ and $\psi$ take the form

[TABLE]

where

[TABLE]

The systematic appearance of the strain $\varepsilon_{B_{\mathrm{2g}}}$ in combination with the odd function $q_{x}q_{y}$ in momentum space implies that all corrections to the unstrained case are (at least) quadratic in $\varepsilon_{B_{\mathrm{2g}}}$. An expansion similar to the earlier ones, with $r\approx\bar{r}+d\phi^{2}+d^{\prime}\alpha^{2}+f_{B_{\mathrm{2g}}}\varepsilon_{B_{\mathrm{2g}}}^{2}$ provides

[TABLE]

from Eq. (62). The positive coefficient $f_{B_{\mathrm{2g}}}$ implies that the correlation length $\xi\!\sim\!(r-3g\phi)^{-1/2}$ is shortened by strain, pointing towards a generic suppression of the nematic transition temperature. The equation of state (62) takes the form $0=a_{B_{\mathrm{2g}}}\!(\varepsilon_{B_{\mathrm{2g}}})\phi+b\phi^{3}$, see Eq. (43), with

[TABLE]

Note that in the limit of quasi two-dimensional systems, the coefficient determining the shift of nematic transition is dominated by the logarithmic term in $J_{3}^{2,0}(\bar{r})$ (constituent of $d^{\prime}$). Its magnitude is also decisive in determining the sign of the shift giving rise to a suppression of the nematic transition temperature.

III.2.4 Discussion

In the following discussion we focus on the (strongly) anisotropic limit $q_{0}\!\to\!0$ relevant for many layered magnetic, or electronic systems. The Landau coefficients (51) and (66) can be expanded in the vicinity of the unstrained value $a$, for simplicity we assume $\alpha=0$ here. This expansion requires using the expressions (76) and (78) listed in the Appendix A as well as the solution $\bar{r}\approx 3gq_{0}/2\pi$ valid in the limit $q_{0}\ll g$. By associating the Landau coefficient with a nematic transition temperature via $a=(1-T_{\mathrm{nem}}/T)$, we identify the phenomenological parameters

[TABLE]

as defined in Eqs. (1) and (2), and thereby quantify the effect of strain on the nematic phase transition.

For $A_{\mathrm{1g}}$ strain the quadratic correction is regular, while that of the $B_{\mathrm{2g}}$ channel is logarithmically divergent in the 2d limit. In light of the experimental finding [Ikeda et al., 2018], that the response of the nematic phase boundary in $\mathrm{Co}$-doped $\mathrm{Ba}\mathrm{Fe}_{2}\mathrm{As}_{2}$ to the symmetric strain333The $B_{\mathrm{2g}}$ strain effects are to be mapped to the $B_{\mathrm{1g}}$ sector in the experimental work [Ikeda et al., 2018], where the nematic axes point along the crystallographic diagonals (1,1) and (-1,1). $B_{\mathrm{2g}}$ as compared to that of the $A_{\mathrm{1g}}$ sector is very large, we compute the ratio $D_{B_{\mathrm{2g}}}/D_{A_{\mathrm{1g}}}$. Figure 4 shows this ratio for fixed $q_{0}$ as a function of the material parameters $u$ and $g$, while the anisotropy-dependence for fixed $u$ and $g$ is shown as an inset. The figure covers a parameter range $-1\lesssim\log(u/g)\lesssim 1$ and $q_{0}\ll 1$ in agreement with qualitative estimates for $\mathrm{Co}$-doped $\mathrm{Ba}\mathrm{Fe}_{2}\mathrm{As}_{2}$. Though a derivation of the phenomenological parameters from a full microscopic treatment exists, e.g. Ref. [Fernandes et al., 2012], the evaluation of an accurate numerical value remains difficult. As $g$ emerges from a perturbation theory Fernandes et al. (2012), the ratio $u/g\approx\{4\mu m_{x}m_{y}/[\epsilon_{0}m(m_{x}-m_{y})]\}^{2}$ [with the chemical potential $\mu$, the offset energy $\epsilon_{0}$, and the hole ($m$) and electronic ($m_{x}$, $m_{y}$) band masses] is expected to be large, yet of order unity. Disorder effects Hoyer et al. (2016) manifestly increase $g$, hence decrease the ratio $u/g$. We estimate the anisotropy parameter $q_{0}$ for this pnictide compound from $q_{0}^{2}\sim J_{ab}/J_{c}\approx 30$ as reported in Ref. [Dai, 2015].

The representations in Fig. 4 support that the $B_{\mathrm{2g}}$ strongly affects the nematic order; possibly to the point of reaching beyond linear effects $\propto C_{A_{\mathrm{1g}}}\varepsilon_{A_{\mathrm{1g}}}$ from the $A_{\mathrm{1g}}$ sector. Although derived in the strongly anisotropic limit, $q_{0}\!\to\!0$, Eq. (67) points towards a sign change of $D_{B_{\mathrm{2g}}}$ when moving away from the two-dimensional limit; a possibility that depends on the specific parameters in the problem.

IV Conclusion

In conclusion, we have investigated the response of a nematic order to strain belonging to different symmetry classes ($\varepsilon_{A_{\mathrm{1g}}}$, $\varepsilon_{B_{\mathrm{1g}}}$, $\varepsilon_{B_{\mathrm{2g}}}$). Hereby special attention was given to the symmetry channel $B_{\mathrm{2g}}$, for which recent experimental work [Ikeda et al., 2018] has found a surprisingly strong suppression of the nematic transition temperature $T_{\mathrm{nem}}$. Our analysis of long- and short-length fluctuations in a $J_{1}$-$J_{2}^{\pm}$ spin model—the latter implements a $B_{\mathrm{2g}}$-strain on a lattice—provides us with clear indications that the nematic transition temperature decreases with strain. This is due to the strain-induced shortening of the magnetic correlation length $\xi$.

These findings, combined with more general symmetry considerations lead to several observations. The degeneracy-lifting $B_{\mathrm{1g}}$ strain is expected to act as an effective source field $h_{\mathrm{eff}}$, thus replacing the nematic transition by a smooth cross-over. In contrast, the symmetry-conserving strains of the $A_{\mathrm{1g}}$ and $B_{\textrm{2g}}$ type preserve the transition and merely yield a shift $\Delta T_{\mathrm{nem}}=T_{\mathrm{nem}}(\varepsilon)-T_{\mathrm{nem}}^{0}$ in the transition temperature. In the $A_{\mathrm{1g}}$ channel, we find $\Delta T_{\mathrm{nem}}=C_{A_{\mathrm{1g}}}\varepsilon_{A_{\mathrm{1g}}}-D_{A_{\mathrm{1g}}}\varepsilon_{A_{\mathrm{1g}}}^{2}$, with a dominant linear contribution. In the $B_{\mathrm{2g}}$ sector, the deviation is quadratic, i.e.  $\Delta T_{\mathrm{nem}}=-D_{B_{\mathrm{2g}}}\varepsilon_{B_{\mathrm{2g}}}^{2}$.

Approaching the problem from a field-theoretical approach we provide a quantitative tool to evaluate the parameters $h_{\mathrm{eff}}$, $C_{A_{\mathrm{1g}}}$, $D_{A_{\mathrm{1g}}}$ and $D_{B_{\mathrm{2g}}}$ from the underlying Lagrangian formalism. In this approach, the strong response of nematic order to symmetric $B_{\mathrm{2g}}$, i.e. $D_{B_{\mathrm{2g}}}\gg D_{A_{\mathrm{1g}}}$, is found to be a generic feature of strongly anisotropic, quasi-2d systems, see Fig. 2; as represented by $\mathrm{Co}$-doped $\mathrm{Ba}\mathrm{Fe}_{2}\mathrm{As}_{2}$.

A related issue is that other emergent vestigial phases are known to appear in the spin-based scenario of iron-based materials [Fernandes et al., 2016, 2018]. Furthermore, charge density wave order is expected to feature similar responses in a situation of $C_{4}$-symmetric collinear double-$\boldsymbol{Q}$ spin order; as observed in a number of systems [Kim et al., 2010; Hassinger et al., 2012; Avci S. et al., 2014; Wang et al., 2016; Böhmer A. E. et al., 2015; Allred et al., 2015; Hassinger et al., 2016; Allred J. M. et al., 2016]. Finally, chiral order is related to $C_{4}$-symmetric non-collinear double-$\boldsymbol{Q}$ spin order that forms spin vortex crystals observed in Ref. [Meier et al., 2018]. This embeds our work in a broader quest for understanding how strain affects vestigial phases of matter.

Acknowledgements.

We are grateful to I. R. Fisher, M. Ikeda, J. C. Palmstrom, and P. Walmsley, for stimulating discussions. J. S.  was funded by the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grant GBMF4302 while visiting the Geballe Laboratory for Advanced Materials at Stanford University.

Appendix A Momentum Integrals

The derivation of the results in the main text involve to compute several integrals of the form $\int_{{\boldsymbol{q}}}f(\bar{r},{\boldsymbol{q}})$, which we define in the following as

[TABLE]

Angular integration in the $xy$ plane provides the relations

[TABLE]

with $z!\equiv\Gamma(z+1)$ and $(2z)!!=2^{z}z!$ defined through the Gamma function $\Gamma(\zeta)\equiv\int_{0}^{\infty}dt\,t^{\zeta-1}e^{-t}$ and $\mathrm{B}(z)=[1+(-1)^{z}]/2$ the Boolean parity function. Furthermore, we have $J_{n}^{0,0}(r)=K_{n}^{0,0}(r)=I_{n}(r)$. The remaining task is to deterine the integrals $I_{n}(r)$.

Integrals in the Anisotropic 3d limit

Next we evaluate the integrals $I_{n}(r)$ in the anisotropic 3d case. We substitute $q_{x}\!\to\!q\cos(\varphi)$, $q_{y}\!\to\!q\sin(\varphi)$, and $q_{z}^{2}\!\to\!2q_{0}^{2}[1-\cos(q_{z}/q_{0})]$ with $q_{0}$ a measure of the uniaxial anisotropy; now the momentum-space integration along $z$ is restricted to $|q_{z}|<\pi/q_{0}$. With an ultraviolet cut-off $\Lambda$ for the in-plane momentum $q=(q_{x}^{2}+q_{y}^{2})^{1/2}$, $q<\Lambda$, the momentum-integral are mapped to

[TABLE]

Within this mapping, the integrals defined in Eq. (70) evaluate to

[TABLE]

with ${}_{2}F_{1}$ the hypergeometric function. The limit $\Lambda\!\to\!\infty$ is taken for all convergent integrals, while only the dominant terms in $\Lambda$ are considered otherwise.

Below, we evaluate some integrals of the form (71) and express them in terms of the series of $I_{n}$ functions.

[TABLE]

More generically we find

[TABLE]

with $C(n-1,m)=(n-1)!/[(n-m-1)!m!]$ the binomial coefficients, $c_{n}=\sum_{k=1}^{n-1}(1/k)=\Gamma^{\prime}(n)/\Gamma(n)+\gamma_{{\scriptscriptstyle\mathrm{E}}}$ positive constants related to the Polygamma function, and $\gamma_{{\scriptscriptstyle\mathrm{E}}}\approx 0.577$ Euler’s constant.