Disorder-Induced Electronic Nematicity

TL;DR

This paper develops a theoretical framework showing how disorder can locally induce nematic order in materials with strong nematic fluctuations, explaining recent experimental observations in iron-based superconductors.

Contribution

It introduces a symmetry-based Ginzburg-Landau model and microscopic calculations to explain disorder-induced nematicity and its enhancement of the transition temperature.

Findings

Single impurities induce local nematicity above the transition temperature.

Disorder causes inhomogeneous nematicity with zero spatial average.

Finite disorder can enhance the nematic transition temperature.

Abstract

We expose the theoretical mechanisms underlying disorder-induced nematicity in systems exhibiting strong fluctuations or ordering in the nematic channel. Our analysis consists of a symmetry-based Ginzburg-Landau approach and associated microscopic calculations. We show that a single featureless point-like impurity induces nematicity locally, already above the critical nematic transition temperature. The persistence of fourfold rotational symmetry constrains the resulting disorder-induced nematicity to be inhomogeneous and spatially average to zero. Going beyond the single impurity case, we discuss the effects of finite disorder concentrations on the appearance of nematicity. We identify the conditions that allow disorder to enhance the nematic transition temperature, and we provide a concrete example. The presented theoretical results can explain a large series of recent experimental discoveries of disorder-induced nematic order in iron-based superconductors.

Equations

$$\mathcal{F}(\mathbf{r})$$

$$\displaystyle\big{[}\alpha(T-T_{\rm nem})-c\bm{\nabla}^{2}\big{]}N(\bm{r})=-g\left(\partial_{x}^{2}-\partial_{y}^{2}\right)V(\bm{r})\,.\quad$$

$$\displaystyle N(\bm{r})=\int\frac{d\bm{q}}{(2\pi)^{2}}e^{i\bm{q}\cdot\bm{r}}\phantom{.}\frac{g}{c}\frac{\big{(}q_{x}^{2}-q_{y}^{2}\big{)}V(\bm{q})}{\bm{q}^{2}+\xi_{\rm nem}^{-2}}\,,$$

$$\int d\mathbf{r}N(\mathbf{r})=N(\mathbf{q}=0)\propto {\int }_0^{2\pi }d\varphi N(r,\varphi )=0.$$

$$\displaystyle\delta{\cal F}(\bm{r})=-\bigg{\{}g^{\prime}V(\bm{r})+g^{\prime\prime}\big{[}V(\bm{r})\big{]}^{2}\bigg{\}}\big{[}N(\bm{r})\big{]}^{2}/2\,.$$

$$\alpha (T-T_{\mathrm{nem}})N(\mathbf{q}=0)=g^{\prime }\int d\mathbf{p}V(\mathbf{p})N(-\mathbf{p})$$

$$+g^{\prime \prime }\iint d\mathbf{p}d{\mathbf{p}}^{\prime }V({\mathbf{p}}^{\prime })V(\mathbf{p}-{\mathbf{p}}^{\prime })N(-\mathbf{p}).$$

$$\left[\alpha (T-T_{\mathrm{nem}})-g^{\prime \prime }\int d\mathbf{p}|V(\mathbf{p}){|}^2\right]N(\mathbf{q}=0)$$

$$\displaystyle=\frac{gg^{\prime}}{c}\int d\bm{p}\phantom{.}\frac{\big{(}p_{x}^{2}-p_{y}^{2}\big{)}|V(\bm{p})|^{2}}{\bm{p}^{2}+\xi_{\rm nem}^{-2}}\,.\quad$$

$$T_{\mathrm{nem}}^{\mathrm{imp}}=T_{\mathrm{nem}}+\frac{g^{\prime \prime }}{\alpha }\int d\mathbf{p}|V(\mathbf{p}){|}^2.$$

$$\displaystyle N_{\bm{R}}=-V_{\rm nem}\sum_{\bm{\delta}}f_{\bm{\delta}}\big{<}c^{{\dagger}}_{\bm{R}+\bm{\delta}}c_{\bm{R}}+c^{{\dagger}}_{\bm{R}}c_{\bm{R}+\bm{\delta}}\big{>}\,,$$

$$\hat{\mathcal{H}}$$

$$\frac{1}{{\tau }_{\mathbf{k}}}=2\pi n_{\mathrm{imp}}{V}^2\frac{1}{\mathcal{N}}\sum _{\mathbf{p}}{\delta }_{{\epsilon }_{\mathbf{p}},{\epsilon }_{\mathbf{k}}}.$$

$$N=-V_{\mathrm{nem}}\frac{1}{\mathcal{N}}\sum _{\mathbf{k}}{f}_{\mathbf{k}}{\int }_{-\infty }^{+\infty }\frac{d\epsilon }{2\pi }\frac{n_F(\epsilon +{\epsilon }_{\mathbf{k}}+N{f}_{\mathbf{k}})}{{\epsilon }^2+1/(2{\tau }_{\mathbf{k}}{)}^2}\frac{1}{{\tau }_{\mathbf{k}}},$$

$$\frac{1}{V_{\mathrm{nem}}}=-\frac{1}{\mathcal{N}}\sum _{\mathbf{k}}{f}_{\mathbf{k}}^2{\int }_{-\infty }^{+\infty }\frac{d\epsilon }{2\pi }\frac{n_F^{\prime }(\epsilon +{\epsilon }_{\mathbf{k}})}{{\epsilon }^2+1/(2{\tau }_{\mathbf{k}}{)}^2}\frac{1}{{\tau }_{\mathbf{k}}}.$$

