Enhanced nematic fluctuations near an antiferromagnetic Mott insulator and possible application to high-$T_{c}$ cuprates

TL;DR

This paper explores how nematic fluctuations are enhanced near a Mott insulator with antiferromagnetic order, suggesting a link to nematicity observed in underdoped cuprate superconductors, through strong-coupling models and simulations.

Contribution

It introduces a strong-coupling model including Cu-$d$ and O-$p$ orbitals that reveals biquadratic spin interactions leading to nematic fluctuations, a novel mechanism for nematicity in cuprates.

Findings

Biquadratic coupling favors stripe-like magnetic fluctuations.

Enhanced nematic susceptibility appears at a temperature scale set by $J$.

Implications for magnetic spectrum and nematicity in cuprates.

Abstract

Motivated by the widespread experimental observations of nematicity in strongly underdoped cuprate superconductors, we investigate the possibility of enhanced nematic fluctuations in the vicinity of a Mott insulator that displays N\'eel-type antiferromagnetic order. By performing a strong-coupling expansion of an effective model that contains both Cu-$d$ and O-$p$ orbitals on the square lattice, we demonstrate that quadrupolar fluctuations in the $p$-orbitals inevitably generate a biquadratic coupling between the spins of the $d$-orbitals. The key point revealed by our classical Monte Carlo simulations and large-$N$ calculations is that the biquadratic term favors local stripe-like magnetic fluctuations, which result in an enhanced nematic susceptibility that onsets at a temperature scale determined by the effective Heisenberg exchange $J$. We discuss the impact of this type of nematic order on the magnetic spectrum and outline possible implications on our understanding of nematicity in the cuprates.

Figures

Equations

$$H_{t-J}=\sum _{ij\alpha }{t}_{ij}{\tilde{d}}_{i\alpha }^{†}{\tilde{d}}_{j\alpha }+J\sum _{⟨ij⟩}\left({\mathbf{S}}_i\cdot {\mathbf{S}}_j-\frac{1}{4}{n}_i{n}_j\right).$$

$$H_K=-K\sum _i{\left[{\mathbf{S}}_i\cdot \left({\mathbf{S}}_{i-\hat{x}}+{\mathbf{S}}_{i+\hat{x}}-{\mathbf{S}}_{i-\hat{y}}-{\mathbf{S}}_{i+\hat{y}}\right)\right]}^2,$$

$$H_{J^{\prime }}=-J^{\prime }\sum _{i,\delta }{n}_{i+\frac{\delta }{2}}^p{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+\delta }$$

$$K=\frac{J^{\prime 2}}{2}\frac{{\int }_{\mathbf{k}}{\Pi }_{\mathbf{k}}^{\eta }}{1-(V_{pp}-\frac{U_{pp}}{8}){\Pi }_{\mathbf{k}=0}^{\eta }}>0.$$

$${\Pi }_{\mathbf{k}}^{\eta }=-{\int }_{\mathbf{q}}\sum _{i,j=1}^4\frac{v_{ij;\mathbf{q},\mathbf{q}+\mathbf{k}}\{n_F({\xi }_{i,\mathbf{q}})-n_F({\xi }_{j,\mathbf{q}+\mathbf{k}})\}}{{ϵ}_{i,\mathbf{q}}-{ϵ}_{j,\mathbf{q}+\mathbf{k}}}$$

$${\phi }_i={\mathbf{S}}_i\cdot \left({\mathbf{S}}_{i-\hat{x}}+{\mathbf{S}}_{i+\hat{x}}-{\mathbf{S}}_{i-\hat{y}}-{\mathbf{S}}_{i+\hat{y}}\right)$$

$${\chi }_{\mathrm{nem}}\left(T\right)={\int }_0^{1/T}d\tau \sum _i⟨{\mathcal{T}}_{\tau }{\phi }_i\left(\tau \right){\phi }_0\left(0\right)⟩$$

$$S$$

$$+{\int }_r\left[r_0{\mathbf{n}}_r^2+\frac{u}{2}{\left({\mathbf{n}}_r\cdot {\mathbf{n}}_r\right)}^2+\frac{{\phi }_r^2}{2g}-h_r{\phi }_r\right]$$

$${\chi }_{\mathrm{nem}}=\frac{{\chi }_{\mathrm{nem}}^{\left(0\right)}}{1-\frac{g}{N}{\chi }_{\mathrm{nem}}^{\left(0\right)}},$$

$${\chi }_{\mathrm{AFM}}\left(\mathbf{Q}+\mathbf{q},\omega \right)=\frac{1}{{\xi }^{-2}+{\mathbf{q}}^2-\phi \left(q_x^2-q_y^2\right)+f\left({\omega }_n\right)},$$

$$H_K|\mathrm{N}\acute{\mathrm{e}}\mathrm{el}⟩=0.$$

$$H_0$$

$$\displaystyle+t_{pp}\sum_{\langle\langle i,j\rangle\rangle}\bigl{(}(-1)^{u^{\prime}_{ij}}p^{\dagger}_{i\sigma}p_{j\sigma}+\text{h.c.}\bigr{)}$$

$$H_U$$

$$H_V$$

$$H_{U_{pp}}+H_{V_{pp}}=-\sum _{\mathbf{k};k_x>0}{\nu }_{\mathbf{k}}^{†}{U}_{\mathbf{k}}^{-1}{\nu }_{\mathbf{k}}$$

$$U_{\mathbf{k}}^{-1}=\frac{2}{N_L}\left(\begin{array}{cc}-\text{Re}{U}_{+,\mathbf{k}}& \frac{i{V}_{pp}}{4}\text{Im}{f}_{\mathbf{k}}\\ -\frac{i{V}_{pp}}{4}\text{Im}{f}_{\mathbf{k}}& \text{Re}{U}_{-,\mathbf{k}}\end{array}\right).$$

$$S$$

$$G_0^{-1}(Q,Q^{\prime })=(i{q}_n-\Delta +\mu )\text{openone}-t_{pp}\text{Re}(h_{\mathbf{q}}){\tau }^x+t_{pp}\text{Im}(h_{\mathbf{q}}){\tau }^y$$

$$[G_{\Phi }^{-1}+G_{J^{\prime }}^{-1}](Q,Q^{\prime })=A_0(q-q^{\prime }){\tau }^0+A_z(q-q^{\prime }){\tau }^z.$$

$$A_0(q-q^{\prime })$$

$$A_z(q-q^{\prime })$$

$$S={\int }_q{\Phi }_q^{†}{U}_q{\Phi }_q-\text{Tr}\ln (-G^{-1}),$$

$$S_2=-\frac{1}{2}{\int }_q\sum _{\alpha ,{\alpha }^{\prime }}{A}_{\alpha }(q)A_{{\alpha }^{\prime }}(-q){\Pi }_{-q}^{\alpha {\alpha }^{\prime }},$$

$$\Pi^{\alpha\alpha^{\prime}}_{q}=-\int_{k}\text{Tr}\bigl{[}G_{0}(k)\tau^{\alpha}G_{0}(k+q)\tau^{\alpha^{\prime}}\bigr{]}\,.$$

