Nematic and spin-charge orders driven by hole-doping a charge-transfer insulator

Recent experimental discoveries have brought a diverse set of broken symmetry states to the center stage of research on cuprate superconductors. Here, we focus on a thematic understanding of the diverse phenomenology by exploring a strong-coupling mechanism of symmetry breaking driven by frustration of antiferromagnetic order. We achieve this through a variational study of a three-band model of the CuO$_2$ plane with Kondo-type exchange couplings between doped oxygen holes and classical copper spins. Two main findings from this strong-coupling multi-band perspective are 1) that the symmetry hierarchy of spin stripe, charge stripe, intra-unit-cell nematic order and isotropic phases are all accessible microscopically within the model, 2) many symmetry-breaking patterns compete with energy differences within a few meV per Cu atom to produce a rich phase diagram. These results indicate that the diverse phenomenology of broken-symmetry states in hole-doped antiferromagnetic charge-transfer insulators may indeed arise from hole-doped frustration of antiferromagnetism.

Closely linked to the question of whether a thematic understanding of the observed phenomena is attainable is a theoretical question of whether to take the weak-coupling Fermi-surface instability perspective or to take the strong-coupling perspective. From a Fermi-surface instability perspective, the simultaneous occurrence of multiple orders requires fine tuning as one usually finds a single dominant instability in one ordering channel. For instance, it was shown that antiferromagnetic exchange interactions [23] or short-range repulsive interactions [24] can drive an instability towards density modulation along an incommensurate vector Q = (Q 0 , Q 0 ) with a dominantly d-form factor. Although a Q = 0 order with the d-form factor will be equivalent to IUC nematic order observed in various experiments [15,16], a modulation at finite Q will not show net IUC nematicity, just as antiferromagnets have no net magnetization.
Although it is well known that motion of doped holes will frustrate the antiferromagnetic background of the cuprate "parent compounds" [45], much of the work on this issue has largely focused on the one-band Hubbard model as a minimal framework to discuss this physics. Nevertheless, a host of experimental observations on cuprate superconductors (e.g., electron-hole doping asymmetry of the phase diagram and unusual broken-symmetry phases with spin and charge order [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] at low hole doping) strongly call for a description which explicitly retains the oxygen orbitals. Such multi-orbital models of the CuO 2 plane, like the Emery model [46], are, however, even less amenable to a theoretical treatment than the one-band Hubbard model.
In this paper, we consider a simplified model which ignores the charge fluctuation on the transition metal sites and treat the spins on those sites as classical local moments that interact with doped itinerant holes living on the oxygen sites through a Kondo type coupling (see Figure 1). This is a tractable three-orbital model that retains the spatial separation between the local moments on the transition metal ion such as Cu or Ni and the doped holes that predominantly go into oxygen sites. We consider several plausible spin-order patterns and exactly solve for the many-hole eigenstates associated with the spin-order patterns. This way we can access not only the lowest energy spinhole configuration, but investigate how magnetism drives charge physics in the cuprates, study the energetic competition between a wide range of spin and charge ordered states, and gain insight into the interplay between spin and intra-unit-cell nematic orders [15].
The rest of the paper is organized as follows. In section 2 we introduce and motivate the model as well as the choice of some parameters. In section 3 we discuss the choice of spin-ordering configuration ansätze. In section 4 we discuss the charge-ordered states . This choice of parameters leads to a minimum at the nodal (±π/2, ±π/2) points in agreement with ARPES experiments on the extremely underdoped cuprate [47].
we obtain associated with each spin-ordering configuration. In section 5 we discuss the energy differences between different ansätze and the phase diagram. In section 6 the effect of a nearest-neighbor oxygen-oxygen interaction is analyzed. Finally, we conclude with a summary and discussion in section 7.

