Enhanced nematic fluctuations near the Mott insulating phase of high-$T_{c}$ cuprates

The complexity of the phase diagram of the cuprates goes well beyond its unique high-$T_{c}$ superconducting state, as it hosts a variety of different electronic phenomena, such as the pseudogap, nematic order, charge order, and strange metallic behavior. The parent compound, however, is well understood as a Mott insulator, displaying quenched charge degrees of freedom and low-energy antiferromagnetic excitations described by the Heisenberg exchange coupling $J$. Here we show that doping holes in the oxygen orbitals inevitably generates another spin interaction - a biquadratic coupling - that must be included in the celebrated $t-J$ model. While this additional interaction does not modify the linear spin wave spectrum, it promotes an enhanced nematic susceptibility that is peaked at a temperature scale determined by $J$. Our results explain several puzzling features of underdoped YBa$_{2}$Cu$_{3}$O$_{7}$, such as the proximity of nematic and antiferromagnetic order, the anisotropic magnetic incommensurability, and the in-plane resistivity anisotropy. Furthermore, it naturally accounts for the absence of nematicity in electron-doped cuprates, and supports the idea that the pseudogap temperature is related to strong local antiferromagnetism.

Introduction. Hole-doped cuprates are susceptible to a variety of different types of electronic order in the underdoped regime [1,2]. Examples include charge order [3][4][5][6], which becomes long-ranged in the presence of large magnetic fields [3,6], and nematic order [7][8][9][10][11], which seems to be stabilized only in compounds whose lattice structures explicitly break the tetragonal symmetry of the system [12,13]. While such a complexity might, at first glance, suggest that there is no universal explanation for the rich physics of hole-doped cuprates, two important guiding principles emerge from the analysis of their phase diagram (see Fig. 1): (i) novel phases other than superconductivity only appear below the ubiquitous pseudogap temperature T * ; (ii) T * (x) monotonically increases as one moves closer to the parent compound (x = 0), approaching values of several hundreds of Kelvin, comparable therefore to the antiferromagnetic (AFM) transition temperature of the parent compound.
This observation suggests a close connection between the pseudogap and AFM correlations [14][15][16]. Near x = 0, these AFM correlations arise from the well-understood Mott (or more precisely, charge-transfer [17]) insulating state of the parent compound. Such a connection between pseudogap and AFM [18] is consistent with the recent observation of a sudden change of the effective carrier number from x to 1 + x near the doping level where T * drops rapidly [19], or at a fixed doping level when T * is crossed [20]. Furthermore, it provides an interesting guiding principle to understand the tendencies towards the different types of electronic order below T * as also arising from strong AFM correlations.
In this context, several recent works have proposed that charge order arises from magnetic fluctuations [21][22][23]. In this paper, we demonstrate that nematic order -an electronic order that breaks the tetragonal symme-try of the CuO 2 plane [24] -is closely related to strong short-range AFM fluctuations as well. The key new ingredient of our approach to underdoped cuprates is the (1) The first two terms comprise the standard t − J model widely employed to describe lightly hole-doped Mott insulators [25], with t ij denoting the hole hopping parameters and J the AFM exchange coupling (weaker, longer range Heisenberg interactions are possibly included as well). The operator S i = 1 2 αβd † iα σ αβdiβ describes the Cu spin and n i = αd † iαd iα the corresponding charge. The strong local Coulomb interaction is incorporated in terms of the Cu-hole creation operator d † iα = (1 − n iᾱ ) d † iα , reflecting the fact that double occupancy of the sites is not allowed. The new term here is the biquadratic exchange term H K (the last term in Eq. (1)) which is governed by the coupling constant K > 0. As we derive below and show in Fig. 2(a), H K is naturally induced by non-critical charge fluctuations on the oxygen orbitals, which in turn affect the interaction between neighboring Cu-spins; K is also enhanced by the oxygenoxygen Coulomb repulsion. We note that the importance of O degrees of freedom for the nematic state of the cuprates has been pointed out in several interesting previous works [26][27][28][29][30].
Using Monte Carlo and analytical methods, we find as shown in Fig. 2(b,c) and Fig. 3, that the main consequence of H K is to enhance the static electronic nematic susceptibility χ nem near the AFM-Mott insulating state. However, χ nem is not found to diverge on its own. 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. It is natural that the CuO chains in YBa 2 Cu 3 O 7 (YBCO) play exactly this role of a symmetry breaking field. Such a scenario is supported by the experiments, which seem to constrain nematic order to compounds in which the tetragonal symmetry is explicitly broken [13]. The enhanced nematic susceptibility is however present in underdoped cuprates with fully intact tetragonal symmetry and can be determined via elastoresistance measurements [31]. Below we make a prediction for the outcome of such a measurement. We note that an enhanced nematic response was previously reported for the singleorbital Hubbard model in both weak-coupling [32] and strong-coupling [33] regimes. Our analysis implies that the inclusion of oxygen states strongly amplifies this behavior.
The nematic order resulting from H K explains the momentum anisotropy and the incommensurability of the spin-spin correlation function observed by neutron scattering in the paramagnetic phase of YBCO [8,34] as a function of the oxygen density np obtained within the three-band model at low temperature T = 10 −2 tpp and fixed n d = 1. Other parameters are t pd = 1, ∆ = 2.5, U dd = 9, Upp = 3, V pd = Vpp. Oxygen quadrupolar fluctuations increase with Vpp and smaller oxygen bandwidth tpp. The inset shows the resulting value of K/J ∝ (J 2 /J)(np/tpp) from which we conclude that enhancement of fluctuations by Vpp is crucial for a significant biquadratic exchange. We take a Hartree mean-field shift of ∆ eff = ∆+2np(Vpp −V pd −Upp/8) into account when calculating J , J. Panels (b-c) show the static nematic susceptibility χnem in Eq. (5) for the Ht−J−K model of Eq. (1) at half-filling as a function of temperature T obtained by Monte Carlo (MC) simulations of classical spins. A non-zero K enhances the response in the nematic B1g channel only. For consistency with the known spin-wave spectrum, we consider a small ferromagnetic next-nearest-neighbor exchange J2 = −0.1J.
well as the direction of the magnetic moment inside the AFM ordered state [35]. Scattering by these anisotropic fluctuations also address another hallmark of the nematic phase, namely, the in-plane resistivity anisotropy [7,36]. Furthermore, because H K is the result of charge fluctuations on the oxygen orbitals, it is only present in the holedoped side of the phase diagram, since electron-doping adds charge carriers directly to the Cu sites [37]. This property is again consistent with experiments, which, to the best of our knowledge, have not reported nematic correlations in electron-doped cuprates [12]. Microscopic model. Our starting point is the threeorbital Hubbard model of the CuO 2 -planes H = H 0 + H U + H V [38,41], which includes the 3d x 2 −y 2 Cu or-bital with creation operator d † iσ at position R i and spin σ as well as the 2p x and 2p y O orbitals (see Fig. 1(b)) with creation operators p † i+x 2 σ and p † i+ŷ 2 σ . The noninteracting part H 0 includes hopping between oxygen orbitals (with amplitude t pp ) and between copper and oxygen orbitals (with amplitude t pd ). The corresponding sign factors of the hopping elements follow from the phases of the orbitals as shown in Fig. 1(b), see Ref. [38].
In addition H 0 contains on-site terms where the energy difference between Cu and O orbitals is given as ∆ = ε p − ε d . We further include on-site interactions The energy scales are the local repulsion U dd at the Cusites and the charge-transfer energy ∆ (with U dd much larger than ∆), suggesting a strong coupling expansion in small t ij ∆, U dd −∆. This yields a description in terms of localized Cu spins S i coupled to mobile O holes. An expansion up to fourth order in the hopping term t pd was performed in Ref. [39]. There appear (for details see the Supplementary Material (SM) [41]) Kondo-like exchange couplings ∝ S i · s jk between Cu and O spin-densities, s jk = 1 2 αβ p † jα σ αβ p kβ , a Heisenberg exchange term J i,j S i · S j , and, most notably for our considerations, a spin exchange term that depends on the occupation of the intermediate oxygen orbital between copper sites: where δ = {±x, ±ŷ} and the coupling constant is 2∆ 3−n (U dd −∆) n and in the large-U dd limit J /J → 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. (1). This follows from decomposing the oxygen densities as n p i+x 2 = n p i + η i and n p i+ŷ 2 = n p i − η i , where η i is the quadrupolar (nematic) component of the oxygen charge density [27]. The combination of O on-site and nearest-neighbor Coulomb interactions leads to a term ( 2 ) and f k = 1+e −ikx + e iky + e i(ky−kx) . Note that the observed charge density order of the cuprates has been interpreted as being dominated by a quadrupolar form factor of this type [3,4], suggesting that it may be related to a condensation of the Fourier transform η q of η i at a finite momentum [43]. In what follows we do not focus on this closely related ordered state, but integrate out the quadrupolar oxygen charge fluctuations. The details of this analysis are summarized in the SM [41] and yield the result for the biquadratic exchange interaction in Eq. (1) Here, Π η k = − q,ω Tr[G p q,ω (τ z σ 0 )G q+k,ω (τ z σ 0 )] is the bare oxygen charge susceptibility in the quadrupolar (i.e. nematic) channel. The Pauli matrices σ i and τ j act in spin and p x,y −orbital space, respectively. Explicit expressions for the oxygen Green's functions G p q,ω are given in [41] and yield Π η Here, q is the (effective) oxygen dispersion and µ the chemical potential. Eq. (3) makes it clear that the biquadratic exchange in Eq. (1) is a natural consequence of quadrupolar oxygen charge fluctuations.
The oxygen quadrupolar susceptibility Π η k=0 (and thus K) is positive and determined by the occupation number difference between the upper and lower oxygen bands. In the relevant regime of small hole fillings n p 1, the response approaches a constant Π η k=0 ∝ n p at low T , peaks around T ≈ |µ| and vanishes as 1/T at large T . In Fig. 2(a) we present results for the renormalized quadrupolar responseΠ η k = 1 2 (Π η k −1 − U k −1 and for K/J within the three-band model, keeping n d = 1 and neglecting the interaction of the oxygen electrons with the magnetic background. We clearly observe that a large nearest-neighbor repulsion V pp and a small bandwidth t pp enhance K. Note that 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 to the physics of underdoped cuprates 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 Q = (π, π), but the striped configuration with Q = (π, 0) or (0, π) that is observed in many ironbased superconductors. 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 demonstrate this 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 stripedmagnetic 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.
A natural disturbance of the tetragonal symmetry is, of course, provided by the CuO chains or double chains in YBa 2 Cu 3 O 7−δ and YBa 2 Cu 4 O 8 , respectively, or by external strain [31].
The behavior described above can be quantified in terms of the composite spin variable: which changes sign under a rotation by π/2. Note that the square of this term, which appears in Eq.(1), is invariant under this transformation, and therefore is fully consistent with the four-fold symmetry of the CuO 2 -plane. While ϕ i = 0 for realistic values of K (and in the absence of external strain), the static nematic susceptibility is a natural measure for the increased relevance of local stripe magnetic configurations. Here, T τ denotes imaginary time ordering. In Fig. 2(b-c), we show classical Monte Carlo results for χ nem for a collection of classical Heisenberg spins that interact according to the H J−K model with an additional small second-neighbor exchange J 2 = −0.1J to obtain consistency with the known spin-wave spectrum. One clearly sees that the biquadratic term K enhances the nematic response in the B 1g (x 2 − y 2 ) channel, corresponding to an inequivalence between the x and y axes. Most importantly, the nematic susceptibility χ nem (T ) is non-monotonic, peaking at a temperature governed by the effective exchange interaction of the Cu-spins, T nem ∼ J, which is independent on K.
The Monte Carlo results also display that χ nem (T → 0) → 0, which is a consequence of the classical nature of the spins in the simulations. As shown in Fig. 3, quantum fluctuations crucially modify this behavior and lead to χ nem (T → 0) > 0. We include the effect of quantum fluctuations by analytically deriving the nematic response within the limit of a large number of (quantum) spin components N . The large-N analysis uses a soft-spin version of the spin degrees of freedom in Eq. (1). 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: where g ∝ K/J > 0, u > g, r ≡ 1/T 0 dτ d 2 r and r = (τ, r) combines imaginary time τ and position r = (x, y). In addition, ϕ r is the nematic order parameter of Eq. (4), h r is an external strain field, and n r is the N = 3 component staggered Néel order parameter, as used in the non-linear sigma model description of the cuprates in Refs. [42,44]. The quantum dynamics of the Néel order parameter is governed by S dyn = q f (ω n ) n q · n −q , where f (ω n ) ∝ ω 2 n at half filling, while f (ω n ) = γ |ω n | was proposed to describe particle-hole excitations. Here, q = (ω n , q) combines Matsubara frequency ω n and momentum q (measured relative to the AFM ordering vector Q = (π, π)) and q ≡ T n d 2 q (2π) 2 . The parameter r 0 measures the distance to the AFM quantum critical point, located at r 0,c . For δr 0 ≡ r 0 − r 0,c < 0, the system has long-range AFM order at T = 0, whereas for δr 0 > 0 it is in the paramagnetic phase (see Fig. 3(c-d)).
In order to make analytic progress we consider the limit where n r is an N -component vector and take the limit of large N . This approach led to important insights in both the description of antiferromagnetic correlations of the cuprate parent compounds [44] and of nematic fluctuations of iron-based superconductors [45]. Computing the nematic susceptibility yields (see [41]): where the bare nematic susceptibility is given by χ Here, ξ is the magnetic correlation length for Néel order, which is also determined self-consistently within large-N for a given distance to the AFM quantum critical point, δr 0 . Despite the similarity between Eq. (7) and the expression for the nematic susceptibility of iron-based superconductors [45], there are very important differences between the two systems. Because the iron pnictides order magnetically in a striped configuration, χ (0) nem diverges when ξ → ∞, which guarantees that a nematic transition takes place already in the paramagnetic state for any g > 0. However, because the cuprates order in a Néel configuration, χ (0) nem remains finite even when ξ → ∞. Although long-range nematic order is not present, nematic fluctuations can be significantly enhanced if the biquadratic exchange K ∝ g is sufficiently large.
In Fig. 3 we show the nematic susceptibility obtained within the large-N approach [41]. Like in the Monte Carlo results (see Fig. 2), we observe a broad maximum at finite temperatures around T ≈ J (the lattice cutoff Λ plays the role of J in the continuum model). This pronounced peak of χ nem is determined by χ (0) nem , which in turn is governed by the magnetic correlation length ξ that is set by T /J. The effect of g, and thus of the biquadratic exchange K, is to determine the amplitude of the peak. At low temperatures, we find that in contrast to the results of the classical Monte Carlo simulation, quantum fluctuations render χ nem finite, which is shown in Fig. 3(b). We further observe that the nematic response increases for an increasing magnetic correlation length, i.e. Néel fluctuations enhance the nematic susceptibility, since χ nem (δr 0 , T ) is an increasing function for decreasing δr 0 . Note that χ nem is a non-universal quantity whose precise shape depends on microscopic details (such as the lattice constant) and will be different for different materials.
Consequences of long-range nematic order. An immediate consequence of this enhanced nematic susceptibility is that a small tetragonal-symmetry breaking field h, caused either by external strain or by the presence of CuO chains, can induce a sizable nematic order parameter ϕ ≈ χ nem h. From the action in Eq. (6), we can readily obtain the dynamic spin susceptibility in the presence of this induced nematic order Therefore, as shown in Fig. 4, non-zero ϕ modifies the In the absence of nematic order (ϕ = 0) the scattering peak is isotropic around the Néel ordering vector Q = (π, π), but nonzero 0 < ϕ < 1 leads to an elliptic deformation of the peak as observed experimentally [34]. For larger values of ϕ > 1 the peak splits and two incommensurate scattering peaks emerge at (π ± δ, π). Note that the inequivalence of x and y direction appears in response to an external (or intrinsic) strain field that explicitly breaks C4 symmetry.
spin-spin structure factor near the Néel ordering vector Q from a circular shape, which preserves tetragonal symmetry, to an elliptical shape, which breaks tetragonal symmetry. In addition, as ϕ increases, it can naturally shift the maximum of χ AFM (Q + q, ω) from the commensurate q = 0 value to an incommensurate q IC = 0 value, with q IC parallel to either the x axis (if ϕ > 0) or to the y axis (if ϕ < 0). Interestingly, both features -the elliptical structure factor and an anisotropic incommensurate peak -have been observed in YBCO by neutron scattering experiments [8,34]. In agreement with our model, these effects onset at a temperature scale comparable to the Néel transition temperature in underdoped compositions. Furthermore, Ref. [36] has shown that scattering of itinerant carriers by these anisotropic spin fluctuations cause an anisotropy in the resistivity anisotropy that is consistent with the transport data in YBCO in the same region of the phase diagram [7]. Note that a somewhat related mechanism for the incommensurate spin order, based on the t − J model, was reported in Refs. [46,47,49]. Previous works have also focused on nematicity arising from a pre-existing incommensurability [48], whereas in our scenario incommensurate magnetic order is a consequence of nematic order. While these effects onset in the paramagnetic phase, the presence of nematic order is also manifested in the Néel ordered state by the direction of the Cu moments, which align parallel to either the x axis or to the y axis [35]. Within our model, this observation can be rationalized by invoking the spin-orbit coupling in the oxygen sites, which convert an imbalance in the charge of the 2p x and 2p y orbitals into a preferred direction for the Cu moment.
At first sight, one might anticipate that K would af-fect the spin-wave dispersion, which, experimentally, is well described in terms of J and a small ferromagnetic nearest-neighbor exchange J 2 only [50]. As we show in the SM [41], however, the biquadratic exchange of Eq. (1) does not modify the linearized 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: Thus, while it seems impossible to gain insight about the biquadratic exchange from measurements of the spin wave dispersion, it is also true that such measurements are fully consistent with the new Hamiltonian Eq. (1).
Comparison with experiments and conclusions. In summary, we showed that charge fluctuations associated with hole-doping in the oxygen orbitals generate a biquadratic exchange coupling between the Cu spins, extending the celebrated t − J model employed to describe lightly-doped cuprates. The main effect of this biquadratic term is to enhance B 1g nematic fluctuations, which however is not translated into a diverging nematic susceptibility. Most importantly, the temperature in which the nematic susceptibility peaks is determined not by the biquadratic coupling K, but by the standard nearest-neighbor exchange coupling J. The amplitude of the peak, however, is controlled by K, which increases for larger values of V pp .
These results are consistent with the current experimental scenario: first, nematic order seems to be observed only in hole-doped compounds whose lattice structures explictly break tetragonal symmetry [12,13]. Within our model, such a small symmetry-breaking field takes advantage of the large (but non-diverging) nematic susceptibility and induces nematic order. Second, the observed temperature scale for the onset of the induced nematic order is comparable to the Néel transition temperature, which in turn is set by J [8,34]. Third, our model naturally addresses the main manifestations of the induced nematic order observed experimentally in the spin spectrum -the elliptical spin-spin structure factor, the anisotropic incommensurate AFM peak, and the direction of the Cu moments inside the AFM state [8,34,35].
To further verify the validity of our scenario, it would be desirable to directly measure the nematic susceptibility in cuprates. Our main result is the non-monotonic temperature dependence of χ nem , with a peak at T ∼ J.
In analogy to what has been done for the iron-pnictides (see Ref. [51]), χ nem is closely related to several observables, such as the elastoresistance, the shear modulus, or electronic Raman scattering. According to our model, performing these measurements in underdoped tetragonal systems, such as HgBa 2 CuO 4 , should reveal an enhanced nematic response without long-range order. Similarly, we predict this behavior to be much weaker in tetragonal electron-doped systems, such as Nd 2 CuO 4 , where oxygen charge fluctuations are expected to be less relevant, since electron-doping involves Cu orbitals.
Finally, an interesting open question is the effect of such enhanced nematic fluctuations on the superconducting state. In particular, the close connection between nematic and Néel spin fluctuations offers an interesting scenario to investigate the formation of Cooper pairs. Indeed, recent theoretical studies have revealed that nematic fluctuations, in conjuction with another pairing mechanism, may provide a sizable boost to the superconducting transition temperature [52].
We gratefully acknowledge helpful discussions with A. Chubukov Here, d † iσ creates a Cu (3d x 2 −y 2 ) hole with spin σ at Bravais lattice site R i . The operators p † i+x 2 σ and p † i+ŷ 2 σ create O (2p x ) and (2p y ) holes in the same unit cell R i , respectively. The vacuum is defined as filled Cu + (d 10 ) and O 2− (p 6 ) states. We define the total number of Cu holes with spin σ in unit cell R i as n d iσ = d † iσ d iσ and the Cu hole density as n d i = σ n d iσ . The corresponding operators for oxygen holes read n p iσ = u=x 2 ,ŷ 2 p † i+uσ p i+uσ and n p i = σ n p iσ . The on-site energies of d(p) orbitals are denoted d ( p ) with ∆ = p − d > 0. The phase factors in the hopping terms arise from the overlap of orbital wavefunctions (see Fig. 1(b) of the main text) and are given by (−1) u = +1 for u = −x 2 ,ŷ 2 , (−1) u = −1 for u =x 2 , −ŷ 2 and (−1) u = +1 for u = ± 1 2 (x +ŷ), (−1) u = −1 for u = ± 1 2 (x −ŷ). The (unrestricted) sum u runs over the four vectors connecting a central Cu site to its four neighboring O sites u ∈ {±x 2 , ±ŷ 2 }, and the sum u runs over the four vectors connecting an O p x orbital to its four p y neighbors u ∈ {± 1 2 (x ±ŷ)}. We consider on-site interactions U dd and U pp on both Cu and O sites as well as nearest-neighbor interactions V pd between Cu and O and V pp between oxygens.

S1.A. Strong-coupling expansion
The hierarchy of energy scales suggests a strong-coupling expansion in small t pd U dd − ∆, ∆, which yields a description in terms of localized Cu spins coupled to mobile O holes. Note that in order to derive the biquadratic exchange interaction term ∝ K in Eq. (1) of the main text, one can focus on the case of singly occupied Cu sites, i.e., all holes reside on the oxygen sites. In the strong coupling expansion we follow Ref. 2

(see also Ref. 3) that contains an expansion up to fourth order in t 4
pd /[∆ n (U dd − ∆) m ] with n + m = 3 within this subspace. At second order one finds a term that renormalizes t pp and a Cu-O Kondo like exchange coupling term H (2) dp = 2t 2 with (non-)local O spin operators s ij = 1 2 τ,τ p † iτ σ τ τ p jτ with σ = (σ x , σ y , σ z ) being a vector of Pauli matrices and Cu spin operators At fourth order, there appear further Cu-O Kondo-like exchange terms, the well-known Heisenberg Cu-Cu spin exchange term and a term that renormalizes hopping and interactions among oxygen sites. Most importantly for our analysis, however, there also appears a Cu spin exchange term that depends on the hole occupation number of the intermediate O orbital and δ = {±x, ±ŷ}. After integration of oxygen density fluctuation, this term will give rise to the biquadratic spin exchange as we show below. Note that in the limit of large-U dd , one finds lim U dd →∞ J /J = 1/2, so this term is of the same order as the Heisenberg spin-exchange.
There also appear non-local manifestations of this term with n p i being replaced by n p ij = σ p † iσ p jσ and term that describes coupling of the (non-)local spin density on the intermediate oxygen site to the cross product of Cu spins H (4) ddp ∝ −2i i,j ,u1,u2 s i+u1,j+u2 · (S i × S j ) , which are, however, not the focus of our analysis.

S1.B. Derivation of biquadratic spin exchange K0 from microscopic Hamiltonian
To consider the effect of charge fluctuations on the oxygen sites, we rewrite the oxygen interaction part of the Hamiltonian H Upp + H Vpp in terms of total n p i and relative oxygen density η i within unit cell R i : This allows to write the oxygen-density-dependent interaction between neighboring Cu spins in Eq. (S.6) as The on-site (U pp ) and nearest-neighbor (V pp , V pd ) oxygen interaction terms in Eq. (S.2) and (S.3) take the form Here, δ ∈ {0,x, −ŷ,x −ŷ} denotes Bravais lattice vectors pointing to unit cells containing the four nearest-neighbor p y oxygen orbitals of a given oxygen p x orbital.
As required, Eq. (S.10) is invariant under the spatial symmetries of the system. In particular, it is invariant under fourfold C 4 rotation C 4 (x i , y i ) = (−y i , x i ), C 4 (k x , k y ) = (−k y , k x ). This follows from the transformation laws of the orbitals p i+x 2 C4 − − → p C4(i)+ŷ 2 and p i+ŷ 2

C4
− − → p C4(i)−x 2 . The transformation laws for total and relative densities follow as η k Since the operators η k and n p k transform into each other under symmetry transformations (such as C 4 ), we need to treat them on equal footing when decoupling the interaction terms using a Hubbard-Stratonovich (HS) transformation. We introduce the vector v k = (n p k , η k ) T and write H Upp (S.11) The HS transformation introduces the fields Φ k = (ψ k , φ k ) and yields the action where k = (ik n , k) combines Matsubara frequency ik n = 2πnT with temperature T and momentum k. We note that the fields Φ k transform identical to v k in order that the action remains invariant under all symmetry transformations.
We have arrived at an action that is quadratic in oxygen hole operators. In a next step, we will perform the exact functional integration over these degrees of freedom. We focus on those terms in the action that are relevant for the derivation of the biquadratic exchange interaction S = S 0 + S Upp+Vpp + S J , which read explicitly where u, u ∈ {x, y} label (p x , p y ) orbitals and σ, σ ∈ {↑, ↓} denote the spin direction. The inverse Green's function G −1 (S i ) in Eq. (S.13) contains the Cu spin operators S i and is a sum of terms Since all Green's functions are diagonal in spin space, G −1 ∝ σ 0 with σ 0 = diag(1, 1), we suppress the spin indices in the following and find Here, τ α are Pauli matrices in orbital (p x , p y ) space, τ 0 = diag(1, 1) and we have defined the lattice function that describes oxygen hopping. We have also introduced the functions Functional integration over p † quσ and p quσ yields the action where the ellipsis stands for the zeroth, first and higher order terms. Focusing on the quadratic term, we write it as where we have introduced the response functions Performing the summation over Matsubara frequencies and setting the external frequency iq n to zero, they read and Π 0z q = −Π z0 q . Here, n F (x) = [exp(x/T ) + 1] −1 is the Fermi function and we have defined∆ = ∆ − µ, q = t pp |h q | and g ±,q,k = 1 2 t 2 pp h * q h q+k ± c.c . We note that the response functions Π αβ are only C 2 symmetric, but their sum as it appears in the action is always fully C 4 symmetric, which we have verified explicitly. Importantly, in the long wavelength limit, one finds that both lim |k|→0 Π 00 > 0 and lim |k|→0 Π zz > 0. This determines the sign of the biquadratic exchange coupling K > 0 as given in Eq. (1) of the main text.
Focusing on the bare biquadratic term arising from the product of the operators S p,q in Eq. (S.23), it is instructive to write it in real space as The bare biquadratic exchange coupling K 0 follows as We note again that the action S S 2 2 is fully C 4 invariant and the interaction terms that involve spin operators in neighboring unit cells i −x and i −ŷ arise from the the off-diagonal components Π z0 ij . In Fig. S.1 we show the response functions Π αβ both in real and momentum space for a realistic choice of parameters t pd = 1, t pp = 0.2, ∆ = 2.5, n p = 0.05, U dd = 9, U pp = 3, V pp = 2.0, V pd = 1.0 and temperature T = 0.02.

S1.C. Analysis of bare biquadratic exchange coupling constant K0
To gain some analytic understanding, we approximate the bare biquadratic coupling constant by are only C2-symmetric, the resulting expression in the action is fully C4-symmetric due to multiplication with the lattice functions hm, hη. We observe that Π zz q peaks at a non-zero wavevector showing that the maximal nematic response occurs at a finite q.
Here, ξ ± = ± q − µ describes the two oxygen bands. Neglecting the interaction with the Néel magnetic background of the localized Cu spins, we obtain the oxygen bandstructure from the lattice function h q = 1 − e iqx − e −iqy + e i(qx−qy) : It is useful to introduce the density of states [complete] elliptic integral of the first kind. We note that the density of states logarithmically diverges as → 0 due to a van-Hove singularity. We can now analyze the low and high-temperature behavior of the nematic response , at low T 2t pp S1.D. Renormalization of biquadratic exchange coupling K by quadrupolar oxygen density fluctuations We now show that quadrupolar oxygen density fluctuations further enhance the biquadratic spin exchange from its bare value K 0 to a renormalized value K > K 0 . These oxygen density fluctuations become stronger for increasing oxygen-oxygen repulsion V pp . They are described by the bosonic Hubbard-Stratonovich fields Φ q = (ψ q , φ q ) defined above Eq. (S.12). In the parameter regime we consider these fluctuations are non-critical and thus remain massive. Nematic order does not develop spontaneously, but occurs only in the presence of a conjugate symmetry-breaking field such as strain or as provided by the CuO chains in YBCO. After integration over the fermionic fields p † quσ and p quσ the action reads (see Eq. (S.22)) where α, β ∈ {0, z} and we identify h 0 (k, q) ≡ h m (k, q) and h z (k, q) ≡ h η (k, q). Performing the Gaussian integration over Φ † q yields where we have defined with U q given in Eq. (S.11). We note that the action in Eq. (S.36) is fully C 4 symmetric after summation over α, β, γ, γ , which we have explicitly verified. We can readily extract the renormalized response functionsΠ αβ q as The renormalized biquadratic exchange interaction is therefore determined by the zz componentΠ zz q = Π zz q + To gain more insight, we approximate the local responseΠ zz ii , which determines the biquadratic coupling K = J 2 2Π zz ii (see Eq. (S.29)), by the q = 0 componentΠ zz q=0 . Using that (Ũ −1 , the biquadratic exchange K, renormalized by quadrupolar oxygen density fluctuations, is given by We show K/J as a function of hole doping n p in Fig. 2(a) of the main text for realistic parameters of the cuprates. The main insight from this result is that while on-site oxygen interactions U pp tend to reduce the biquadratic exchange, repulsive oxygen-oxygen interactions V pp enhance the biquadratic spin coupling K > K 0 . For realistic parameters of the cuprates it holds that V pp U pp /8 and the enhancement due to V pp is the dominant effect.

S2. Nematic susceptibility within a soft-spin description of the half-filled t − J − K-model
In order to obtain an analytic understanding of the nematic response in the presence of a biquadratic exchange term ∝ K close to a Néel ordered state, we investigate a soft-spin version of the two-dimensional t-J-K-model at half-filling. We have also analyzed the nematic susceptibility in spatial dimensions 2 < d ≤ 3 and found that the results for the nematic response from d = 2 remain qualitatively unchanged. After decoupling the biquadratic term at the expense of introducing the Hubbard-Stratonovich field ϕ r , the action reads Here, M q = (M 1 q , M 2 q , . . . , M N q ) with q = (iω n , q) combining Matsubara frequency iω n = 2πinT and momentum q = (q x , q y ) denotes the (dimensionless) N -component staggered Néel magnetization. The integrations run over (2π) 2 with dimensionless momentum and frequency cutoffs Λ and γΛ ω and r = β 0 dτ d 2 r. The coupling constant g ∝ K/J is proportional to the ratio of biquadratic exchange K to nearest-neighbor Heisenberg exchange J in the spin model. The (bare) mass parameter r 0 controls the distance to the quantum critical point between Néel ordered and a paramagnetic T = 0 phases, andũ = u/γ is a dimensionless interaction constant. We have rescaled the interaction term u/N to obtain a well-defined large-N limit. In the following we set the dynamic critical exponent to z = 2, which describes damping due to particle-hole excitations in the presence of doped holes. We have added a source field h ϕ (denoted h r ≡ h ϕ in the main text) that couples to homogeneous nematic order with inverse temperature β = 1/T and bare nematic susceptibility (S.42) To obtain this expression, we have shifted the fieldφ r = ϕ r + h ϕ in Eq. (S.40) before taking the derivatives with respect to h ϕ and assumed the absence of an external field h ϕ = 0 so that φ r 2 = 0. We focus on static and homogeneous Hubbard-Stratonovich fields ϕ r andφ r , i.e., both fields are independent of r. To analytically calculate the expectation value φ 2 r , we decouple the quartic term ∝ u in Eq. (S.40) by defining the (dimensionless) density ρ ψ = 1 N M x · M x and introducing a factor of unity as 1 = D(ρ ψ , ψ)e − 1 γ r (M 2 x −N ρ ψ )ψ [4], to arrive at Here, the dimensionless field ψ describes the renormalization of the mass from r 0 → r ≡ r 0 +ψ. Separating longitudinal and transverse components M r = ( √ N M, π r ), where we restrict to homogeneous magnetic order M , and integrating over the (N − 1) transverse components, we arrive at the action density s: Next, we expand the logarithm in smallφ to find where q = |q|(cos θ, sin θ) and r q = r + q 2 + γ|ω n | with r = r 0 + ψ. The generating functional ofφ r then reads 46) and the bare nematic susceptibility χ nem,0 in Eq. (S.42) follows to .
(S. 47) In the final step we have gone from summation over Matsubara frequencies to integration along the real frequency axis, and performed the angular integration over θ. Note that while ω n has units of energy, the integration variable ω is dimensionless by expressing energies in units of γ −1 . While we must keep both momentum and frequency cutoff Λ and γΛ ω finite when we solve for r(r 0 , T ) in the following section (using the large-N approximation), we may take the limit γΛ ω → ∞ in Eq. (S.47). This allows us to exactly perform the momentum and frequency integrations and completely absorb the cutoff Λ by expressing χ nem,0 in terms of the (dimensionless) variablesr = r/Λ 2 andT = γT /Λ 2 as χ nem,0 = N Λ 4 64π 2 γ 4πT 2r + 1 r + 1 + 2r logr r + 1 + 2ψ 1 +r + 1 2πT In Fig. S.2, we show χ nem,0 /N , where N = N Λ 4 /(64π 2 γ) as a function ofr = r/Λ 2 andT = T /Λ 2 . Note that N ≡ χ nem,0 (r = 0, T = 0) is the value of the susceptibility at the quantum critical point. Along a path of constant temperature, the nematic susceptibility increases as r decreases, which is also shown in the upper right panel of Fig. S.2. Decreasing r ∝ ξ −2 implies an increasing magnetic Néel correlation length as one approaches the T = 0 quantum critical point or the renormalized classical regime with exponentially large magnetic correlation length at T > 0. The nematic susceptibility thus increases as a result of larger magnetic Néel fluctuations. Our analysis also reveals that while for classical spins χ nem,0 vanishes as T → 0, quantum fluctuations render the zero temperature limit of χ nem,0 finite. To plot χ nem,0 along a path of constant r 0 , which controls the distance to the quantum critical point beyond which Néel order disappears, one needs to solve for the renormalized mass parameter r(T, r 0 ). At finite temperatures above the quantum critical point, one finds r(T, r 0,c ) = aγT ∝ T with a non-universal proportionality constant a that depends on microscopic details of the system. As shown in the lower right panel of Fig. S.2, the shape of χ nem,0 crucially depends on the value of the slope parameter a, which controls the relative importance of quantum and thermal fluctuations. For small values of a < π, χ nem,0 develops a pronounced finite-temperature peak, whose amplitude increases with decreasing a. This behavior of χ nem,0 closely resembles the behavior found within the classical Monte-Carlo simulations (see Fig. 2 of the main text), and show that thermal fluctuations are dominant for a < π. In contrast, for larger values a > π, χ nem,0 peaks at T = 0 and is a monotonically decreasing function for finite T , showing the dominance of quantum fluctuations in this case.

S2.B. Large-N analysis of the nematic susceptibility
In order to determine the effective mass parameter r(T, δr 0 ) as a function of temperature T and distance to the quantum critical point δr 0 = r 0 − r 0,c , we consider the limit of large-N , where the partition function is governed by the saddle point of the action in Eq. (S.40). Finding the saddle-point of the action in Eq. (S.44) in the absence of nematic order, leads to the well-known large-N self-consistency equations [4,5]: rM = 0 with r = r 0 + ψ; ρ ψ = γψ/u and The equation requires a finite frequency cutoff Λ ω . For the results of χ nem,0 in Fig. 2 of the main text, we have therefore numerically solved these equations for finite momentum Λ and (dimensionless) frequency cutoffs γΛ ω , which yields r(r 0 , T ) shown in Fig. S.3. The qualitative behavior of χ nem,0 as discussed in the previous section does not depend on the exact value of γΛ ω and Λ as long as both cutoffs are much larger than r, γT Λ, γΛ. This parameter controls the distance to criticality. The resulting nematic susceptibilities are shown in Fig. 3 of the main text. Panels (a) and (b) correspond to same panels in Fig. 3 of the main text. The left panel (a) shows r as a function of temperature at fixed distance to criticality δr0, as depicted in the schematic diagram with vertical paths corresponding to different values of δr0. For δr0 > 0 (red), r approaches a finite value as T → 0 corresponding to a finite magnetic correlation length in the quantum disordered phase. In contrast, for δr0 ≤ 0 (yellow, green), r → 0 as T → 0. Approaching the quantum critical point (yellow), one finds r = aγT with non-universal slope a that depends among others on the size of the interactions u/γ. Larger u/γ yields a larger a. Approaching the magnetically ordered phase (green), we observe a change in functional behavior of r(T ) as we cross from the quantum critical to the renormalized classical regime. In the quantum critical regime at higher T , one finds r(T ) ∝ T . In contrast, in the renormalized classical regime at lower T , it holds r(T ) ∝ T exp(−∆(δr0)/T ), where ∆ depends on the spin stiffness in the ordered phase (rigidity towards magnetic fluctuations) [5]. The right panel (b) corresponds to changing δr0 at fixed temperature T . In an experiment, this corresponds, for example, to tuning the chemical composition or pressure. Importantly, we observe that r ∝ ξ −2 is a monotonously decreasing function for decreasing δr0, i.e., approaching the quantum critical point or the renormalized classical region. While in two dimensions the Hohenberg-Mermin-Wagner theorem ensures a finite correlation length r > 0 at any finite temperature T , one finds that r(T ) becomes exponentially small, because ξ becomes exponentially large, in the renormalized classical regime to the left of the dashed line. From Fig. 2 of the main text, we conclude that the nematic susceptibility χnem increases if r decreases or ξ increases, which shows that larger Néel fluctuations enhance the nematic response.
S3. Spin-wave treatment of t − J − K model at half-filling At half-filling the Hamiltonian in Eq. (1) of the main text describes localized Cu spins interacting via nearestneighbor Heisenberg exchange interaction J and a biquadratic exchange interaction K. Additional weaker next-nearest neighboring and ring-exchange terms may be added to obtain agreement with the experimental spin-wave spectra [6]. Since these terms do not change our conclusions, we do not explicitly consider them below. Note that we include a realistic ferromagnetic exchange coupling J 2 = −0.1J in our Monte-Carlo simulations (see Fig. 2 of the main text).
We now show that the biquadratic spin exchange term ∝ K does not modify the spin-wave spectrum. Adding such a term in the Hamiltonian is thus fully consistent with previous experimental results of the spin-wave spectrum. We derive our results starting from the J-K model spin Hamiltonian in Eq. (1) of the main text (S.50) Here, {δ ν } = {±x, ±ŷ} connect nearest-neighbors Cu sites. Let us calculate the spin-wave spectrum around the Néel state. As this corresponds to a large-S limit, we have rescaled the biquadratic term. We follow the standard procedure of Holstein-Primakoff spin-wave calculations [7] and begin with defining local triads n 1,Ri = cos(Q · R i ), 0, − sin(Q · R i ) , n 2 = 0, 1, 0 and n 3.Ri = sin(Q · R i ), 0, cos(Q · R i ) with Néel ordering wavevector Q = (π, π). Expressing the spins in terms of this local coordinate system S i = αS α i n α,Ri , the Hamiltonian takes the form (suppressing the tilde) with n α β,δν = δ αβ (−δ xα + δ yα − δ zα ). A transformation to momentum space via S i = 1 p,q,k α,β S α p+k S α −p S β q−k S β −q n α α,δν n β β,δν cos p x − cos p y )(cos q x − cos q y ) , (S.52) where we have defined the lattice function f p = 1 4 δν e −ipδν = 1 2 (cos p x + cos p y ). To obtain the spin-wave spectrum, we now express spin operators in terms of Holstein-Primakoff bosons S x q = S 2 b † −q + b q , S y q = i S 2 b † −q − b q and S z q = √ N L Sδ q,0 − 1 The classical ground state energy follows as the O(S 2 )-term to H (S 2 ) /(N L S 2 ) = −2J. Note that the energy of the biquadratic term vanishes in the Néel state.
The next lowest order in S is quadratic in the bosons and yields upon diagonalization the spin-wave spectrum. Keeping the quadratic terms, yields Most importantly, according to Eq. (S.54), the biquadratic term does not contribute to the spin-wave spectrum at order O(S). Intuitively, the vanishing of the biquadratic contribution to the spin-wave spectrum follows from the observation that inserting the classical spin state into one of the factors S i (S i+x + S i−x − S i+ŷ − S i−ŷ ) gives zero, because the term in the brackets vanishes in the Néel state. More explicit, this result can be seen already from Eq. (S.52): to obtain a term of O(S) the term in the bracket has to be of order S 3 in order to combine with the prefactor K/S 2 to a term of O(S). There are two possibilities, which both vanish: (i) zzzz terms, i.e., selecting 3 Kronecker symbols and selecting one term of the form k b † k−q b k (from S z q ). The Kronecker symbols enforce p = q = k = 0 and therefore cos p x − cos p y = 0 and cos q x − cos q y = 0. The other possibility are (ii) terms of the form (zzxx + zzyy+xxzz+yyzz), i.e., selecting two Kronecker symbols and two S x , S y terms. The Kronecker symbols give either p = k = 0 or q = k = 0 and thus one of the cos-terms vanishes: cos p x − cos p y = 0 or cos q x − cos q y = 0.
We finally want to note an interesting observation. If one uses the well-known relations for spin-1/2 operators S i · S j 2 = 3 16 − 1 2 S i · S j and S i · S j S i · S k = 1 4 S j · S k + i 2 S i · (S j × S k ), one may rewrite the biquadratic term as a sum of three spin exchange terms such that the spin Hamiltonian in Eq. (S.50) takes the form (S.55)