$$K_0=\frac{J^{\prime 2}}{2}{\int }_{\mathbf{q}}{\Pi }_{\mathbf{q}}^{zz}.$$

$$\tilde{\Pi}^{\alpha\alpha^{\prime}}_{{\bm{q}}}=\Pi^{\alpha\alpha^{\prime}}_{{\bm{q}}}+\frac{1}{2}\sum_{\beta,\beta^{\prime}}\bigl{(}\tilde{U}^{-1}_{{\bm{q}}}\bigr{)}_{\beta\beta^{\prime}}\Pi^{\beta\alpha}_{\bm{q}}\Pi^{\beta^{\prime}\alpha^{\prime}}_{{\bm{q}}}$$

$$S$$

$$+\frac{{\gamma }^{2}u}{2N}{\int }_{q_1,q_2,q_3}{M}_{q_1}^{\alpha }{M}_{q_2}^{\alpha }{M}_{q_3}^{\beta }{M}_{-q_1-q_2-q_3}^{\beta }+{\int }_r\frac{N{\phi }_r}{2g}$$

$$\chi_{\text{nem}}=\frac{T}{L^{2}}\frac{\partial^{2}Z}{\partial h_{\varphi}^{2}}\Bigr{|}_{h_{\varphi}=0}=\frac{\chi_{\text{nem},0}}{1-\frac{g}{N}\chi_{\text{nem},0}}$$

$${\chi }_{\text{nem},0}=\frac{N}{g}-\frac{T}{L^2⟨{\bar{\phi }}_r^2⟩},$$

$$s$$

Full text

Introduction

Hole-doped cuprates are susceptible to a variety of different types of electronic order in the underdoped regime. Examples include tendencies towards charge order Wu-Nature-2011 ; Ghiringhelli-2012 ; Chang2012 ; LeBoeuf2012 , which becomes long-ranged in the presence of large magnetic fields Wu-Nature-2011 ; LeBoeuf2012 , and tendencies towards nematic order Ando2002 ; Hinkov-Science-2008 ; Daou-Nature-2010 ; Lawler-Nature-2010 ; Cyr-PRB-2015 ; Ramshaw-npj-2017 , characterized by the breaking of the tetragonal symmetry of the system Kivelson-RMP-2003 ; Vojta-AdvPhys-2009 . The fact that these tendencies appear in the region of the phase diagram where a pseudogap is also observed (see schematic Fig. 1) suggest a close interplay between these seemingly different phenomena, a topic that remains widely debated in the field (for recent reviews, see KeimerKivelson-Nature-2015 ; Fradkin2015 ).

Although the microscopic mechanisms behind these different ordering tendencies, and particularly of nematicity, remain unsettled, they have been the subject of many different theoretical proposals (see, for instance, Kivelson-Nature-98 ; Yamase-JPhysSocJpn-2000 ; Kivelson-PRB-2004 ; Yamase-PRB-2006 ; Yamase-PRB-2009 ; Okamoto-PRB-2010 ; Fischer-PRB-2011 ; Andersen-EPL-2012 ; Bulut-PRB-2013 ; Fischer-NJP-2014 ; Volkov-PRB-2016 ; Yamase-PRB-2009 ; Chubukov-PRB-2014 ; Schuett-PRL-2015 ; Nie-PRB-2017 ; Scheurer-PRL-17 ; Tsuchiizu-PRB-2018 , and also the reviews Kivelson-RMP-2003 ; Vojta-AdvPhys-2009 ). While a complete theory for nematicity in the cuprates is beyond the scope of our work, here we show that an important contribution to the nematic susceptibility arises already near the Mott (or more precisely, charge-transfer Zaanen-PRL-1985 ) insulating state of the parent compound. For the rest of the paper, thus, we focus only on the spin correlations near the Mott state, and neglect other phenomena that are certainly important for a complete description of the hole-doped cuprates, and which may also be important to describe nematicity, such as charge order, pseudogap, time-reversal symmetry-breaking, pair-density waves, and superconductivity KeimerKivelson-Nature-2015 ; Fradkin2015 .

To be more specific, we consider the so-called Emery model Emery-PRL-1987 , an effective model that attempts to capture both Cu and O low-energy degrees of freedom by introducing $d_{x^{2}-y^{2}}$ orbitals on the sites of the square lattice and $p_{x}$ ($p_{y}$) orbitals on the horizontal (vertical) bonds. Consider first the case where only $d$-orbitals are present. In the half-filled Mott insulating state, the charge degrees of freedom are quenched, and the low-energy physics is described completely in terms of an AFM Heisenberg interaction $J$ between the $d$-orbital spins, which ultimately gives rise to a Néel AFM ground state. Upon light hole-doping, the effective Hamiltonian is known as the $t-J$ model Lee-RMP-2006 :

[TABLE]

Here, $t_{ij}$ denotes the hole hopping parameters and $J$ the AFM exchange coupling. The operator $\mathbf{S}_{i}=\frac{1}{2}\sum_{\alpha\beta}\tilde{d}_{i\alpha}^{\dagger}\bm{\sigma}_{\alpha\beta}\tilde{d}_{i\beta}$ describes the $d$-orbital spin and $n_{i}=\sum_{\alpha}\tilde{d}_{i\alpha}^{\dagger}\tilde{d}_{i\alpha}$ the corresponding charge. The strong local Coulomb interaction is incorporated in terms of the hole creation operator $\tilde{d}_{i\alpha}^{\dagger}=\left(1-n_{i\bar{\alpha}}\right)d_{i\alpha}^{\dagger}$, reflecting the fact that double occupancy of the sites is not allowed near the Mott insulating state.

As we demonstrate below via a strong coupling expansion of the Emery model, the inclusion of the $p$-orbitals leads to an important additional term in the $t-J$ Hamiltonian. While the two terms in Eq. (1) remain the same, albeit with a different microscopic expression for $J$, non-critical quadrupolar fluctuations of the $p$-orbitals, enhanced by the repulsion between $p$-orbitals, generate a positive biquadratic coupling $K>0$ between the $d$-orbital spins:

[TABLE]

resulting in an effective $t-J-K$ Hamiltonian, $H_{t-J-K}=H_{t-J}+H_{K}$.

Using classical Monte Carlo and large-$N$ analytical methods, we find that the main consequence of $H_{K}$ is to enhance the static electronic nematic susceptibility $\chi_{\mathrm{nem}}$ near the AFM-Mott insulating state. However, $\chi_{\mathrm{nem}}$ is not found to diverge on its own – instead, it peaks at a temperature scale proportional to $J$, instead of $K$. The location of the peak depends on the relative strength of quantum and thermal fluctuations and shifts towards smaller temperatures for larger quantum fluctuations. As illustrated in Fig. 1(c), the enhancement of nematic fluctuations promoted by $H_{K}$ has its origins on the short-ranged magnetic stripe ordered regions that this term favors within the (much longer ranged) Néel ordered background. Consequently, within the $t-J-K$ model, the onset of nematic order requires an additional symmetry breaking field that can take advantage of the enhanced susceptibility. While a more detailed discussion of the application of these results to the cuprates is left to the end of this paper, we note that this mechanism for enhanced nematic susceptibility can in principle be combined with other mechanisms proposed in the literature to yield long-range nematic order. Detailed reviews on the proposed mechanisms for nematicity in cuprates, both within weak and strong coupling regimes, can be found for instance in Refs. Kivelson-RMP-2003 ; Vojta-AdvPhys-2009 .