Model
In this paper we restrict our attention to the Hilbert space with singly occupied Cu d x 2 −y 2 orbitals represented as local moments that interact antiferromagnetically. The doped holes will be assumed to go into the oxygen p x or p y orbitals [46]. These holes couple to the local moments on Cu sites through both (spin-dependent) hopping between O sites mediated by Cu sites and spin-spin interactions. Treating the moments on the Cu sites classically, we arrive at a model Hamiltonian wherep Is are the oxygen hole annihilation operators at oxygen site I with spin s, S i is a classical spin vector at Cu-site i with | S i | = 1/2, and σ are the Pauli matrices. The first two terms represent two hopping processes through the Cu sites, the t a (t b ) process with hole spin anti-parallel (parallel) to S i . These hopping processes amount to nearest-and next-nearest-neighbor hopping on the lattice of oxygen sites. Note that the second term in equation (1) hence introduces coupling between the motion of the holes on O-sites and the spin on the Cu sites. The third term represents direct hopping between O-sites, which amounts to nearest-neighbor hopping on the lattice of oxygen sites. The last two terms represent exchange coupling between Cu spins (the J 1 term) and between a Cu spin and the hole spin on a neighboring O site (the J 2 term). The model Hamiltonian of equation (1) is a simplified form of a perturbative expansion in Cu-O hopping t pd (to order t 4 pd ) in the large Coulomb interaction limit of the Emery model Hamiltonian for the cuprates [46]: However, the actual model obtained through such perturbative expansion is still highly non-trivial as it does not allow double occupancy at oxygen sites and the spins on Cusites should be quantum mechanical spins. Frenkel et al [48] studied the exact quantum ground state of such a model in the presence of a single hole in a small cluster. Despite the reduction of the Hilbert space due to the constraint of singly occupied Cu d x 2 −y 2 orbitals, the largest cluster they could study was a 4 × 4 Cu-O cluster which would be too small to see the observed stripe phenomena. Such studies have been extended to larger clusters [49] with the aim of understanding spin polaron formation in a 3-band model. However, this exact diagonalization study is restricted to 1-hole and 2-hole states in the undoped insulator and thus cannot address doping dependent competing orders. By treating the local moments on Cu sites to be classical and relaxing the no-double-occupancy constraint on oxygen sites, we arrive at the more tractable model Hamiltonian of equation (1). Our simple model has several features. Firstly, it is a solvable model that retains much of the microscopic details and non-trivial interactions of under-doped cuprates. The solvability of the model allows us to study a wide variety of states including incommensurate orderings. Secondly, we expect similar effective models are applicable to other doped strongly-correlated charge-transfer insulators such as the nickelates. In particular, a classical approximation of the local moment on the transition-metal ion is likely to be a better approximation in Ni given its larger moment. Finally, through spatial separation between spins and holes, the model offers a rich playground for strongcoupling-driven spin and charge orders co-existing with intra-unit-cell nematic order.
Earlier consideration of the strong-coupling limit of the Emery model by Kivelson et al. [50] focused on the limit of vanishing inter-oxygen site hopping. In that limit, the dynamics is strictly one-dimensional and hence the ground state is a nematic phase. Three key differences between the limit considered in Ref. [50] and the limit captured by our model equation (1) are that 1) we have taken U d to be so strong to the extent that we suppressed the charge fluctuation in the Cu sites, 2) we consider inter-oxygen hopping to be comparable to exchange interactions, 3) we take the Cu-O exchange interaction into account. As a result, the charge dynamics in our model is strictly two dimensional. However we make the approximation of leaving out U p and V pp which yields an exactly solvable model for the hole motion on a system size that can exhibit translational symmetry breaking. This approximation may not affect the conclusions qualitatively in the limit of dilute hole density.
Although the model Hamiltonian has quite a few parameters, many of them are constrained by experiments. We use a Cu-Cu exchange interaction J 1 ∼ 125meV, close to the value estimated for La 2 CuO 4 from neutron scattering studies of the spin-wave dispersion [51]. In order to constrain hopping strengths, we examine the dispersion of a single hole doped in an ordered Néel antiferromagnet. We can diagonalize the hole Hamiltonian assuming an out of plane spin moment on Cu atoms in an antiferromagnetic arrangement, i.e. S z i = 1 2 (−1) ix+iy , and obtain four dispersing bands (see Appendix A). The resulting lowest dispersion can be compared to the ARPES measurements on extremely underdoped cuprates [47]. Although the ARPES energy distribution curves are broad, we can infer the following energy scales from the peak positions. The binding energy at (π/2, π/2) is lower than that at (π, 0) by about 300meV [47], which corresponds to 2t pp ; this fixes t pp = 150meV. Shen et al. [47] also found that the spectral peaks disperse by about 1.4eV going from the Γ-point to (π/2, π/2). In our model, this energy difference is just 4t b + 2t pp ; this fixes t b = 275meV. Cluster diagonalization calculations in Ref. [48,52] indicate that t a is the same sign as We therefore examine a range of values t b /2 < t a < 3t b /2. The above choice of hopping parameters yields a Fermi surface consisting of hole-pockets centered at (±π/2, ±π/2) for thep Is holes (see Fig. 1(b)), consistent with the energy minimum of a single hole in the Néel ordered antiferromagnet being centered at (±π/2, ±π/2). We expect the antiferromagnetic Cu-O exchange interaction J 2 , which is only present at finite doping, to be stronger than the Cu-Cu exchange J 1 . Hence we explore a range of J 1 < J 2 < 3J 1 . We note that our values for the parameters in the effective Hamiltonian for doped holes differ slightly from those used in previous work [48,52,53] -the values we use should be viewed as effective couplings given our assumption of a classical copper spin. They also place the regime of interest in the strong coupling limit, since t a + t b and J 2 are much greater than the free fermion bandwidth set by |t a − t b | ∼ t pp and motivate the use of a variational approach as discussed in the next and following sections.