$$\frac{\partial \mathcal{F}}{\partial N}={\partial }_x\frac{\partial \mathcal{F}}{\partial ({\partial }_{x}N)}+{\partial }_y\frac{\partial \mathcal{F}}{\partial ({\partial }_{y}N)}$$

$$\displaystyle\big{[}\alpha(T-T_{\rm nem})-c\bm{\nabla}^{2}\big{]}N(\bm{r})+\beta[N(\bm{r})]^{3}\quad$$

$$=-g\left({\partial }_x^2-{\partial }_y^2\right)V(\mathbf{r}).$$

$$N(r,\varphi )=\cos (2\varphi ){\int }_{+\infty }^0\frac{q\mathrm{d}q}{2\pi }\frac{g}{c}\frac{q^{2}V(q)}{q^2+{\xi }_{\mathrm{nem}}^{-2}}{J}_2(qr),$$

$$N(r,\varphi )=\frac{\gamma }{c}\left[I_2\left(\frac{r}{{\xi }_{\mathrm{nem}}}\right)-L_{-2}\left(\frac{r}{{\xi }_{\mathrm{nem}}}\right)\right]\cos (2\varphi ),$$

$$\displaystyle\int{\rm d}\bm{r}\phantom{.}[N(\bm{r})]^{n}[V(\bm{r})]^{m}\big{(}\partial_{x}^{2}-\partial_{y}^{2}\big{)}^{\ell}V(\bm{r})$$

$$\displaystyle\int{\rm d}\bm{r}\phantom{.}[N(\bm{r})]^{2n}\phantom{{\dagger}}{\rm and}\phantom{{\dagger}}\int{\rm d}\bm{r}\phantom{.}[N(\bm{r})]^{2n+1}\big{(}\partial_{x}^{2}-\partial_{y}^{2}\big{)}V(\bm{r})\,.$$

$$\displaystyle N(\bm{r})\propto[N(\bm{r})]^{2n-1}\phantom{.}{\rm and}\phantom{.}N(\bm{r})\propto[N(\bm{r})]^{2n}\big{(}\partial_{x}^{2}-\partial_{y}^{2}\big{)}V(\bm{r})\,.$$

$$N(\mathbf{q}=0)\propto$$

$$\int \mathrm{d}{\mathbf{p}}_1\dots \mathrm{d}{\mathbf{p}}_{2n-1}N({\mathbf{p}}_1)\dots N({\mathbf{p}}_{2n-1})\delta \left(\sum _s^{2n-1}{\mathbf{p}}_s\right).$$

$$\displaystyle\bar{N}(\bm{q})=\frac{g}{c}\frac{\big{(}q_{x}^{2}-q_{y}^{2}\big{)}V(\bm{q})}{\bm{q}^{2}+\xi_{\rm nem}^{-2}}\equiv\cos(2\theta)\phantom{.}\frac{g}{c}\frac{q^{2}V(q)}{q^{2}+\xi_{\rm nem}^{-2}}\,,$$

$${\int }_0^{2\pi }\mathrm{d}{\theta }_1\dots {\int }_0^{2\pi }\mathrm{d}{\theta }_{2n-2}\cos (2{\theta }_1)\dots \cos (2{\theta }_{2n-2})\cdot$$

$$\displaystyle\left[\sum_{s=1}^{2n-2}p_{s}^{2}\cos\big{(}2\theta_{s}\big{)}+\sum_{s\neq\ell}^{2n-2}p_{s}p_{\ell}\cos\big{(}\theta_{s}+\theta_{\ell}\big{)}\right]=0\,,\qquad$$

$$\displaystyle\quad\big{[}\alpha(T-T_{\rm nem})-c\bm{\nabla}^{2}\big{]}N(\bm{r})+\beta[N(\bm{r})]^{3}\quad$$

$$\displaystyle\phantom{{\dagger}}=-g\big{(}\partial_{x}^{2}-\partial_{y}^{2}\big{)}V(\bm{r})+g^{\prime}N(\bm{r})V(\bm{r})\,.\qquad\qquad$$

$$\displaystyle\big{[}2\alpha(T_{\rm nem}-T)/c-\bm{\nabla}^{2}\big{]}\delta N(\bm{r})\qquad$$

$$=-\frac{g}{c}\left({\partial }_x^2-{\partial }_y^2\right)V(\mathbf{r})+\frac{g^{\prime }}{c}{N}_{\mathrm{B}}V(\mathbf{r}).$$

$$\delta N(r,\varphi )$$

Full text

The study of electronic nematic quantum phases NemReview is becoming increasingly important in condensed matter systems due to a growing class of recently discovered materials exhibiting nematic behavior, i.e. spontaneous generation of spatial anisotropy. Nematicity has been experimentally identified in a number of correlated quantum materials NemReview , including quantum Hall states in higher Landau levels of 2D electron gases QHEa ; QHEb , bilayer strontium ruthenates Ruthenates , cuprate high-temperature superconductors (SCs) Cuprates , doped Bi2Se3 SCs Bi2Se3a ; Bi2Se3b ; yonezawa , Fe-based SCs (FeSCs) tanatar ; chu ; ying ; chu2 ; blomberg ; ishida13 ; kuo14 ; yi11 ; kostin ; zhao ; nakajima ; dusza ; LGreene ; kasahara and, possibly, twisted bilayer graphene Pashupathy18 . Thus, nematicity begins to establish as a universal electronic state of matter, which motivates further theoretical studies of its distinct properties.