Results

Microscopic model. Our starting point is the interacting three-orbital Emery model $H=H_{0}+H_{U}+H_{V}$ Emery-PRL-1987 . As depicted in Fig. 1(b), it includes the $d_{x^{2}-y^{2}}$ Cu orbital with creation operator $d_{i,\sigma}^{{\dagger}}$ at Bravais lattice position ${\bm{R}}_{i}$ and spin $\sigma$ as well as the $p_{x}$ and $p_{y}$ O orbitals with creation operators $p_{i+\frac{\hat{x}}{2},\sigma}^{{\dagger}}$ and $p_{i+\frac{\hat{y}}{2},\sigma}^{{\dagger}}$. The non-interacting part $H_{0}$ includes hopping between $p$-orbitals with (amplitude $t_{pp}$) and between $d$- and $p$-orbitals (with amplitude $t_{pd}$). The corresponding sign factors of the hopping elements follow from the phases of the orbitals (see Fig. 1(b)) Emery-PRL-1987 . In addition, $H_{0}$ contains on-site terms where the energy difference between Cu and O orbitals is given as $\Delta=\varepsilon_{p}-\varepsilon_{d}$. Interactions are included on-site $H_{U}=U_{dd}\sum_{i}n_{i,\uparrow}^{d}n_{i,\downarrow}^{d}+\frac{U_{pp}}{2}\sum_{i,u}n_{i+u,\uparrow}^{p}n_{i+u,\downarrow}^{p}$ with $u\in\{\frac{\hat{x}}{2},\frac{\hat{y}}{2}\}$, and number operators $n_{i,\sigma}^{d}=d_{i,\sigma}^{{\dagger}}d_{i,\sigma}$ and $n_{i+u,\sigma}^{p}=p_{i+u,\sigma}^{{\dagger}}p_{i+u,\sigma}$. We also consider nearest-neighbor interactions $H_{V}=\frac{V_{pp}}{2}\sum_{i,u,u^{\prime}}n_{i+u}^{p}n_{i+u+u^{\prime}}^{p}+V_{pd}\sum_{i,u}n_{i}^{d}n_{i+u}^{p}$ with $u^{\prime}\in\{\pm\frac{1}{2}(\hat{x}+\hat{y}),\pm\frac{1}{2}(\hat{x}-\hat{y})\}$ and $n_{i}^{d}=\sum_{\sigma}n_{i,\sigma}^{d}$, $n_{i+u}^{p}=\sum_{\sigma}n_{i+u,\sigma}^{p}$.

The largest energy scales are the local repulsion $U_{dd}$ between $d$-orbitals and the charge-transfer energy $\Delta$ (with $U_{dd}$ much larger than $\Delta$), suggesting a strong coupling expansion in small $t_{ij}\ll\Delta,U_{dd}-\Delta$. This yields a description in terms of localized $d$-orbital spins ${\bm{S}}_{i}$ coupled to mobile $p$-orbital holes. An expansion up to fourth order in the hopping term $t_{pd}$ was performed in Ref. Zaanen-PRB-1988 ; Kolley-JPhysC-1992 . There appear Kondo-like exchange couplings $\propto{\bm{S}}_{i}\cdot{\bm{s}}_{i+u_{1},i+u_{2}}$ between the $d$- and $p$-orbital spin-densities ${\bm{s}}_{i+u_{1},i+u_{2}}=\frac{1}{2}\sum_{\tau,\tau^{\prime}}p_{i+u_{1},\tau}^{\dagger}\bm{\sigma}_{\tau\tau^{\prime}}p_{i+u_{2},\tau^{\prime}}$ ZhangRice-PRB-1988 , the familiar Heisenberg spin exchange term $J\sum_{\langle i,j\rangle}{\bm{S}}_{i}\cdot{\bm{S}}_{j}$ and terms that renormalize the $p$-orbital hole dispersion. For details we refer to the Methods section and the Supplementary Information. The Kondo-like terms also modify the hole dispersion as the tunneling process of holes through a $d$-orbital becomes spin dependent. For example, tunneling through a background of Néel ordered $d$-orbital spins leads to the spin dependent hopping parameters $t_{a}=\frac{t_{pd}^{2}}{2}(\frac{1}{\Delta}+\frac{3}{U_{dd}-\Delta})$ and $t_{b}=\frac{t_{pd}^{2}}{2}(\frac{3}{\Delta}+\frac{1}{U_{dd}-\Delta})$ for holes with spin parallel and antiparallel to the central $d$-orbital spin Fischer-NJP-2014 . The hole Fermi surface thus appears at momenta ${\bm{k}}=(\pm\frac{\pi}{2},\pm\frac{\pi}{2})$ for small doping $n_{p}$.

Most notably for our considerations, the strong coupling expansion also yields a spin exchange term that depends on the occupation of the intermediate $p$-orbital between $d$-orbital sites:

[TABLE]

where $\delta\in\{\pm\hat{x},\pm\hat{y}\}$ and the spin exchange coupling constant is given by $J^{\prime}=\sum_{n=0}^{3}\frac{t_{pd}^{4}{\rm sign}\left(3-2n\right)}{\Delta^{3-n}\left(U_{dd}-\Delta\right)^{n}}$. Note that $J=\sum_{n=0}^{2}\frac{t_{pd}^{4}(4-n^{2}-\delta_{n,2})}{2\Delta^{3-n}(U_{dd}-\Delta)^{n}}$ and in the large-$U_{dd}$ limit both are of the same order $J^{\prime}/J\rightarrow 1/2$. Oxygen charge fluctuations thus not only renormalize the Heisenberg exchange via the Kondo coupling terms, but, as we show now, also lead to the biquadratic spin exchange interaction $K$ in Eq. (2).

We derive the biquadratic exchange $K$ by first decomposing the $p$-orbital densities as $n_{i+\frac{\hat{x}}{2}}^{p}=n_{i}^{p}+\eta_{i}$ and $n_{i+\frac{\hat{y}}{2}}^{p}=n_{i}^{p}-\eta_{i}$, where $\eta_{i}$ is the quadrupolar (nematic) component of the oxygen charge density Fischer-PRB-2011 . The combination of on-site and nearest-neighbor Coulomb interactions between $p$-orbital leads to a term $(-2\sum U_{{\bm{k}}}\eta_{{\bm{k}}}\eta_{-{\bm{k}}})$ in the Hamiltonian, where $U_{{\bm{k}}}=\frac{1}{4}(V_{pp}\text{Re}f_{{\bm{k}}}-\frac{U_{pp}}{2})$ and $f_{{\bm{k}}}=1+e^{-ik_{x}}+e^{ik_{y}}+e^{i(k_{y}-k_{x})}$. Integrating out the quadrupolar charge fluctuations associated with the $p$-orbitals (details of this analysis are presented in the Methods section and the Supplementary Information) yields the result for the biquadratic exchange interaction in Eq. (2) with:

[TABLE]

Here, $\Pi_{\mathbf{k}}^{\eta}=-\int_{\mathbf{q},\omega^{\prime}}\text{Tr}[G_{\mathbf{\mathbf{q},\omega^{\prime}}}^{p}(\tau^{z}\sigma^{0}s^{0})G_{\mathbf{q}+\mathbf{k},\omega^{\prime}}(\tau^{z}\sigma^{0}s^{0})]$ is the bare $p$-orbital charge susceptibility in the quadrupolar (i.e. nematic) channel We have used the long-wavelength approximation in the denominator for simplicity (the full expression can be found in the Supplementary Information) and yields qualitatively identical results). The Pauli matrices $\tau^{j}$, $s^{j}$ and $\sigma^{i}$ act in orbital $(p_{x},p_{y})$, spin and reduced wavevector \bigl{(}{\bm{k}},{\bm{k}}+(\pi,\pi)\bigr{)} space, respectively. Note that the presence of the AF background of $d$-orbital spins leads to a doubling of the unit cell, and thus $\int_{{\bm{q}}}\equiv\int_{\text{mBZ}}\frac{d^{2}q}{2\pi^{2}}$ is an integration over the magnetic Brillouin zone (mBZ). Explicit expressions for the $p$-orbital Green’s functions $G_{\mathbf{q},\omega}^{p}$ are given in the Methods section and the Supplementary Information and yield

[TABLE]

with $\epsilon_{i,{\bm{q}}}$ being the renormalized $p$-orbital dispersion and $\xi_{i,{\bm{q}}}=\epsilon_{i,{\bm{q}}}-\mu$ with the chemical potential $\mu$. The matrix elements $v_{ij;{\bm{q}},{\bm{q}}+{\bm{k}}}=U_{{\bm{q}}}^{{\dagger}}(\tau^{z}\sigma^{0})U_{{\bm{q}}-{\bm{k}}}$ contain $U_{{\bm{q}}}$, which are unitary matrices that transform between orbital/reduced ${\bm{k}}$-space ($\sigma^{i}\otimes\tau^{j})$ and band space. Most importantly, Eq. (4) makes it clear that the biquadratic exchange $K$ in Eq. (2) is a direct consequence of quadrupolar oxygen charge fluctuations. These are a generic feature of the model and exist even if the $p$-orbital holes are not dressed by $d$-orbital spins.

The oxygen quadrupolar susceptibility $\Pi_{{\bm{k}}}^{\eta}$ (and thus $K$) is strictly positive for all ${\bm{k}}$ and is determined by the occupation number difference between the different oxygen bands. In the relevant regime of small hole fillings $n_{p}\ll 1$, the response approaches a value $\Pi_{{\bm{k}}=0}^{\eta}\propto n_{p}$ at low $T$, peaks around $T\approx|\mu|$ and vanishes as $1/T$ at large $T$. The response increases for smaller bandwidth, e.g., smaller $t_{pp}$. This is derived explicitly in the Supplementary Information for a simpler two-band model that neglects the interaction with the AF background. It also holds true numerically for the full four-band model, as shown in Fig. 2(a), where we present results for the renormalized quadrupolar response \tilde{\Pi}_{\bm{k}=0}^{\eta}=\frac{1}{2}\bigl{[}(\Pi_{\mathbf{k}=0}^{\eta}\bigr{)}^{-1}-U_{{\bm{k}}=0}\bigr{]}^{-1} and for the resulting $K/J$ within the microscopic four-band model. In the calculation we keep $n_{d}=1$, assuming that holes are doped into the $p$-orbitals, but we take the interaction of the mobile holes with the AF background of $d$-orbital spins fully into account. We clearly observe that a large nearest-neighbor repulsion $V_{pp}$ and a small bandwidth $t_{pp}$ enhance $K$ (see Eq. (4)). Our results also indicate that an enhancement of the quadrupolar density fluctuations by $V_{pp}$ is necessary for a significant biquadratic exchange coupling. This follows from $K/J\propto(J^{\prime 2}/J)(n_{p}/t_{pp})$ at small $n_{p}$ where $J^{\prime}=0.08$, $J=0.15$ for the parameters in Fig. 2. Finally, while phonon modes in the same channel are, by symmetry, allowed to give rise to similar behavior, the electronic mechanism for biquadratic exchange is expected to be quantitatively much stronger.

Enhanced nematic susceptibility. The implications of $H_{K}$ can be better understood in the limit of $K>J$. In this case, the AFM ground state is no longer the Néel configuration with ordering vector $\mathbf{Q}=\left(\pi,\pi\right)$, but the striped configuration with $\mathbf{Q}=\left(\pi,0\right)$ or $(0,\pi)$. While the limit of large $K/J$ is clearly not realized in the cuprates, it reveals that $H_{K}$ supports quantum and classical fluctuations with local striped-magnetic order that have significant statistical weight.

We qualitatively demonstrate this behavior in Fig. 1(c) by showing typical spin configurations of a Monte Carlo analysis of $H_{t-J-K}$ in the limit of classical spins and where the kinetic energy of the holes is ignored. One clearly sees local striped-magnetic fluctuations (light red background) in an environment of Néel ordered spins (yellow background). Configurations with parallel spins along the $x$-axis and along the $y$-axis occur with equal probability, hence preserving the tetragonal symmetry of the system. If one, however, weakly disturbs tetragonal symmetry, e.g. by straining one of the axes, this balance is disturbed and one favors striped configurations of one type over the other.

The behavior described above can be quantified in terms of the composite spin variable:

[TABLE]

which changes sign under a rotation by $\pi/2$. Note that the square of this term, which appears in Eq.(2), is invariant under this transformation, and therefore is fully consistent with the four-fold symmetry of the Emery model. While $\left\langle\varphi_{i}\right\rangle=0$ for realistic values of $K$ (and in the absence of external strain), the static nematic susceptibility

[TABLE]

is a measure for the increased relevance of local stripe magnetic configurations. Here, $\mathcal{T}_{\tau}$ denotes imaginary time ordering.

We present a quantitative demonstration that the biquadratic exchange $K$ yields an enhanced nematic susceptibility in the $B_{1g}$ ($x^{2}-y^{2}$) symmetry channel in Fig. 2(b-c). It contains Monte Carlo results for $\chi_{{\rm nem}}$ for a collection of classical Heisenberg spins that interact according to the $H_{J-K}$ model. For consistency with the known spin-wave spectrum, we have included an additional small second-neighbor exchange $J_{2}=-0.1J$ in the simulation. One clearly sees that the biquadratic term $K$ enhances the nematic response in the $B_{1g}$ channel, corresponding to an inequivalence between the $x$ and $y$ axes. In the limit of classical spins, the nematic susceptibility $\chi_{{\rm nem}}\left(T\right)$ is non-monotonic, peaking at a temperature governed by the effective exchange interaction of the spins, $T_{{\rm nem}}\sim J$, which is independent on $K$.