Ansätze for spin-order configurations
The key mechanism by which the model Hamiltonian of equation (1) drives spin and charge order at finite doping is the frustration of antiferromagnetic order in the parent compound (i.e., zero doping). The undoped cuprate and nickelate insulators are ordered antiferromagnets, possessing long-range Néel order in planes which are stacked along the c axis. The planar Néel order has a wavevector (π, π) and is a collinear state with S i = Sn(−1) ix+iy withn being a unit vector. Doped oxygen holes frustrate the antiferromagnetic alignment of Cu moments driven by the Cu-Cu exchange J 1 through a strong antiferromagnetic Cu-O exchange interaction J 2 . Thus, doped holes promote a ferromagnetic arrangement of the two neighboring local Cu-moments and drive spatial symmetry breaking in the spin configuration. Interestingly, such frustration is not expected in the electron-doped cuprates, where doped electrons go onto Cu sites. This naturally explains why the Néel order is present to much higher doping in electron doped cuprates whereas the hole-doped cuprates exhibit various broken-symmetry states.
While ultimately the model equation (1) can be exactly solved by combining Monte Carlo simulation of the spin problem with the exact diagonalization of the fermion problem, we here consider two classes of spin-order configurations and solve the fermion problem within each class. This approach offers us an understanding of the parameterspace landscape. Moreover, although our understanding is limited to the set of ansatz spin patterns we consider, this approach has the advantage over the explicit solution in that it allows us to compare energetics of different candidate states. The two classes we consider are coplanar spin spirals and collinear spin stripes with different wave vector Q.
As there is much literature on both of these candidate states we defer in-depth discussion to review articles (see e.g., [54,55]). Rather, we discuss aspects of these states that are directly relevant to our study and the rationale for their consideration below.
(i) Coplanar spin spirals: While a doped hole promotes ferromagnetic correlations of the local moments in its vicinity, the spins far from the hole maintain their Néel correlations. Hence, in the absence of magnetic anisotropy a coplanar spiral state may be expected as a solution that interpolates between the Néel order and ferromagnetism when holes are doped into a quantum antiferromagnet as it has been argued in Ref. [56]. However, spiral states are accompanied by translationally invariant charge distribution and they are unstable against magnetic anisotropy as well as phase separation. Further, susceptibility studies [57] in spin-ordered cuprates are consistent with collinear order. However, the smoking-gun polarized neutron scattering experiment has not been performed to date and local spiral distortions of the spin background may be an appropriate picture in the insulating regime at low doping [58].
Below, we pick the spiral to lie in the S x − S y plane, setting S i = S(cos Q · r i , sin Q · r i , 0), and consider different propagation wavevectors Q running parallel to the Cu-O bonds (labelled Psp1, Psp2) or along the diagonal direction (Dsp), see Fig. 2, and compute their physical properties and energies. We show that within a three-orbital picture, the spin spirals with a wavevector parallel to the Cu-O bonds are naturally accompanied by a translationally invariant oxygen hole density, but with charge nematic order on the oxygen sites as has been observed in STS experiments [15]. This opens up a possibility of disordered spiral leaving the discrete and robust nematic order as the only observable effect as it has been proposed for various frustrated magnets [59].
(ii) Stripes: Two theoretical approaches discussed stripes prior to their experimental observations. The first approach was that of Hartree-Fock mean-field theory [27][28][29][30] which predicted insulating (fully filled) period-8 charge stripes at x = 1/8 doping. The other approach was built on the observation that phase separation should be expected in t-J models for J t, where 'slow' doped holes cluster together to avoid disrupting the Néel order of the local moments [60]. This led to the proposal that stripe orders could emerge as periodic density modulations due to frustration of such phase separation by long-range Coulomb interactions [44]. Following the experimental observation of metallic period-4 charge stripe at x = 1/8 doping in La-based cuprates [2], theoretical efforts focused on going beyond Hartree-Fock mean-field theory to obtain the observed metallic stripe [31,[39][40][41][42][43]. A particularly convincing case was made by a density matrix renormalization group study by White and Scalapino [31] which not only obtained the observed metallic stripe ground states at x = 1/8 but also demonstrated that anti-phase domain walls between collinear antiferromagnetic domains attract charge stripes. This issue of charge stripe periodicity (period-4 and metallic v.s. period-8 and insulating) bears particular importance in the discussion of the role of charge stripes in superconductivity [32]. However, most previous efforts at going beyond Hartree-Fock mean-field theory and incorporating strong-coupling physics relied on sophisticated numerics. The spatial separation of Cu-moment and doped holes in our model allows us to investigate the complex spin-charge interplay in the strong-coupling regime in a transparent manner. In our model, spin-antiphase domain walls can form upon hole-doping to relieve the frustration between the AFM order favored by Cu-Cu antiferromagnetic superexchange interactions and the ferromagnetic order favored by an antiferromagnetic Cu-O exchange.
We model the spin stripes as antiphase domain walls in the collinear antiferromagnet and consider several different uni-directional domain-wall configurations depending on the direction of periodicity as well as the period. Note that these spin stripes will be bond-centered by construction. We consider parallel stripes with wave vectors Q along the Cu-O bond directions as well as diagonal stripes with Q at 45 o angle with respect to Cu-O bond direction. At fixed hole density of x = 1/8, different periods determine whether the stripe can support metallic transport along the stripe: period-4 parallel charge stripe (labeled PM, see Fig 2) will be half-filled and metallic while period-8 parallel charge stripe (labeled PI, see Fig 2) will be fully-filled and insulating. Similarly, one can consider metallic (labeled DM, see Fig 2) and insulating (labeled DI, see Fig 2) diagonal stripes.

Charge orders associated with different spin-order ansätze
For each spin-order candidate, we diagonalize the (quadratic) hole Hamiltonian and use the resulting lowest-energy many-hole wave function and the spin configuration to evaluate the total energy for the spin-charge-ordered state. The total energy for a given spin configuration { S} at a given doping x = n/N , with N the number of Cu-O unit cells, is given by where ξ l ({ S}) are the (energy-ordered) eigenvalues of the Hamiltonian (1) for the given configuration { S}. In equation (4) the last term is the interaction energy due to Cu-Cu  is a diagonal spiral. The parallel spiral is accompanied by IUC nematic order with η = (n x − n y )/(n x + n y ) ≈ 5% for the above configuration. Note that the length of the (spin) arrow on the oxygen site is scaled by a factor of 5 as compared to the Cu sites for illustration purposes.
exchange (the J 1 term in equation (1)). For the spiral configurations given by S i = S(cos Q·r i , sin Q·r i , 0) we take advantage of a closed form of the Hamiltonian equation (1) in momentum space (See Appendix A) to find the wave vector Q that minimizes the energy for each spiral. Fig. 3 shows the lowest-energy charge distribution for a parallel spiral (Fig. 3(a)) and a diagonal spiral (Fig. 3(b)) spin configurations. Though it is well known that spiral order of spins would not be accompanied by any charge modulation, Fig. 3(a) shows that parallelspiral tendency drives IUC nematic charge order. Note that even when the spiral order of the spins is disordered due to thermal or quantum fluctuations, the discrete symmetry breaking of IUC nematic in the charge sector can be more robust.
For the collinear stripe configurations, we consider a lattice of 32 × 32 unit cells (2048 oxygen sites) and diagonalize the quadratic Hamiltonian of the holes only living on the oxygen site, equation (1). The resulting charge-order patterns for different trial spin  configurations (see Fig. 4) clearly show that the holes are attracted to the anti-phase domain walls in the spin configurations driven by the Cu-O exchange-coupling J 2 . The kinetic terms broaden the hole distribution and favor the metallic charge stripe. It is remarkable that our simple model can readily access the spin and charge striped groundstate configuration reminiscent of those obtained in density matrix renormalization group studies of the 1/8-doped t-J model [31]. Moreover, as this model incorporates the mostly-oxygen character of doped holes, the charge stripes centered at the antiphase domain wall of the antiferromagnetic background are naturally coupled to IUC nematic. The parallel stripe configurations obtained in Fig. 4 clearly demonstrate a coupling between the (Ising) IUC nematic order and the stripe-ordering wave vector.
In the sense of Ginzburg-Landau theory of order parameters, Fig. 4

Phase diagram
While our variational study is limited by the choice of states we consider, it has the advantage of allowing the energetic comparison between different candidate states over Hartree-Fock or variational Monte Carlo studies. Surprisingly, we find the entire collection of states we consider to show close energetic competition. For the most part of the parameter space we consider (varying hopping through Cu t a and Cu-O exchange J 2 ) different states differ in energy only by a few meV. The close energetic competition is clear in Figure 5 which shows energy differences as a function of t a at fixed J 2 = 170meV. Such close competition indicates small perturbations to our model Hamiltonian such as spin-orbit coupling or lattice anisotropies could change the phase diagram significantly.   Fig. 2) compared to the diagonal spiral state (Dsp) as a function of t a for J 2 = 170meV (dashed line in the phase diagram in Figure 6(b)). Shown are the antiferromagnetic (AFM), the parallel spiral (Psp1), the diagonal spiral (Dsp), as well as the diagonal metallic stripe (DM), diagonal insulating stripe (DI) parallel metallic stripe (PM), and parallel insulating stripe (PI). Notice the parallel orders are preferred for larger t a .
It is well known that incorporating quantum fluctuations of spins strongly favors collinear order over coplanar or non-coplanar states [62,63], Similarly, spin anisotropy due to weak spin-orbit coupling could also hinder coplanar spirals. We therefore present two phase diagrams here: a phase diagram that includes co-planar spiral and a phase diagram that only concerns collinear orders. Figure 6(a) shows the phase diagram that allows for co-planar spiral as a function of hopping through Cu-sites t a and the Cu-O exchange J 2 . We have limited the plotting range of J 2 to be a reasonable range of 100meV < J 2 < 300meV. Figure 6(a) shows that the lowest-energy configuration with the parameters motivated from cuprates is dominated by diagonal orders such as diagonal spiral or diagonal metallic stripe. In the limit of large J 2 (not shown), parallel spiral with IUC nematic appears before Cu spins align ferromagnetically. Figure 6(b) shows the phase diagram restricted to collinear orders. Notably metallic diagonal and parallel stripe orders appear in a substantial region of the phase diagram. The close competition between diagonal stripe order and parallel stripe order is remarkable in light of the experimentally known transition from diagonal to parallel stripe upon doping [64]. The parallel metallic charge stripe order (PM) at x = 1/8 is consistent with observations in La-based compounds [64,65] as well as those in Bi-based compounds [5,6,15,61]. The microscopic tie between spin-ordering patterns and charge-ordering patterns resulting from our simple model depicted in Figs. 3 and 4 hints at a microscopic mechanism of charge order without spin order at finite temperature [66]. The spin orders under consideration require breaking of spin rotational symmetry as well as spatial translation and point-group symmetries while charge orders only break spatial symmetries. While the driving force for charge order in our model is the frustration of antiferromagnetic order upon holes entering oxygen sites, fluctuations in spin space will suppress spin order. It has been shown in the context of Fe-based superconductors [67] that spin fluctuations | S Q | 2 = 0 in the absence of spin order ( S Q = 0) can drive nematicity. Similarly, since the cuprates are quasi-2D systems, it is quite plausible that thermal fluctuations may prevent spin order in the stripe state, while | S Q | 2 = 0 may still leave the charge order visible. This can be easily seen by the fact that the charge stripe and IUC nematic in the parallel stripe phases in Fig. 4 as well as IUC nematic in the parallel spiral phase in Fig. 3 are insensitive to the spin-orientation, so (thermal) averaging over the global spin orientation will leave these orders unaffected. Hence, charge order in Fig. 3 can onset without detectable spin order or with a lower temperature transition into a state with coexisting spin order. This would amount to a microscopic realization of a so-called "charge-driven transition" in the Landau theory study Ref. [68]. Experimental observations of charge and/or spin order in cuprate families broadly show such preferential visibility of charge order [5,6,[8][9][10]15,61,64,65].