Nematic phases are particularly prevalent in FeSCs, where experimental evidence for electronic nematicity comes from a wide range of techniques, including transport measurements tanatar ; chu ; ying ; chu2 ; blomberg ; ishida13 ; kuo14 , angle-resolved photoemission spectroscopy yi11 , scanning tunneling spectroscopy kostin , neutron scattering zhao , light spectroscopy nakajima ; dusza ; gallais2013 , Andreev-point-contact measurements LGreene and torque magnetometry kasahara . In this case, the emergence of nematicity refers to the spontaneous breaking of fourfold (C4) rotational symmetry. Notably, the identification of the driving mechanism of nematicity in these systems is complicated, due to the coupling of spin, orbital, and lattice degrees of freedom at temperatures ($T$) below the tetragonal-to-orthorhombic structural phase transition occurring at $T=T_{s}$ Fernandes14 . Particularly, the origin of nematicity in FeSe remains controversial at present AnnaAndreas .

The growing ubiquity of nematic correlated electronic systems, that are scarcely free from impurities, calls for resolving the influence of disorder on the emergence of nematicity. In fact, the strong relevance of disorder to the nematic ordering is also supported from a notable number of experiments detecting local C4-symmetry breaking around impurities chuang10 ; song11 ; zhou11 ; grothe12 ; allan13 ; Rosenthal14 ; kostin . While some of these results may be attributable to, for instance, residual sample strain which explicitly breaks the C4 symmetry locally Ren15 ; Rosenthal14 ; Baek16 , or the presence of stripe-ordered antiferromagnetism chuang10 , the possible pinning of nematic fluctuations due to the presence of disorder appears as a promising and, at the moment, poorly-explored alternative chen09 ; Navarro_cigar ; kontani12 ; Gastiasoro14 ; YWang15 . Even more, there are strong indications for disorder-pinned static local nematicity in the bulk tetragonal phase, i.e. above $T_{s}$ Kasahara12 ; Iye15 ; Zhou16 ; Wiecki17 ; Wang17 ; Toyoda18 . For example, two recent NMR experiments on FeSe Wiecki17 ; Wang17 , found a clear splitting and broadening of the NMR lineshape above $T_{s}$. The presence of short-range nematic order above the bulk $T_{s}$ in FeSe has also been inferred from ARPES and optical-pump conductivity measurements Zhang15 ; Luo17 . Finally, two very recent pair distribution function (PDF) measurements of FeSe found clear evidence of pronounced local orthorhombicity at the length scale of a few nm well above $T_{s}$ koch19 ; frandsen19 , thus providing additional evidence for disorder-induced local nematicity in these systems.

In this Letter, we perform a detailed theoretical study of the role of disorder in systems with D4h point-group symmetry, which additionally feature strong fluctuations or ordering in the nematic channel. The emergence of nematicity is reflected in a non-zero field $N$, which transforms according to the B1g irreducible representation (IR) of D4h. We mainly focus on $T$ above the respective $T_{\rm nem}$ (same as $T_{s}$), at which, the spontaneous thermodynamic C4-symmetry breaking takes place. By employing both phenomenological and microscopic approaches, we address the following three questions: 1) Under what circumstances can disorder generate nematicity locally? 2) What is the spatial profile of the resulting nematic-defect structure? 3) How do finite disorder concentrations influence the nematic transition?

Our main results can be summarized as follows. For $T>T_{\rm nem}$, (i) an impurity potential of arbitrary strength with a spatial profile which respects the C4 symmetry, generates a local nematic field $N(\bm{r})$ with a spatial profile belonging to the B1g IR. By transferring to a polar coordinate system $(x,y)\mapsto(r,\phi)$, this yields the spatial profile $N(r,\phi)\propto\cos(2\phi)$. (ii) This further implies, that, a potential with a C4-symmetric profile does not induce net nematicity, i.e., $\int d\bm{r}N(\bm{r})=0$, but local probes may still detect evidence of C4 symmetry breaking. (iii) However, we show that such a potential can still drive a nematic transition already at $T>T_{\rm nem}$, since it modifies the Stoner criterion for the nematic instability. (iv) A C4-symmetry-breaking impurity potential can induce nonzero net nematicity and, thus, stabilize long-range nematic order. For $T<T_{\rm nem}$, an impurity potential with a spatial profile which respects C4 symmetry modifies the bulk nematicity ($N_{\rm B}$) locally, and results in an inhomogeneous nematic field $N(\bm{r})=N_{\rm B}+\delta N(\bm{r})$, with a $\delta N(r,\phi)$ which is generally not proportional to $\cos(2\phi)$.

We first examine the implications of disorder within a continuum Ginzburg-Landau (GL) approach, that allows exposing generic features of the induced nematic field, i.e. independent of the origin of the electronic nematicity. In fact, our GL results also apply to situations where the nematic field originates from the spontaneous mixing of superconducting order parameters belonging to the A1g and B1g IRs Livanas ; Fernandes-Millis . However, there, one has to further include the possible influence of disorder on the pairing.