The Monte Carlo results also display that $\chi_{{\rm nem}}\left(T\rightarrow 0\right)\rightarrow 0$, which is a consequence of the classical nature of the spins in the simulations and a resulting absence of (thermal) fluctuations in the zero temperature limit. Quantum fluctuations crucially modify this behavior and lead to $\chi_{{\rm nem}}\left(T\rightarrow 0\right)>0$. This is demonstrated in Fig. 3, where we present results of an analytical calculation of the nematic response $\chi_{\text{nem}}$ that includes the effect of quantum fluctuations within a soft-spin field-theoretical version of the spin degrees of freedom in Eq. (2). After decoupling the biquadratic exchange term $K$ in the nematic channel and taking the long wavelength limit, which is appropriate to study the low-energy excitations, we obtain the effective action:

[TABLE]

where $\mathbf{n}_{r}$ is the $N=3$ component staggered Néel order parameter, as used in the non-linear sigma model of Refs. ChakravartyHalperinNelson-PRB-1989 ; Chubukov-PRB-1994 . The parameter $r_{0}$ controls the distance to the AFM Néel quantum critical point located at $r_{0,c}$ For $\delta r_{0}\equiv r_{0}-r_{0,c}<0$, the system has long-range AFM order at $T=0$, whereas for $\delta r_{0}>0$ it is in the paramagnetic phase (see sketches at the bottom of Fig. 3), and the interaction parameters are $g\propto K/J>0$, $u>g$. The integrations are over $\int_{r}\equiv\int_{0}^{1/T}d\tau\int d^{2}r$, where $r=\left(\tau,\mathbf{r}\right)$ combines imaginary time $\tau$ and position $\mathbf{r}=(x,y)$. In addition, $\varphi_{r}$ is the nematic order parameter of Eq. (6) and $h_{r}$ is an external strain field. The quantum dynamics of the Néel order parameter is governed by $S_{{\rm dyn}}=\int_{q}f\left(\omega_{n}\right)\mathbf{n}_{q}\cdot\mathbf{n}_{-q}$, where $f\left(\omega_{n}\right)\propto\omega_{n}^{2}$ at half filling, while $f\left(\omega_{n}\right)=\gamma\left|\omega_{n}\right|$ was proposed to describe particle-hole excitations, and will be used below as we describe the system away from half-filling $n_{p}>0$. Here, $q=\left(\omega_{n},\mathbf{q}\right)$ combines Matsubara frequency $\omega_{n}$ and momentum $\mathbf{q}$ (measured relative to the AFM ordering vector $\mathbf{Q}=\left(\pi,\pi\right)$) and $\int_{q}\equiv T\sum_{n}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}$.

The nematic susceptibility in Eq. (7) can be obtained for general $N$ and reads (see Methods section and Supplementary Information):

[TABLE]

where the bare nematic susceptibility is given by $\chi_{{\rm nem}}^{\left(0\right)}=\frac{N}{2}\int_{q}\frac{|\bm{q}|^{4}\cos^{2}(2\theta)}{\left(\xi^{-2}+|\bm{q}|^{2}+f(\omega_{n})\right)^{2}}$ with ${\bm{q}}=|\bm{q}|(\cos\theta,\sin\theta)$. Here, $\xi$ is the magnetic correlation length for Néel order, that includes interaction corrections and diverges at the AFM phase transition. Right above the quantum critical point at $\delta r_{0}=0$, one finds $\xi^{-2}=a\gamma T$ with non-universal constant $a$. As shown in Fig. 3(a), the exact shape of $\chi_{\text{nem}}^{(0)}(T)$ depends on this non-universal parameter $a$, which depends, for example, on the interaction parameter $u$ or the lattice constant. While $\chi_{\text{nem}}^{(0)}$ peaks at finite temperatures for $a<\pi$, which is similar to the classical case, the maximum occurs at $T=0$ for $a>\pi$. Note that nematic correlations remain finite ranged at the AFM quantum critical point and universal behavior of $\chi_{\text{nem}}^{(0)}$ is not guaranteed (in contrast to the AFM susceptibility, which is universal). Being a non-universal quantity, we thus expect that the precise shape of $\chi_{\text{nem}}^{(0)}$ can be different for different systems.

In order to make analytic progress and calculate $a(u)$, or more generally $\xi(T,\delta r_{0})$, we consider the limit of large $N$. This approach led to important insights in both the description of antiferromagnetic correlations of the cuprate parent compounds Chubukov-PRB-1994 and of nematic fluctuations of iron-based superconductors Fernandes-PRB-2012 . The magnetic correlation length $\xi$ is determined self-consistently within large-$N$ for a given distance to the AFM quantum critical point $\delta r_{0}$. Despite the similarity between Eq. (9) and the expression for the nematic susceptibility of iron-based superconductors Fernandes-PRB-2012 , there are very important differences between the two systems. Because the iron pnictides order magnetically in a striped configuration, $\chi_{{\rm nem}}^{\left(0\right)}$ diverges when $\xi\rightarrow\infty$, which guarantees that a nematic transition takes place already in the paramagnetic state for any $g>0$. However, because our model orders in a Néel configuration, $\chi_{{\rm nem}}^{\left(0\right)}$ remains finite even when $\xi\rightarrow\infty$. Although long-range nematic order is not present, nematic fluctuations can be significantly enhanced if the biquadratic exchange $K\propto g$ is sufficiently large.

In Fig. 3(b,c) we show the nematic susceptibility obtained within the large-$N$ approach (see Methods section and Supplementary Information). Like in the Monte Carlo results (see Fig. 2), we observe, in Fig. 3(b), a broad maximum at finite temperatures around $T\approx J$, corresponding at $\delta r_{0}=0$ to $a<\pi$. The lattice cutoff $\Lambda$ plays the role of $J$ in the continuum model. The effect of $g$, and thus of the biquadratic exchange $K$, is to enhance the amplitude of the peak (comparing dashed and solid lines). The pronounced peak of $\chi_{{\rm nem}}$ originates from the bare susceptibility $\chi_{{\rm nem}}^{\left(0\right)}$. As discussed above, the bare response is in turn governed by the magnetic correlation length $\xi$ that is set by $T/J$. Notably, at low temperatures, quantum fluctuations render $\chi_{{\rm nem}}$ (and $\chi_{\text{nem}}^{(0)}$) finite, in stark contrast to our MC results for classical spins. Keeping $T$ fixed and varying the non-thermal tuning parameter $\delta r_{0}$, we observe in Fig. 3(c) that the nematic response increases for an increasing magnetic correlation length, i.e. Néel fluctuations enhance the nematic susceptibility. This follows from the observation that \chi_{{\rm nem}}$$\left(\delta r_{0},T\right) is an increasing function for decreasing $\delta r_{0}$.

Consequences of long-range nematic order. As we discussed above, the nematic susceptibility does not diverge within our $H_{J-K}$ model. Nevertheless, it is interesting to study what happens to the magnetic spectrum if nematic order is induced – either by the presence of a small tetragonal-symmetry breaking field $h$, which can induce a sizable nematic order parameter $\varphi\approx\chi_{{\rm nem}}h$, or by combination with other microscopic mechanisms for nematicity. From the action in Eq. (8), we can readily obtain the dynamic spin susceptibility in the presence of nematic order

[TABLE]