Effects of inter-oxygen repulsion V pp
As IUC nematic naturally accompanies all modulations along the Cu-O bond directions within our model, it is natural to ask what would be the effect of the inter-oxygen repulsion V pp that was found to drive IUC nematic [38,50]. Recently Bulut et al. [24] extended the mean-field study by two of us [38] to include the charge-density wave instability as well as IUC nematic instability. They found the d-form factor component ψ d (r), defined on oxygen sites using the notation of Ref. [23] via where f s (r) (f s (r)) is zero (one) on copper sites and one (zero) on oxygen sites and f d (r) is zero on copper sites, one on x oriented oxygen sites and minus one on yoriented oxygen sites, to dominate the density wave instability along the Brillouin zone diagonal. While different form-factor components of density waves are not symmetry distinct, experimental observation of the d-form factor [5,69] hints at the importance of microscopic interactions promoting an anti-phase relation between the two oxygen sites [69]. As charge fluctuations on Cu-sites have been projected out in our model (ψ s (r) = 0), charge stripes found in our model consist only of s'-and d-form factor components. In order to study the effect of V pp on the charge configurations associated with candidate spin configurations, we treat the V pp term at the level of self-consistent Hartree approximation. As experiments observe Q = 0 IUC nematic simultaneously with shortranged d-form factor density waves [69] we focused the study to the parameter space where a Cu-O bond-direction charge stripe is present in the absence of V pp .
On symmetry grounds, ψ s (r), ψ d (r) and IUC nematic order parameter can form a cubic coupling within Landau theory. Hence, we generically expect both s'-and dform factor components to be present given the robust IUC nematicity coexisting with charge stripes in our charge-ordering patterns. For t a = 350meV, J 2 = 170meV, figure  7 shows the resulting ratio of the d-and s'-form factor components. While the density waves we obtained for the PM configuration have predominantly s'-form factor, Fig. 7 clearly shows that V pp promotes a d-form factor for the charge stripe along the Cu-O bond direction with the wave vector consistent with the experimental observations.