The free energy density ${\cal F}(\bm{r})$ is a functional of $N(\bm{r})$ and the disorder potential $V(\bm{r})$. Its invariance under D4h point group operations and time reversal, leads to:

[TABLE]

with $\alpha,\beta,c>0$. Here, we restricted to the lowest-order possible coupling between $V(\bm{r})$ and $N(\bm{r})$. Later on, we consider effects of higher-order terms. For further details on the GL approach, we refer to the Supplementary Material (SM) Suppl . From Eq. (1), one observes that the nematic field couples to the second derivatives of the disorder potential and, thus, to a particular linear combination of the electric field gradients (EFGs). The nematic field is proportional to the quadrupolar electronic charge density defined as Q_{x^{2}-y^{2}}(\bm{r})=\big{(}x^{2}-y^{2}\big{)}\rho(\bm{r}), which transforms according to the B1g IR of D4h, i.e. similar to $N(\bm{r})$. In the above, $\rho(\bm{r})$ denotes the electric charge density, which belongs to the trivial (A1g) IR of D4h. The appearance of a nonzero $N(\bm{r})$, solely due to the presence of disorder, is a consequence of the broken translational invariance, and can be viewed as a result of linear response, since the EFG $\left(\partial_{x}^{2}-\partial_{y}^{2}\right)V(\bm{r})$ acts as a quadrupolar source field, which leads to a nonzero and necessarily inhomogeneous $Q_{x^{2}-y^{2}}(\bm{r})$ and thus $N(\bm{r})$.

For the remainder, we consider $T>T_{\rm nem}$ (unless explicitly stated), which implies that the system resides in the C4-symmetric phase in the absence of disorder. In this case, we can drop the quartic nematic term, since $N(\bm{r})$ is generally small. Thus, for $T>T_{\rm nem}$, the Euler-Lagrange equation of motion (EOM) for Eq. (1) reads:

[TABLE]

The above EOM provides the proportionality relation between the EFG and the resulting nematic field, i.e.:

[TABLE]

where we introduced the nematic coherence length in the tetragonal phase $\xi_{\rm nem}^{-1}=\sqrt{\alpha(T-T_{\rm nem})/c}$.

For a C4-symmetric impurity potential we integrate the angular part of the rhs in Eq. (3), and find the earlier-announced spatial profile $N(r,\phi)\propto\cos(2\phi)$. This profile decays away from the impurity within a range given by $\xi_{\rm nem}$, and this length scale diverges as $T\rightarrow T_{\rm nem}$. Both results are depicted in Fig. 1. The angular dependence transforms exactly according to the B1g IR of D4h. This constraint on the spatial profile of $N(\bm{r})$ is a directly consequence of the featureless (A1g) nature of the disorder potential $V$ itself. Thus, the net electronic nematicity and quadrupolar charge are zero, since:

[TABLE]

Nonetheless, probes like NMR and PDF pick up a signal from atoms in the lobes of the induced $N(\bm{r})$, and do therefore detect clear evidence for local nematicity and orthorhombicity even though global effects are absent.

Equation (4) also reveals that the linear coupling of the nematic field to the EFG cannot stabilize a net thermodynamic nematicity which emerges when $N(\bm{q}=\bm{0})\neq 0$. Therefore, the quadrupolar coupling can neither preempt nor smear out the bulk nematic phase transition. A nonzero $N(\bm{q}=\bm{0})$ can, however, be induced when the spatial profile of the disorder potential explicitly breaks C4 symmetry. This can be seen by including higher-order couplings in the GL free energy (see also SM Suppl ):

[TABLE]

The above terms provide couplings between $V(\bm{q}\neq\bm{0})$ and $N(\bm{q}=\bm{0})$, as well as the $N(\bm{q}\neq\bm{0})$ nematic-field components. These couplings are essential to describe a disorder-driven preemptive nematic transition above $T_{\rm nem}$, as well as the emergence of net nematicity when the potential breaks C4 symmetry. To demonstrate both aspects, we derive the modified EOM for the $N(\bm{q}=\bm{0})$ component after adding the contribution of Eq. (5) to the free energy of Eq. (1). We find the following EOM:

[TABLE]

Thus, a nonzero $N(\bm{q}=\bm{0})$ can only emerge when components with $\bm{q}\neq\bm{0}$ are already nonzero. By assuming that the potential is weak, the $N(\bm{q}\neq\bm{0})$ components remain given by the Fourier transform of Eq. (2). Therefore, we obtain the following equation up to second order in $V(\bm{r})$:

[TABLE]

Eq. (7) implies that a C4-symmetric configuration of impurities cannot source a homogeneous component for the nematic field, since the rhs is zero. As we prove in SM Suppl , this holds even after including all the symmetry-allowed higher-order GL terms. In fact, this result is also recovered in the case of a large number of randomly-distributed and uncorrelated impurities, in which situation, translational and rotational symmetries are preserved on average. Thus, a C4-symmetric disorder potential solely modifies the nematic Stoner criterion, i.e.:

[TABLE]

Depending on the microscopic details which control the sign of the coupling constant $g^{\prime\prime}$, the nematic transition temperature can be enhanced. Note, however, that such an enhancement tends to zero in the thermodynamic limit, unless a critical density of impurities $n_{\rm imp}$ is present. This is because the $g^{\prime\prime}$ coefficient is inversely proportional to the system size. Interestingly, a detailed transport study with controlled disorder by electron irradiation found cases where the critical nematic transition temperature increased slightly with disorder Timmons2019 .