Therefore, as shown in Fig. 4, non-zero $\varphi$ modifies the spin-spin structure factor near the Néel ordering vector $\mathbf{Q}$ from a circular shape, which preserves tetragonal symmetry, to an elliptical shape, which breaks tetragonal symmetry. In addition, as $\varphi$ increases, it shifts the maximum of $\chi_{{\rm AFM}}\left(\mathbf{Q}+\mathbf{q},\omega\right)$ from the commensurate $\mathbf{q}=0$ value to an incommensurate wavevector $\mathbf{q}_{\mathrm{IC}}\neq 0$, with $\mathbf{q}_{\mathrm{IC}}$ parallel to either the $x$ axis (if $\varphi>0$) or to the $y$ axis (if $\varphi<0$). Note that a somewhat related mechanism for the incommensurate spin order, based on the $t-J$ model, was reported in Refs. Shraiman-PRL-1988 ; Sushkov-PRL-2005 ; Gabay-PhysicaC-1989 . Previous works have also focused on nematicity arising from a pre-existing incommensurability Nie-PRB-2017 , whereas in our scenario incommensurate magnetic order is a consequence of nematic order, caused by an enhanced nematic susceptibility in the presence of an external symmetry breaking field.

While these effects onset in the paramagnetic phase, the presence of nematic order should also be manifested in the Néel ordered state by the direction of the $d$-orbital moments, which would align parallel to either the $x$ axis or to the $y$ axis Nafradi-PRL-2016 . Within our model, we argue that such an effect would arise in the presence of spin-orbit coupling in the $p$-orbitals, which convert an imbalance in the charge of the $p_{x}$ and $p_{y}$ orbitals into a preferred direction for the $d$-orbital moment.

At first sight, one might anticipate that $K$ would affect the spin-wave dispersion of the AFM Néel ground state Coldea-PRL-2001 . As we show in the Supplementary Information, however, the biquadratic exchange of Eq. (2) does not modify the linearized classical spin-wave spectrum. The reason for this peculiar behavior is that the biquadratic exchange annihilates the classical Néel state, i.e. the vacuum of the linear spin wave excitations:

[TABLE]

It is important to point out that all results discussed here were obtained considering that the spins of the $H_{J-K}$ Hamiltonian are treated as vectors, either classical or in the large-$N$ regime. It is interesting to ask what happens if one considers the quantum spin-$1/2$ case. It turns out that, for spin-$1/2$, the biquadratic term $K>0$ transforms into an AFM next-nearest neighbor bilinear exchange coupling, which certainly changes the spin-wave spectrum. It remains an open question whether the results presented here remain unchanged if one performs this transformation from biquadratic to bilinear exchange in the microscopic model. Importantly, however, we note that a large AFM next-nearest neighbor exchange also favors a stripe magnetic state over a Néel state. Thus, the main ingredient that enhances the nematic susceptibility in the classical spin case seems also to be present in the spin-$1/2$ case.

Discussion

In summary, we showed via a strong coupling expansion of the Emery model that quadrupolar charge fluctuations in the $p$-orbitals generate a biquadratic exchange coupling between the $d$-orbital spins, extending the celebrated $t-J$ model employed to describe lightly-doped Mott insulators. The main effect of this biquadratic term is to enhance $B_{1g}$ nematic fluctuations, which however is not translated into a diverging nematic susceptibility. Importantly, the temperature at which the nematic susceptibility peaks is set not by the biquadratic coupling $K$, but by the standard nearest-neighbor exchange coupling $J$. The peak position is controlled by the relative strength of thermal versus quantum fluctuations, and moves from a temperature of order $J$ for dominant thermal fluctuations towards zero for dominant quantum fluctuations. The biquadratic exchange $K$, however, sets the amplitude of the peak, and both increase for larger values of the repulsion $V_{pp}$ between nearest-neighbor $p$-orbitals.

Thus, our main result is that magnetic correlations associated with the Mott insulating state generate an enhanced nematic susceptibility, which is driven by quadrupolar oxygen density fluctuations. In the remainder of this section, we discuss the possible applications of these results to the nematic tendencies observed in hole-doped cuprates Kivelson-RMP-2003 ; Vojta-AdvPhys-2009 . The mechanism discussed here does not lead to long-range nematic order on its own. However, given the enhanced nematic susceptibility, it is expected that a small tetragonal symmetry-breaking field would lead to a sizable nematic order parameter. Such a symmetry-breaking field is naturally provided by the CuO chains or double chains in YBa2Cu3O7-δ and YBa2Cu4O8, respectively. Interestingly, in YBCO, several experimental observations are consistent with the existence of an electronic nematic order parameter Ando2002 ; Hinkov-Science-2008 ; Daou-Nature-2010 ; Ramshaw-npj-2017 . Whether the observed nematicity is the result of the intrinsic small symmetry-breaking field combined with a large nematic susceptibility, or the consequence of true long-range order that would onset even if the chains were absent, remains to be determined.

Still in what concerns YBCO, it is interesting to note that nematic order is observed already at rather small doping levels, in the vicinity of the Mott insulating Néel state Hinkov-Science-2008 ; Haug-NJP-2010 . In this region of the phase diagram, where our results are the most relevant, the experimental nematic onset temperature is comparable to that of the Néel transition temperature, which in turn is set by $J$. Of course, since nematicity is not restricted only to the vicinity of the Néel state, it is possible that there are different mechanisms responsible for nematicity in different regions of the phase diagram Nie-PRB-2017 .

Neutron scattering experiments in YBCO also reveal a strong feedback of nematic order on the magnetic spectrum Hinkov-Science-2008 ; Haug-NJP-2010 ; Nafradi-PRL-2016 . In particular, nematic order is manifested as an elliptical spin structure factor centered at the Néel ordering vector. Upon lowering the temperature, the peak splits and gives rise to two unidirectional incommensurate peaks. These observations are qualitatively consistent with our results for the effect of nematicity on the AFM magnetic spectrum (see Fig. 4).

To further test the applicability of the effect discussed here on the physics of the cuprates, it would be desirable to directly measure the nematic susceptibility in tetragonal cuprates. In analogy to what has been done for the iron-pnictides (see Ref. Fernandes-NatPhys-2014 ), $\chi_{\mathrm{nem}}$ is closely related to several observables, such as the elastoresistance Chu-Science-2012 , the shear modulus, or electronic Raman scattering. If the biquadratic term found here was to govern the nematic properties of tetragonal cuprates, such as HgBa2CuO4, $\chi_{\mathrm{nem}}$ should be enhanced but not divergent – possibly displaying a peak at a temperature comparable to $J$. Furthermore, the temperature dependences of $\chi_{\mathrm{nem}}$ in the $B_{1g}$ and $B_{2g}$ channels would be similar, although the former would be larger.

Because the biquadratic term is the result of charge fluctuations on the oxygen $p$-orbitals, it is only present in the hole-doped side of the phase diagram, since electron-doping adds charge carriers directly to the Cu sites Armitage-RMP-2010 . To the best of our knowledge, nematic tendencies have not been reported in electron-doped cuprates Kivelson-RMP-2003 . It would be interesting to verify this effect by experimentally determining $\chi_{\mathrm{nem}}$ in tetragonal electron-doped systems, such as Nd2CuO4.