Discussion
To summarize, we considered a three-orbital model for underdoped cuprates, which incorporates strong-coupling physics through singly occupied Cu-sites hosting local moments and which reflects the charge-transfer energy through constraining holes to live on the oxygen sites. We took a variational approach for the spin configurations and solved for many-hole states exactly for each configuration at x = 1/8 doping. We found that the balance between the Cu-Cu exchange interaction, the kinetic energy of holes and the Kondo type coupling between the Cu spins and the spin of the holes conspire to a rich phase diagram featuring several experimentally observed phases. What is more interesting is the close energetic competition between the candidate states such as coplanar spirals, diagonal stripes and parallel stripes. Within this model, IUC nematic accompanies any parallel order including co-planar spiral orders, providing a microscopic mechanism for robustness of a IUC nematic.
Though we focused on the parameter space constrained by various experiments on the cuprates, our model can be applied to other transition metal oxides like the doped nickelates which exhibit similar spin-charge ordered states as cuprates [1]. In the nickelates, a local moment at Ni sites may arise from the effect of strong Hund's coupling between two electrons (or two holes) in the two e g orbitals, leading to an orbital singlet and a spin triplet S = 1 state. This spin-triplet nature of the local moment suppresses hopping through Ni sites which renders the system insulating at all doping (x < 0.9) without ever exhibiting superconductivity. Antiferromagnetism is more robust in Nickelates and extends to hole concentration of n h = x + 2δ ∼ 0.2 before it is replaced by fully filled diagonal stripe order [1]. Indeed, we find the diagonal insulating stripe as most favored co-linear order in the limit t a J 1 < J 2 . Recently, several model studies investigating the Fermi-surface instability towards formation of charge-density waves [23,24,[70][71][72][73] have found charge-density waves with predominantly d-form factor. Some theories focused on the "hot spot" regions of the Fermi surface, where antiferromagnetic scattering will be especially strong near the antiferromagnetic quantum critical point [23,70], while others considered one-band models with d-form factor interactions [71,72] or a three-band model [24]. Though the prominence of the d-form factor component in the density waves appear compatible with experiments [5,69], these models universally obtained charge-density waves with the wave vector along the Brillouin zone diagonal, i.e., diagonal stripes. This direction and the magnitude of the wave vector is at odds with experimental findings of parallel stripes. More recently, Atkinson et al [73] showed that the wave vector found in Ref. [24] can be rotated to lie in the Cu-O bond direction by assuming large staggered moments which reconstruct the Fermi surface. However, experimentally it is known that antiferromagnetic fluctuations at (π, π) are replaced by incommensurate magnetic peaks at temperatures well above the charge ordering temperature in inelastic neutron scattering at low energies [64,74].
In this paper, we took a strong-coupling approach of projecting out charge fluctuation at the Cu-sites yet dealing with an exactly solvable model by ignoring quantum fluctuations of spins at the Cu-sites. In this approach, we found Q = 0 IUC nematic to naturally coexist with any modulation along the Cu-O direction be it co-planar spiral or collinear spin stripes. This is in contrast with weak-coupling Fermi-surface-instability approaches finding competition between Q = 0 order and Q = 0 density waves. Also, we find various experimentally observed charge-ordered states (diagonal insulating stripe, diagonal metallic stripe, parallel metallic stripe) to dominate at different parts of the physically-motivated parameter space exhibiting a close competition. The understanding of the parameter space of our model we gained from the present variational approach can guide us in the attempt for a more explicit exact solution of the model combining Monte Carlo simulations with exact diagonalization.