Before proceeding, we point out that the first term of Eq. (5) also allows to describe the induced net nematicity when the disorder potential breaks C4 symmetry. To exemplify this, we consider a dimer impurity potential $V(\bm{r})=V\left[\delta(\bm{r}-\hat{\bm{x}})+\delta(\bm{r}+\hat{\bm{x}})\right]$, which yields $V(\bm{p})=V(\cos p_{x}+\cos p_{y})+V(\cos p_{x}-\cos p_{y})$, for a lattice constant $a=1$. The breaking of C4 symmetry is ensured by the combined presence of the A1g and B1g IRs. In general, a nonzero $N(\bm{q}=\bm{0})$ arises whenever $|V(\bm{q})|^{2}$ contains at least one B1g term.

To support the above GL findings, we employ a microscopic tight-binding model of electrons coupled to disorder. This analysis not only verifies the above GL results, but more importantly, allows to uncover further microscopic details which control the emergence of nematicity. In the absence of disorder, the electrons are described by the dispersion $\varepsilon_{\bm{k}}=-2t(\cos k_{x}+\cos k_{y})-\mu$. The spin degree of freedom is neglected throughout this work, since it merely introduces a factor of 2 in all thermodynamic averages. We assume that the electrons feel an attractive effective interaction in the Pomeranchuk nematic channel as in Ref. Gallais, , which, after mean-field decoupling, yields the nematic order parameter (for details see SM Suppl ):

[TABLE]

i.e. the lattice analog of $N(\bm{r})$. This introduces a local or global C4-breaking to the electron-hopping matrix elements. In the above, $c_{\bm{R}}$ denotes the annihilation operator of an electron at position $\bm{R}=(n,m)$ of the lattice, with $n,m\in\mathbb{Z}$. In addition, $\hat{\bm{x}}$ ($\hat{\bm{y}}$) corresponds to the unit vector in the $x$ ($y$) direction. The nematic form factor is nonzero for $\bm{\delta}=\pm\hat{\bm{x}},\pm\hat{\bm{y}}$, and reads $f_{\pm\hat{\bm{x}}}=-f_{\pm\hat{\bm{y}}}=1/4$. Disorder is considered in the form of point-like identical impurities. The total mean-field Hamiltonian becomes:

[TABLE]

For a single delta-function impurity potential $V_{\bm{R}}=V\delta_{\bm{R},\bm{0}}$, we evaluate the nematic order parameter in Eq. (9) self-consistently for a fixed electron density $\langle n\rangle$ (see SM Suppl ). The resulting nematic order is displayed in Fig. 2(a), and possesses the same spatial profile as those shown in Fig. 1. In the case of a dimer impurity potential $V_{\bm{R}}=V(\delta_{\bm{R},\hat{\bm{x}}}+\delta_{\bm{R},-\hat{\bm{x}}})$, which explicitly breaks C4 symmetry, we obtain the profile shown in Fig. 2(c). Its Fourier transform, see Fig. 2(d), exhibits $N_{\bm{q}=\bm{0}}\neq 0$, which originates from the rhs of Eq. (7). We stress that, the fact that the induced clover pattern in Fig. 2(a) is directly sourced by the EFG, makes it distinct from other microscopic studies of impurity-induced local order andersen07 ; romer2012 ; Navarro_LiFeAs . There, the impurity potential induces a spontaneous symmetry breaking locally, by means of a local fulfillment of the Stoner criterion, i.e. analogously to Eq. (8).

We proceed by studying the effects of a single impurity for $T<T_{\rm nem}$, where the system resides in the bulk phase with a value $N_{\rm B}$ for the nematic order parameter. In this case, the order parameter assumes the form $N(\bm{r})=N_{\rm B}+\delta N(\bm{r})$, where $\delta N(\bm{r})$ incorporates the spatial variation of the nematic order parameter near the impurity. For a weak impurity potential, we expand the EOM stemming from Eq. (1) up to linear order in $\delta N(\bm{r})$ (see SM Suppl ). We find that $\delta N(\bm{r})$ possesses the spatial profile of Eq. (3), with the difference that the coherence length is now given by $\xi_{\rm nem}^{-1}=\sqrt{2\alpha(T_{\rm nem}-T)/c}$ due to an additional contribution of the quartic term which has to be taken into account for $T<T_{\rm nem}$. From a microscopic calculation, we obtain the spatial profile for the nematic order parameter which is shown in Fig. 2(b), exhibiting an anisotropic local structure which is slightly elongated along the $y$ direction. To lowest order in $V(\bm{r})$, this asymmetry found via the microscopic model can be reproduced in the GL theory by including the first term of Eq. (5). The presence of this term yields $\delta N(\bm{r})\propto f(r)\cos(2\phi)+h(r)N_{\rm B}$, where $f(r)$ and $h(r)$ are decaying functions of the radial coordinate, transforming according to the A1g IR. Note that additional higher order terms, e.g. $\propto V(\bm{r})(\partial_{x}^{2}+\partial_{y}^{2})[N(\bm{r})]^{2}$, can further contribute to this anisotropy by modifying $h(r)$. In general, we find that depending on the sign of the impurity potential, point-like disorder at $T<T_{\rm nem}$ may either locally enhance or decrease the nematic order.