Methods

Derivation of $t-J-K$-model from microscopic three-band model

We derive the biquadratic $K$ spin exchange term in Eq. (2) from a microscopic interacting three-band model $H=H_{0}+H_{U}+H_{V}$ that takes Cu $d_{x^{2}-y^{2}}$ and O $p_{x},p_{y}$ orbitals into account and reads

[TABLE]

Here, $d^{\dagger}_{i\sigma}$ creates a hole in the $d_{x^{2}-y^{2}}$ orbital at Bravais lattice site ${\bm{R}}_{i}$, and $p^{\dagger}_{i\sigma}$ creates a hole in the O $p_{x}$ and $p_{y}$ orbital $i\equiv{\bm{R}}_{i}+\frac{\hat{{\bm{x}}}}{2}$ and $i\equiv{\bm{R}}_{i}+\frac{\hat{{\bm{y}}}}{2}$, respectively (see Fig. 1(b)). The parameters in the Hamiltonian are the on-site orbital energy difference $\Delta=\epsilon_{p}-\epsilon_{d}$, hoppings $t_{pp}$, $t_{pd}$ (see Fig. 1), on-site interactions $U_{pp}$, $U_{dd}$ and nearest-neighbor interactions $V_{pp}$, $V_{pd}$.

As $U_{dd}$ is the largest energy scale, we perform a strong coupling expansion which yields a description in terms of localized Cu-site spins ${\bm{S}}_{i}$ and mobile oxygen holes. The second order terms contain direct O hopping terms and Cu-O Kondo coupling terms. We consider these terms in two complementary limits: (i) assuming an antiferromagnetically ordered background of Cu spins that renormalizes the oxygen bandstructure (main text), and (ii) disregarding both terms which yields a free oxygen dispersion (SI). Both calculations yields qualitatively identical results for quadrupolar response $\tilde{\Pi}^{\eta}$ (see Fig. 2(a)) and biquadratic exchange $K$. To fourth order appears the Heisenberg exchange interaction term and the exchange interaction term in Eq. (3) that includes the density of the intermediate O orbital. Upon integrating out the O holes this term yields the biquadratic exchange $K$ term of Eq. (2).

As the theory is quartic in O hole operators, we must first decouple the interaction terms. We perform the decoupling in the channel of the total and relative density of O atoms in a unit cell: $n^{p}_{i}=n^{p}_{i+\frac{\hat{x}}{2}}+n^{p}_{i+\frac{\hat{y}}{2}}$ and $\eta^{p}_{i}=n^{p}_{i+\frac{\hat{x}}{2}}-n^{p}_{i+\frac{\hat{y}}{2}}$. Introducing the vector $\nu_{\bm{k}}=(n^{p}_{\bm{k}},\eta^{p}_{{\bm{k}}})^{T}$, where $n_{i}^{p}=\frac{1}{N_{L}}\sum_{{\bm{k}}}n^{p}_{k}e^{i{\bm{k}}\cdot{\bm{R}}_{i}}$, the interaction terms read

[TABLE]

with interaction matrix

[TABLE]

Here, U_{\pm}=\frac{1}{8}\bigl{(}2V_{pp}f_{{\bm{k}}}\pm U_{pp}\bigr{)} and $f_{\bm{k}}=1+e^{-ik_{x}}+e^{ik_{y}}+e^{i(k_{y}-k_{x})}$. We decouple the interactions using a Hubbard-Stratonovich (HS) transformation, which yields an action that is purely quadratic in O operators, but contains the HS fields $\Phi_{\bm{k}}=(\psi_{\bm{k}},\phi_{\bm{k}})$:

[TABLE]

Here, $Q=(q,\alpha,\sigma)$ combines the Matsubara frequency-momentum index $q=(iq_{n},{\bm{q}})$ with index $\alpha=x,y$, which runs over the O orbital index $p_{x},p_{y}$, and $\sigma$, which denotes spin. The Green’s function $G^{-1}$ contains the Cu spin operators due to the coupling term $\propto J^{\prime}$ in Eq. (3), the O hopping terms and the HS variables: $G^{-1}=G_{0}^{-1}+G^{-1}_{\Phi}+G^{-1}_{J^{\prime}}$. It is diagonal in spin space, and is in orbital space $(\tau^{\alpha})$ given by

[TABLE]

with $h_{\bm{q}}=1-e^{iq_{x}}-e^{-iq_{y}}+e^{i(q_{x}-q_{y})}$. Here we suppress the second-order terms of the strong coupling expansion, which simply renormalize the dispersion. If one considers the motion of O holes in the background of AF ordered Cu spins, as we have done to calculate the results of Fig. 2(a), the dispersion entering $G_{0}^{-1}$ is modified accordingly Fischer-NJP-2014 ; SupplMat-CuprateNematicity . The other terms in the Green’s function read

[TABLE]

Here, $\tau^{0}\equiv\openone$ and we have defined the two functions

[TABLE]

with Cu spin bilinear $\mathcal{S}_{{\bm{p}},{\bm{q}}}={\bm{S}}_{\bm{p}}\cdot{\bm{S}}_{-{\bm{p}}-{\bm{q}}}$ and lattice functions $h_{n/\eta}({\bm{p}},{\bm{k}})=\frac{J^{\prime}}{2}(e^{ip_{x}}\pm e^{ip_{y}}+e^{-i(p_{x}+k_{x})}\pm e^{-i(p_{y}+k_{y})})$, where the upper (lower) sign relates to $h_{n}$ ($h_{\eta}$).

Integrating over the O degrees of freedom results in an action of the form

[TABLE]

We expand this expression to second order $S_{2}=\frac{1}{2}\text{Tr}[\{G_{0}(G^{-1}_{\Phi}+G^{-1}_{J^{\prime}})\}^{2}]$ in order to find

[TABLE]

which includes the biquadratic exchange $K$ term. Here, we have introduced the oxygen density response function

[TABLE]

The bare biquadratic exchange constant $K_{0}$ is given by the $zz$-component of this response function as

[TABLE]

Note that we write $\Pi^{zz}_{\bm{q}}\equiv\Pi^{\eta}_{\bm{q}}$ in the main text. It is straightforward to obtain the biquadratic exchange renormalized by O density fluctuations by performing the Gaussian integration over the HS fields in Eq. (22), which yields the renormalized response function

[TABLE]

where we have defined $(\tilde{U}_{q})=(U_{\bm{q}})_{\alpha\alpha^{\prime}}-\frac{1}{2}\Pi^{\alpha\alpha^{\prime}}_{q}$. Approximating the local response by the long-wavelength ${\bm{q}}=0$ component $\tilde{\Pi}^{zz}_{ii}\approx\tilde{\Pi}^{zz}_{{\bm{q}}=0}$ yields for renormalized biquadratic exchange constant $K$ as given in Eq. (4) of the main text.

Nematic susceptibility within soft-spin quantum field theory