Finally, we verify the possibility of disorder-enhanced $T_{\rm nem}$ within the microscopic model. We assume random and dilute disorder of density $n_{\rm imp}$, that on average preserves the C4 symmetry. Within the 1st order Born approximation Bruus-Flensberg , the quasiparticle lifetime is:

[TABLE]

By use of Eq. (11), we can evaluate the microscopic coefficients which enter the modified Stoner criterion of Eq. (8), brought about by the impurities. Starting from Eq. (9), we find that the self-consistency equation for the $\bm{q}=\bm{0}$ component of the nematic mean-field order parameter, corresponding to net nematicity $N\equiv\sum_{\bm{R}}N_{\bm{R}}/{\cal N}=N_{\bm{q}=\bm{0}}/{\cal N}$, reads

[TABLE]

with $f_{\bm{k}}=\cos k_{x}-\cos k_{y}$ and $n_{F}(\varepsilon)$ the Fermi-Dirac distribution function. Linearizing the rhs with respect to $N$, yields the modified Stoner criterion

[TABLE]

In the absence of disorder, i.e. $\tau_{\bm{k}}\rightarrow\infty$, the integration yields the derivative of the Fermi-Dirac distribution $n_{F}^{\prime}(\varepsilon_{\bm{k}})$. However, for finite $\tau_{\bm{k}}$, each $\bm{k}$ state is broadened, and the density of states (DOS) for the $\bm{k}$ points mainly contributing to the nematic instability may be enhanced. To explore this effect, we numerically calculate the nematic susceptibility $\chi_{\rm nem}^{\rm imp}$ in the presence of disorder, which is identified with the rhs of Eq. (12). Fig. 3 shows how this quantity changes versus the disorder concentration $n_{\rm imp}$, relative to the disorder-free case. For an impurity density of $n_{\rm imp}\approx 5\%$, we obtain the maximal relative enhancement of $\chi_{\rm nem}^{\rm imp}$ leading to a corresponding small enhancement of $T_{\rm nem}$. It is tempting to assign the similar small enhancement of the nematic transition temperature measured experimentally in Ref. Timmons2019, to the effect demonstrated in Fig. 3.

The origin of the enhancement effect shown in Fig. 3 is the presence of a nearby van Hove singularity whose spectral weight can be utilized to boost $\chi_{\rm nem}^{\rm imp}$ in the presence of finite $\tau_{\bm{k}}$. Without favorable DOS conditions, disorder generally suppresses the nematic susceptibility and hence also $T_{\rm nem}$. Such a suppression tendency has been previously found in Ref. Schofield, and is also demonstrated in our SM Suppl . Further, we remark that even in the disorder-free case, the presence of a van Hove singularity is pivotal for the stabilization of an electron nematic phase of the Pomeranchuk type. For more details see Refs. Kee ; Yamase .

In summary, we have elucidated the coupling of nematicity to disorder from both a phenomenological GL approach and microscopic calculations. Importantly, disorder is always locally relevant for inducing nematicity since the EFG $\left(\partial_{x}^{2}-\partial_{y}^{2}\right)V(\bm{r})$ directly acts as a quadrupolar source field for nematicity. This explains the detection of local nematicity/orthorhombicity in experimental probes sensitive to different atomic environments within materials. At the global scale, however, disorder does not generally generate long-range nematicity at $T>T_{\rm nem}$ where the system remains tetragonal. Finite disorder concentrations may, however, under favorable circumstances enhance nematic order.

We thank Andreas Kreisel and Karsten Flensberg for useful discussions. D. S. and B. M. A. acknowledge financial support from the Carlsberg Foundation. P. K. and B. M. A. acknowledge support from the Independent Research Fund Denmark grant number DFF-6108-00096. I. P. is supported by ANR grant “IRONIC” (ANR-15-CE30-0025).

I Ginzburg-Landau Theory: Phenomenological Analysis

Eq. (1) of the main text is obtained by demanding that the free energy density is a real functional transforming according to the trivial irreducible representation (IR) of the ensuing point group. Here, we assume a system with tetragonal and inversion symmetries present, which is described by the D4h point group symmetry. The free energy density transforms according to the A1g IR of D4h and is here also assumed invariant under time reversal.

I.1 Equation of Motion

The equation describing the nematic field is found via the Euler-Lagrange equation of motion (EOM):

[TABLE]

and reads:

[TABLE]

From the above, one notes that if the potential $V(\bm{r})$ is homogeneous, i.e. $V(\bm{r})=V$, the EOM includes only derivatives of $N$ and no other spatially-dependent functions or source terms. Thus, for an infinite (bulk) system $N(\bm{r})=N$. When $T>T_{\rm nem}$, the appearance of nematicity is disfavored and, thus, $N=0$ in the bulk. In contrast, the presence of an inhomogeneous potential functions as a source of nematicity and allows for non-zero inhomogeneous solutions of $N(\bm{r})$.

I.2 Case Study: Single Impurity for $\bm{T>T_{\rm nem}}$

Above $T_{\rm nem}$, we drop the cubic term in the EOM in Eq. (14) above, and obtain Eq. (2) of the manuscript. For a potential satisfying $V(\bm{q})=V(|\bm{q}|)$, we set $q_{x}=q\cos\theta$, $q_{y}=q\sin\theta$, $x=r\cos\phi$ and $y=r\sin\phi$, and find:

[TABLE]

with $J_{2}(z)$ the respective Bessel function of the first kind. One notes the distinctive angular dependence of the spatial profile of the induced nematic order, which is fixed by the B1g IR of $N$, the $A_{\rm 1g}$ IR of $V$, and the fourfold-symmetric impurity profile. We resolve the radial dependence in the case $V(\bm{r})=V/r$, and find:

[TABLE]

where we introduced the modified Bessel and Struve functions, and defined $\gamma=-\pi gV/(2\xi_{\rm nem})$. Notably, the decaying function in the brackets yields $\approx 1/2$ for $r=\xi_{\rm nem}$.

I.3 Non-Induction of Net Nematicity by a C4-Symmetric Potential

In this section we explore whether there exists a term in the Ginzburg-Landau expansion which can induce a nonzero $N(\bm{q}=\bm{0})$ for an impurity-potential profile which preserves C4 symmetry. Consider the most general term:

[TABLE]

where, if $\ell$ is odd, then $n=\ell+2\mathbb{N}$. We fix the spatial profile of $V$ to be C4-symmetric. The above general term can be mapped to two distinct types of couplings:

[TABLE]

The respective equations of motion read:

[TABLE]

We Fourier transform the first equation and find:

[TABLE]

Assuming that the components appearing on the rhs are given by:

[TABLE]

where we set $q_{x}=q\cos\theta$ and $q_{y}=q\sin\theta$, we find that the angular part of the integral is proportional to:

[TABLE]

where we set $\cos\theta_{s}=p_{s,x}/p_{s}$ and $\sin\theta_{s}=p_{s,y}/p_{s}$, with $p_{s}=|\bm{p}_{s}|$. A similar treatment for the remaining equation also yields zero. This result naturally confirms, that, a C4-symmetric spatial profile for the impurity potential cannot lead to net nematicity.

I.4 Case Study: Single Impurity for $\bm{T<T_{\rm nem}}$

In order to explain the elongated clover-like spatial profile induced by the impurity in the bulk nematic phase ($T<T_{\rm nem}$), we need to include higher order terms in the free energy described by Eq. (14) of the SM. To demonstrate how this elongation comes about, it is sufficient to solely retain the first term of Eq. (5) presented in the main text. The modified EOM has the form:

[TABLE]

We separate the nematic field into two parts, i.e. $N(\bm{r})=N_{\rm B}+\delta N(\bm{r})$. Here, $N_{\rm B}$ denotes the value of the bulk nematic order parameter and is given by $\beta N_{\rm B}^{2}=\alpha(T_{\rm nem}-T)$ for $T<T_{\rm nem}$, while $\delta N(\bm{r})$ denotes the contribution stemming from the presence of the impurity. For $|\delta N(\bm{r})|\ll|N_{\rm B}|$ we linearize the above EOM and obtain:

[TABLE]

In the above, we retained the terms which lead to a $\delta N(\bm{r})$ which is linear in terms of the strength of the impurity potential. Within this assumption, we dropped the term $\delta N(\bm{r})V(\bm{r})$ which leads to higher-order contributions with respect to the strength of the impurity potential. In the same line of arguments as the ones leading to Eq. (15), we obtain a constant angular profile superimposed on the usual $\cos(2\phi)$-form:

[TABLE]

with the coherence length being given now by $\xi_{\rm nem}^{-2}=2\alpha(T_{\rm nem}-T)$ due to the contribution of the quartic term of the free energy. In connection to Eq. (16) of the SM, we find that for $V(\bm{r})=V/r$:

[TABLE]

with $\gamma^{\prime}=-\pi g^{\prime}V\xi_{\rm nem}/2$. From the above, one can read off the decaying functions $f(r)$ and $h(r)$ discussed in the main text. This spatial profile does indeed lead to a profile on the same form as the anisotropic induced order in Fig. 2(b) of the main text. Furthermore, note that it is the presence of the nonzero bulk nematic order parameter $N_{\rm B}$, that induces the anisotropy.

II Interaction in the Nematic Channel and Mean-Field Theory Decoupling

We assume the presence of the interaction

[TABLE]

which contributes to the desired nematic channel. In the above, we have introduced:

[TABLE]

where we have defined the form factor $f_{\pm\hat{\bm{x}}}=-f_{\pm\hat{\bm{y}}}=1/4$. Note that the lattice constant has been set to unity. We perform a mean-field decoupling of the interaction in the direct channel by introducing the nematic order parameter N_{\bm{R}}=-V_{\rm nem}\big{<}\widehat{{\cal O}}_{\bm{R}}\big{>}. The latter steps led to Eq. (9) of the main text.

In wavevector space we have $N_{\bm{q}}=\sum_{\bm{R}}N_{\bm{R}}e^{-i\bm{q}\cdot\bm{R}}$ and the complete mean-field Hamiltonian reads:

[TABLE]

with ${\cal N}$ being the number of lattice sites, while the nematic form factor in wavevector space takes the form:

[TABLE]

The mean-field Hamiltonian has to be supplemented with the self-constistency equation for the nematic order parameter, which reads

[TABLE]

where we introduced the full single-particle fermionic Matsubara Green function:

[TABLE]

In the above, $k_{n}=(2n+1)\pi T$ ($k_{\rm B}=1$) and the Matsubara Green function for the free electrons has the form $G^{0}_{\bm{k},k_{n}}=1/(ik_{n}-\varepsilon_{\bm{k}})$. The above construction allows us to employ Dyson’s equation in order to perform an expansion of the rhs of the self-consistency equation with respect to the nematic order parameter and/or the impurity potential.

III Ginzburg-Landau Theory: Microscopic Analysis

Given the above, here we show how the electro-nematic coefficient $g$ relates to the microscopic parameters for the disorder-free microscopic model under consideration. We employ a perturbative expansion by employing the Dyson equation for the full Matsubara Green function which reads:

[TABLE]

where we introduced U_{\bm{q};\bm{k}}=\big{(}V_{\bm{q}}+N_{\bm{q}}f_{\bm{q},\bm{k}}\big{)}/{\cal N}. We obtain the lowest order contribution of $U$ by replacing the full Green function on the rhs by the bare one. We find:

[TABLE]

To facilitate the calculations, we consider the continuum limit of our model and assume spinless single-band electrons with a parabolic dispersion \varepsilon(\bm{k})=E_{F}\big{[}\left(k/k_{F}\right)^{2}-1\big{]} with $\bm{k}=(k_{x},k_{y})$, $k=|\bm{k}|$ and set f(\bm{k})=\big{(}k_{x}^{2}-k_{y}^{2}\big{)}/k_{F}^{2}. The quantity of interest, after taking into account the symmetries of $\varepsilon(\bm{k}),f(\bm{k})$ and restricting up to second order terms in $\bm{q}$, reads:

[TABLE]

IV Self-Consistent Calculation of the Nematic Order Parameter

By means of the microscopic Hamiltonian in Eq. (10) of the main text, we calculate the nematic order parameter self-consistently until we reach an accuracy of $10^{-6}$, while keeping the electron density fixed. The expectation values entering in the order parameter and the electron density are calculated by expressing the fermionic field operators in the diagonal basis of the Hamiltonian $c_{\bm{R}}=\sum_{m}\gamma_{m}\langle m|\bm{R}\rangle$ with the defining equation $\widehat{\cal H}\gamma^{\dagger}_{m}|0\rangle=E_{m}|m\rangle$. This leads to the following simplified expressions for the order parameter, and electron density, respectively:

[TABLE]

V Disorder-Modified Stoner Criterion and the Resulting $\bm{T_{\rm nem}^{\rm imp}}$

In the presence of dilute and uncorrelated identical impurities, the disorder may enhance the $T_{\rm nem}$. This was shown in the main text by investigating the modified nematic Stoner criterion. In Fig. 1 of the SM, we provide additional results for other electron-density values. The electron density is calculated via:

[TABLE]

which recovers its usual form $\langle n\rangle=\sum_{\bm{k}}n_{F}(\varepsilon_{\bm{k}})/\mathcal{N}$ in the disorder-free case, i.e. $\tau_{\bm{k}}\rightarrow\infty$. For these calculations finite size effects are diminishing for $\mathcal{N}\sim 40\times 10^{3}$.

In Fig. 1 we demonstrate two typical situations, in which, $T_{\rm nem}$ becomes either enhanced or reduced. This is reflected in the behavior of the quantity $\delta\chi_{\rm nem}/\chi^{0}_{\rm nem}\equiv(\chi_{\rm nem}^{\rm imp}-\chi^{0}_{\rm nem})/\chi^{0}_{\rm nem}$ which is depicted. We first focus on $n_{\rm imp}$ in the vicinity of $5\%$, i.e. the optimal value discussed in the main text.

For the value $\langle n\rangle=0.51$ of the electron density, the Fermi energy is tuned very near the van Hove singularity (see Figs. 1(a,b)), which constitutes the sweet spot for the development of the nematic order parameter in the absence of disorder, since there, $\chi^{0}_{\rm nem}$ obtains its maximum value. From Fig. 1(c) we find that introducing disorder worsens the tendency of the system to develop a nematic order parameter as reflected in the negative values of $\delta\chi_{\rm nem}/\chi^{0}_{\rm nem}$. The addition of disorder broadens the spectral function, and the density of states (DOS) unavoidably becomes lowered, since contributions from low DOS $\bm{k}$ points are taken into account. In contrast, in the case $\langle n\rangle=0.53$ discussed in the main text, and also shown here, the broadening allows the DOS to increase by picking up contributions from the van Hove singularity, while at the same time avoiding significant contributions from other low DOS $\bm{k}$ points. Increasing the electron density to $\langle n\rangle=0.55$ shifts the Fermi level further away from the van Hove singularity and thus reduces its favorable impact on the DOS. As a result, the nematic susceptibility drops and $\delta\chi_{\rm nem}/\chi^{0}_{\rm nem}$ is negative.

The balance between the contributions to the DOS originating from the van Hove singularity and the low DOS $\bm{k}$ points is controlled by the concentration of impurities. Varying $n_{\rm imp}$ leads to a modification of the relative strength of the two competing contributions and allows the sign changes of $\delta\chi_{\rm nem}/\chi^{0}_{\rm nem}$ which are shown in Fig. 1(c) for $\langle n\rangle=0.55$.