In the main text, we analyze the nematic susceptibility $\chi_{\text{nem}}$ in the $t$-$J$-$K$-model using a soft-spin quantum field theory. This allows us to investigate the effect of quantum fluctutations on the nematic response. Our main results are shown in Fig. 3. After decoupling the biquadratic $K$ term using HS variable $\varphi_{r}$ the soft-spin action reads

[TABLE]

Here, ${\bm{M}}_{q}$ denotes an $N$-component Néel magnetization order parameter ($N=3$ in the physical system) and summation over repeated indices $\alpha,\beta$ is implied. The integrations are over $\int_{r}=\int_{0}^{1/T}d\tau\int d^{2}r$ and $\int_{q}=T\sum_{\omega_{n}}^{\Lambda_{\omega}}\int^{\Lambda}\frac{d^{2}{\bm{q}}}{(2\pi)^{2}}$ up to some dimensionless momentum and frequency cutoffs $\Lambda$ and $\gamma\Lambda_{\omega}$. The parameter $r_{0}$ controls the distance to the quantum critical point separating a Néel ordered regime from a quantum disordered paramagnetic regime, $u$ is an interaction constant and the coupling constant $g\propto K/J$ is proportional to the biquadratic exchange. We have added a source field $h_{\varphi}\equiv h_{r}$ that couples to homogeneous nematic order. We use a dynamic critical exponent of $z=2$ in the following, which describes damping due to particle-hole excitations in the presence of mobile holes.

The nematic susceptibility in Eq. (7) can be calculated from the partition function $Z=\int\mathcal{D}({\bm{M}}_{q},\varphi_{r})e^{-S}$ as

[TABLE]

where $L$ is the linear system size and the bare nematic susceptibility is given by

[TABLE]

with $\bar{\varphi}_{r}=\varphi_{r}+h_{\varphi}$. In the following, we consider homogeneous HS fields $\varphi_{r}$, $\bar{\varphi_{r}}$. To calculate the expectation value $\langle\bar{\varphi}^{2}_{r}\rangle$, we first decouple the quartic $u$-term in Eq. (27) using HS field $\psi$, then separate longitudinal and transverse components ${\bm{M}}_{r}=(\sqrt{N}M,\bm{\pi}_{r})$ and integrate over the transverse ones to arrive at the (dimensionless) action $s\equiv S/[L^{2}(\gamma T)^{-1}]$ given by

[TABLE]

Here, $r_{q}=r+{\bm{q}}^{2}+\gamma|\omega_{n}|$, $r=r_{0}+\psi$ and we have defined dimensionless interaction constants $\tilde{g}=g/\gamma$ and $\tilde{u}=u/\gamma$. Next, we expand the logarithm in small $\bar{\varphi}_{r}$ up to second order and obtain Eq. (29) by differentiation as

[TABLE]

We can exactly perform the summation over Matsubara frequencies, the momentum integration and then absorb the cutoff $\Lambda$ by expressing $\chi_{\text{nem,0}}$ in terms of the dimensionless variables $\tilde{T}=\gamma T/\Lambda^{2}$ and $\tilde{r}=r/\Lambda^{2}$. The lengthy expression is given in the Supplemental Material SupplMat-CuprateNematicity together with a three-dimensional plot as a function of $\tilde{T}$ and $\tilde{r}$. Cuts for different functional behaviors of the magnetic Néel correlation length on temperature $r(T)\equiv\xi^{-2}(T)$ are shown in Fig. 3(a).

We can derive the functional behavior of $r(T)$ within a large-$N$ approach, where we need to solve the following well-known self-consistency equation

[TABLE]

Solving this equation requires us to introduce a finite frequency cutoff $\Lambda_{\omega}$, but the qualitative behavior of $\chi_{\text{nem},0}$ and $\chi_{\text{nem}}$ does not depend on the cutoff choice as long as $\Lambda,\Lambda_{\omega}\gg r,\gamma T$. The results for $\chi_{\text{nem},0}$ and $\chi_{\text{nem}}$ shown in Fig. 3(b, c) are obtained from the large-$N$ solution of $r(T)$ for fixed parameters $\tilde{u},\Lambda,\Lambda_{\omega}$ and distance to the quantum-critical point $\delta r_{0}=r-r_{0,c}$.

Details on the classical Monte-Carlo simulations

The Monte Carlo simulations were carried out at $100$ equally spaced temperature points in the interval $0.001<T/J<2.971$. We applied a combination of single-move Metropolis Monte Carlo steps and non-local parallel-tempering-exchange steps between neighboring temperature configurations. The simulations shown in Fig. 2(b, c) of the main text were carried out for systems of $40\times 40$ spins and biquadratic exchange couplings $K/J=\{0.0,0,35,0.45\}$. We consider a ferromagnetic next-nearest-neighbor exchange coupling $J_{2}=-0.1J$ as well. Note that the ground state phase transition in the classical model between Néel and collinear order occurs at $J/2=J_{2}+K$. Following thermalization, the averages were computed for each temperature with at least $4.5\times 10^{6}$ Monte Carlo sweeps (MCS). The error bars were estimated by using the well-known Jackknife procedure.

Finally, we mention that we have performed Monte-Carlo simulations also for the purely bilinear spin Hamiltonian that is obtained from $H_{J-K}$ by using the well-known relations valid for spin-$1/2$ operators: \bigl{(}{\bm{S}}_{i}\cdot{\bm{S}}_{j}\bigr{)}^{2}=\frac{3}{16}-\frac{1}{2}{\bm{S}}_{i}\cdot{\bm{S}}_{j} and \bigl{(}{\bm{S}}_{i}\cdot{\bm{S}}_{j}\bigr{)}\bigl{(}{\bm{S}}_{i}\cdot{\bm{S}}_{k}\bigr{)}=\frac{1}{4}{\bm{S}}_{j}\cdot{\bm{S}}_{k}+\frac{i}{2}{\bm{S}}_{i}\cdot({\bm{S}}_{j}\times{\bm{S}}_{k}). These allow rewriting the biquadratic $K$ term as a sum of three bilinear spin exchange terms

[TABLE]

Here, $\delta(\delta^{\prime})$ runs over the (next-)nearest neighbors of the square lattice and $\delta^{\prime\prime}$ runs over the second-neighbors along the bonds. Importantly, classical Monte-Carlo simulation results for this Hamiltonian $\tilde{H}_{J-K}$ show the same enhancement of the nematic susceptibility $\chi_{\text{nem}}$ as a function of $K$ as results for the original Hamiltonian $H_{J-K}$ that includes the biquadratic exchange term.

Acknowledgments

We gratefully acknowledge helpful discussions with A. V. Chubukov, M.-H. Julien, B. Keimer, M. Le Tacon, and L. Taillefer.

Competing interests

The authors declare no competing interests.

Author contributions

P.P.O., B.J., R.M.F. and J.S. contributed extensively to the calculations, prepared the figures and wrote the paper.

Funding

P.P.O. acknowledges support from Iowa State University Startup Funds. J.S. acknowledges financial support by the Deutsche Forschungsgemeinschaft through Grant No. SCHM 1031/7-1. This work was carried out using the computational resource bwUniCluster funded by the Ministry of Science, Research and Arts and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHPC.

Data Availability

The data that support the findings of this study are available from the authors upon request.