# Monte Carlo study of cuprate superconductors in a four-band $d$-$p$ model:
Role of orbital degrees of freedom

Hiroshi Watanabe1<EMAIL_ADDRESS>Tomonori Shirakawa2,3 Kazuhiro
Seki3 Hirofumi Sakakibara4,5,6 Takao Kotani4,5 Hiroaki Ikeda7 Seiji
Yunoki2,3,6,8 1Research Organization of Science and Technology, Ritsumeikan
University, Shiga 525-8577, Japan
2Computational Materials Science Research Team, RIKEN Center for Computational
Science (R-CCS), Hyogo 650-0047, Japan
3Quantum Computational Science Research Team, RIKEN Center for Quantum
Computing (RQC), Saitama 351-0198, Japan
4Advanced Mechanical and Electronic System Research Center (AMES), Faculty of
Engineering, Tottori University, Tottori 680-8552, Japan
5Center of Spintronics Research Network (CSRN), Graduate School of Engineering
Science, Osaka University, Osaka, 560-8531, Japan
6Computational Condensed Matter Physics Laboratory, RIKEN Cluster for
Pioneering Research (CPR), Saitama 351-0198, Japan
7Department of Physics, Ritsumeikan University, Shiga 525-8577, Japan
8Computational Quantum Matter Research Team, RIKEN Center for Emergent Matter
Science (CEMS), Saitama 351-0198, Japan

###### Abstract

Understanding the complex phase diagram of cuprate superconductors is a long-
standing challenging problem. Recent studies have shown that orbital degrees
of freedom, both Cu $e_{g}$ orbitals and O $p$ orbitals, are a key ingredient
for a unified understanding of cuprate superconductors, including the material
dependence. Here we investigate a four-band $d$-$p$ model derived from the
first-principles calculations with the variational Monte Carlo method, which
allows us to elucidate competing orders on an equal footing. The obtained
results can consistently explain the doping dependence of superconductivity,
antiferromagnetic and stripe phases, phase separation in the underdoped
region, and also novel magnetism in the heavily-overdoped region. Our four-
band $d$-$p$ model with neighbouring intersite interactions is a minimal model
to describe the phase diagram comprehensively. The presence of $p$ orbitals is
critical to the charge-stripe features, which induce two types of stripe
phases with $s^{\prime}$-wave and $d$-wave bond stripe. On the other hand, the
presence of $d_{z^{2}}$ orbital is indispensable to material dependence of the
superconducting transition temperature ($T_{\mathrm{c}}$), and enhances local
magnetic moment as a source of novel magnetism in the heavily-overdoped region
as well. These findings beyond one-band description could provide a major step
toward a full explanation of unconventional normal state and high
$T_{\mathrm{c}}$ in cuprate supercondutors.

##  I INTRODUCTION

Over 35 years since its discovery [1], cuprate superconductors have
continuously challenged our conventional understandings, such as the recent
discovery of nematic transitions in the pseudogap region [2]. It has not yet
been achieved to consistently explain the whole phase diagram and the
correlation between competing orders and the high transition temperature
($T_{\mathrm{c}}$). The key to explain the features is considered to be
multiorbital effects. The importance of the orbital degrees of freedom is of
great interest in modern condensed matter physics as a source of emergent
phenomena such as spin currents [3, 4] and, in the field of superconductivity,
as a source of novel pairing states and high $T_{\mathrm{c}}$ [5, 6].

In cuprate superconductors, for a long time, the anomalous features have been
investigated as the physics of an effective single band crossing a Fermi
surface rather than multiorbital effects [7]. The effective one-band models,
such as the Hubbard model and the $t$-$J$ model, successfully predicted
$d$-wave superconductivity, but were insufficient to describe the material
dependence of $T_{\mathrm{c}}$ and the unconventional competing orders. The
$d$-$p$ model or the Emery model [8], which consists of Cu $d_{x^{2}-y^{2}}$
orbital and O $p_{x}/p_{y}$ orbitals, was studied early on as a model
involving the multiorbital effect. These models have been intensively studied
in terms of the material dependence of $T_{\mathrm{c}}$, pseudogap phenomena,
stripe features, and so on [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. Although these models
have partially captured the unconventional features, it could not fully
explain the anomalous features of cuprates. For example, the large difference
in $T_{\mathrm{c}}$ between $\sim$40 K for a La-based system and $\sim$90 K
for a Hg-based system remained unclear. Recently, one of the authors,
Sakakibara and co-authors suggested the importance of the $d_{z^{2}}$ orbital
based on the analysis of a two-orbital model [36]. The importance of the
$d_{z^{2}}$ orbital has also been supported by the latest angle-resolved
photoemission spectroscopy experiment [37].

These studies indicate that both $p$\- and $d$-orbital degrees of freedom
should be properly taken into account for the understanding of cuprate
superconductors. In the previous work, we have proposed a four-band $d$-$p$
model composed of the Cu $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals and the O
$p_{x}$ and $p_{y}$ orbitals as a minimal model to obtain the unified
description of cuprate superconductors [38]. On the basis of the variational
Monte Carlo (VMC) method, we have shown that this model explains well two key
factors about the material dependence of $T_{\mathrm{c}}$: the contribution of
the Cu $d_{z^{2}}$ orbital to the Fermi surface and the site-energy difference
$\Delta_{dp}$ between the Cu $d_{x^{2}-y^{2}}$ and O $p$ orbitals [39].

In this paper, we investigate the four-band $d$-$p$ model with the VMC method
in more detail. We take the La2CuO4 and HgBa2CuO4 systems as typical examples,
and especially, elucidate the competing orders of superconductivity,
antiferromagnetic (AF) and stripe phases, phase separation in the underdoped
region, and novel magnetism in the heavily-overdoped region.

Our major findings are as follows. First, charge/spin stripe state is stable
over a wide doping range, and its period decreases as the hole doping rate $x$
increases. At $x=1/8$, we obtain the same stripe phase as observed
experimentally, which is extremely robust in the La-based system, but fragile
in the Hg-based system. The AF phase is confined to a narrow doping range near
$x\sim 0$. In the underdoped region $0<x<1/8$, the phase separation (PS)
occurs between the AF and the $x=1/8$ stripe phase. Suppression of the AF
correlation by finite intersite $d$-$d$ repulsion causes the transition to a
$d$-wave bond stripe observed in non-La-based systems. Second, concerning the
superconductivity, the material dependence of the $T_{\mathrm{c}}$ dome as a
function of $x$ is consistently explained. The dome shape is shown to be
strongly correlated with whether the undoped AF insulator is Slater- or Mott-
type. Finally, we show that novel magnetism observed in the heavily-overdoped
region comes from the development of Cu local moments via the Hund’s coupling
due to the rapid increase of $d_{z^{2}}$ component with doping.

The rest of this paper is organized as follows. In Sec. II, we introduce the
model and the numerical method used in this paper. The four-band $d$-$p$ model
on the two-dimensional square lattice is introduced in Sec. II.1. The VMC
method and the Gutzwiller-Jastrow type trial wave function are explained in
Sec. II.2. The tight-binding energy bands for the La2CuO4 and HgBa2CuO4
systems, obtained on the basis of the first-principles calculation, are shown
in Sec. II.3. The numerical results are provided in Sec. III. The energy
competition between AF and stripe phases are studied in detail and the ground
state phase diagram including the PS is shown in Sec. III.1. The symmetry
change within the stripe phase is also discussed. The material and Coulomb
interaction dependence of superconductivity is shown in Sec. III.2 through the
behavior of superconducting correlation functions. The difference between La-
and Hg-based systems is discussed from the band structure and electron
correlation effects. The behavior of $d_{z^{2}}$ hole density and the local
magnetic moment is studied in Sec. III.3. The novel magnetism in the heavily-
overdoped region is also discussed. Finally, the paper concludes with a
summary in Sec. IV. The details of the variational wave functions are shown in
Appendix A.

Figure 1:  (a) Schematic lattice structure of the four-band $d$-$p$ model,
forming the two-dimensional square lattice. Each Cu site contains
$d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals, while there is either $p_{x}$ or
$p_{y}$ orbital on O sites. The lattice constant between the nearest-neighbor
Cu sites is set to be one, i.e., primitive translation vectors
$|\textbf{e}_{x}|=|\textbf{e}_{y}|=1$. The dotted square represents a unit
cell. (b) The Coulomb interaction parameters between nearest-neighbor
orbitals.

##  II MODEL AND METHOD

###  II.1 Four-band $d$-$p$ model

As mentioned in Introduction, it is important to incorporate the four
orbitals, the Cu $e_{g}$ orbitals and the O $p_{x}/p_{y}$ orbitals, to provide
a unified description of the cuprate superconductors. Therefore, as an
effective low-energy model, we consider a four-band $d$-$p$ model on the two-
dimensional square lattice [see Fig. 1(a)] defined by the following
Hamiltonian:

$H=H_{\text{kin}}+H_{\text{int}}-H_{\text{dc}}.$ (1)

First, the kinetic term $H_{\text{kin}}$ is described by

$\displaystyle H_{\text{kin}}$
$\displaystyle=\sum_{i,j,\sigma}\sum_{\alpha,\beta}t^{\alpha\beta}_{ij}c^{\dagger}_{i\alpha\sigma}c_{j\beta\sigma}$
(2)
$\displaystyle=\sum_{\textbf{k},\sigma}\sum_{m}E_{m}(\textbf{k})a^{\dagger}_{\textbf{k}m\sigma}a_{\textbf{k}m\sigma},$
(3)

where Eq. (2) is the kinetic term in an orbital representation and Eq. (3) is
in a band representation. $c^{\dagger}_{i\alpha\sigma}$ ($c_{i\alpha\sigma}$)
is a creation (annihilation) operator of an electron at site $i$ with spin
$\sigma\,(=\uparrow,\downarrow)$ and orbital $\alpha\,(=1,2,3,4)$
corresponding to ($d_{x^{2}-y^{2}}$, $d_{z^{2}}$, $p_{x}$, $p_{y}$),
respectively. $t^{\alpha\beta}_{ij}$ denotes a hopping integral between
orbital $\alpha$ at site $i$ and orbital $\beta$ at site $j$.
$t^{\alpha\alpha}_{ii}$ is a site energy $\varepsilon_{\alpha}$ for orbital
$\alpha$ at site $i$. These hopping integrals and site energies are determined
from the first-principles calculations (see Sec. II.3). Equation (3) is
obtained by diagonalizing Eq. (2), and the energy eigenvalue
$E_{m}(\textbf{k})$ is characterized by the wave vector k and the energy band
index $m\,(=1,2,3,4)$. $a^{\dagger}_{\textbf{k}m\sigma}$
($a_{\textbf{k}m\sigma}$) is a creation (annihilation) operator of the
corresponding energy band with spin $\sigma$. The undoped parent compounds of
cuprates correspond to one hole (i.e., seven electrons) per unit cell in this
model, which is conventionally referred to as half filling, and hereafter we
denote the carrier density as the number $x$ of holes per unit cell that are
introduced into the system at half filling.

Second, the Coulomb interaction term $H_{\text{int}}$ is composed of eight
terms,

$\displaystyle H_{\text{int}}$
$\displaystyle=U_{d}\sum_{i}\left(n^{d_{1}}_{i\uparrow}n^{d_{1}}_{i\downarrow}+n^{d_{2}}_{i\uparrow}n^{d_{2}}_{i\downarrow}\right)$
$\displaystyle+\left(U^{\prime}_{d}-\frac{J}{2}\right)\sum_{i}n^{d_{1}}_{i}n^{d_{2}}_{i}-2J\sum_{i}\textbf{S}^{d_{1}}_{i}\cdot\textbf{S}^{d_{2}}_{i}$
$\displaystyle-J^{\prime}\sum_{i}\left(c^{\dagger}_{i1\uparrow}c^{\dagger}_{i1\downarrow}c_{i2\uparrow}c_{i2\downarrow}+c^{\dagger}_{i2\uparrow}c^{\dagger}_{i2\downarrow}c_{i1\uparrow}c_{i1\downarrow}\right)$
$\displaystyle+U_{p}\sum_{i}\left(n^{p_{x}}_{i\uparrow}n^{p_{x}}_{i\downarrow}+n^{p_{y}}_{i\uparrow}n^{p_{y}}_{i\downarrow}\right)+V_{dp}\sum_{\left<i,j\right>}n^{d}_{i}n^{p_{x/y}}_{j}$
$\displaystyle+V_{pp}\sum_{\left<i,j\right>}n^{p_{x}}_{i}n^{p_{y}}_{j}+V_{dd}\sum_{\left<i,j\right>}n^{d}_{i}n^{d}_{j}.$
(4)

Here, $n^{\alpha}_{i}=n^{\alpha}_{i\uparrow}+n^{\alpha}_{i\downarrow}$ with
$n^{\alpha}_{i\sigma}=c^{\dagger}_{i\alpha\sigma}c_{i\alpha\sigma}$ is the
number operator and $\textbf{S}^{\alpha}_{i}$ is the spin angular momentum
operator at site $i$ with orbital $\alpha$. $d_{1}$ and $d_{2}$ are
abbreviations for $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals, respectively,
and $n^{d}_{i}=n^{d_{1}}_{i}+n^{d_{2}}_{i}$. $U_{d},U_{d}^{\prime},J,$ and
$J^{\prime}$ represent on-site intraorbital, interorbital, Hund’s coupling,
and pair-hopping interactions between $d$ orbitals, respectively. In this
study, we set $J^{\prime}=J$ and $U_{d}=U^{\prime}_{d}+2J$ [40]. The on-site
intraorbital Coulomb interaction within $p$ orbitals, $U_{p}$, is also
introduced. The last three terms in $H_{\text{int}}$ take into account the
intersite Coulomb interactions between nearest-neighbor orbitals,
$V_{dp},V_{pp},$ and $V_{dd}$, as shown in Fig. 1(b), where the sum
$\sum_{\left<i,j\right>}$ runs over all pairs of nearest-neighbor orbitals
located at site $i$ and $j$. These Coulomb interactions are estimated from the
first-principles calculations (see Sec. II.3).

Finally, the following double counting correction term $H_{\text{dc}}$ is
introduced,

$H_{\text{dc}}=\biggl{[}\left\\{U_{d}+2\left(U^{\prime}_{d}-\frac{J}{2}\right)+16V_{dd}\right\\}\bar{n}^{d}\
\\\ +8V_{dp}\bar{n}^{p}\biggr{]}\sum_{i}n^{d}_{i}\\\
+\left\\{(U_{p}+8V_{pp})\bar{n}^{p}+8V_{dp}\bar{n}^{d}\right\\}\sum_{i}(n^{p_{x}}_{i}+n^{p_{y}}_{i}),$
(5)

where $\bar{n}^{d}$ and $\bar{n}^{p}$ are the average electron densities of
the $d$ and $p$ orbitals per spin per orbital obtained from the first-
principles calculation. When we apply a many-body calculation method to a
multiorbital model, the site energy of each orbital is shifted due to the
interaction effect. However, such energy shifts have already been included in
the energy band of the tight-binding model constructed from the first-
principles calculation. This is a so-called double counting problem and should
be treated with care especially in the $d$-$p$ model [41]. Here, we subtract
the term $H_{\text{dc}}$ from the Hamiltonian to correct the energy shift.
This is one of the reasonable treatments to avoid the double counting.

###  II.2 VMC method

In general, it is very difficult to treat the effect of Coulomb interactions
correctly. Here, the VMC method [42, 43, 44] was employed as a powerful
computational method. The trial wave function considered here is a Gutzwiller-
Jastrow type composed of four parts,

$\left|\Psi\right>=P^{(2)}_{\text{G}}P_{\text{J}_{\text{c}}}P_{\text{J}_{\text{s}}}\left|\Phi\right>.$
(6)

$\left|\Phi\right>$ is a one-body part constructed by diagonalizing the one-
body Hamiltonian including the off-diagonal elements
$\\{\rho^{\alpha}_{i}\\}$, $\\{m^{\alpha}_{i}\\}$, and
$\\{\Delta^{\alpha\beta}_{ij}\\}$, which induce long-range ordering of charge,
spin, and superconductivity, respectively. The renormalized hopping integrals
$\\{\tilde{t}^{\alpha\beta}_{ij}\\}$ are also included in $\left|\Phi\right>$
as variational parameters. The explicit forms of them are described in
Appendix.

The Gutzwiller factor

$P^{(2)}_{\text{G}}=\prod_{i,\gamma}\bigl{[}\text{e}^{-g_{\gamma}}\left|\gamma\right>\left<\gamma\right|_{i}\bigr{]}$
(7)

is extended to two-orbital systems [45, 46]. In $P^{(2)}_{\text{G}}$, possible
16 patterns of charge and spin configuration of the $d_{x^{2}-y^{2}}$ and
$d_{z^{2}}$ orbitals at each site $\left|\gamma\right>$, i.e.,
$\left|0\right>=\left|0\;0\right>$, $\left|1\right>=\left|0\uparrow\right>$,
$\cdots$,
$\left|15\right>=\left|\uparrow\downarrow\;\uparrow\downarrow\right>$, are
differently weighted with $\text{e}^{-g_{\gamma}}$ and $\\{g_{\gamma}\\}$ are
optimized as variational parameters. The remaining operators

$P_{\text{J}_{\text{c}}}=\exp\Bigl{[}-\sum_{i,j}\sum_{\alpha,\beta}v^{\text{c}}_{ij\alpha\beta}n^{\alpha}_{i}n^{\beta}_{j}\Bigr{]}$
(8)

and

$P_{\text{J}_{\text{s}}}=\exp\Bigl{[}-\sum_{i,j}\sum_{\alpha,\beta}v^{\text{s}}_{ij\alpha\beta}s^{z}_{i\alpha}s^{z}_{j\beta}\Bigr{]}$
(9)

are charge and spin Jastrow factors, which control long-range charge and spin
correlations, respectively. $\alpha$ and $\beta$ run over four orbitals.
$s^{z}_{i\alpha}$ is the $z$ component of the spin angular momentum operator
at site $i$ with orbital $\alpha$. We set
$v^{\text{c}}_{ii\alpha\beta}=v^{\text{s}}_{ii\alpha\beta}=0$ for
$\alpha,\beta=d_{1},d_{2}$ because the on-site correlation between the
$d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals are already taken into account in
$P^{(2)}_{\text{G}}$. The variational parameters in $\left|\Psi\right>$ are
therefore $\\{\tilde{t}^{\alpha\beta}_{ij}\\}$, $\\{\rho^{\alpha}_{i}\\}$,
$\\{m^{\alpha}_{i}\\}$, $\\{\Delta^{\alpha\beta}_{ij}\\}$, $\\{g_{\gamma}\\}$,
$\\{v^{\text{c}}_{ij\alpha\beta}\\}$, and
$\\{v^{\text{s}}_{ij\alpha\beta}\\}$. They are simultaneously optimized using
stochastic reconfiguration method [47]. We show results for
$N_{\text{S}}$=$L\times L$=24$\times$24=576 unit cells (and thus
576$\times$4=2304 orbitals in total), which is large enough to avoid finite
size effects. The antiperiodic boundary conditions are set for both $x$ and
$y$ directions of the primitive lattice vectors.

###  II.3 Band structures

As typical examples of single-layer hole-doped cuprates, we study the La2CuO4
and HgBa2CuO4 systems. We construct maximally localized Wannier orbitals [48,
49] from the first-principles calculation in the local-density approximation
(LDA) with ecalj package [50] and fit them with the hopping integrals $t_{i}$
($i=1-6$) and the site energy of each orbital $\varepsilon_{\alpha}$. The
parameter sets determined for these systems are listed in Table 1. The
explicit form of the tight-binding model is described in Appendix A.1. Note
that $\varepsilon_{d_{z^{2}}}$ for the La-based system is corrected to a lower
value with reference to the quasiparticle self-consistent $GW$ (QSGW) method
[53, 54, 55, 56], which gives realistic correction to the LDA. We have checked
its validity in our previous study [38].

Figure 2 shows the noninteracting tight-binding energy bands for La- and Hg-
based systems. We can summarize the difference between them as follows. (i)
The density of states (DOS) of the $d_{z^{2}}$ component is extended from -1
to -2 eV in the La-based system [Fig. 2(a)], while it is almost localized
around -2 eV in the Hg-based system [Fig. 2(b)]. In addition, a small but
finite peak structure of the $d_{z^{2}}$ component exists around the Fermi
energy (0 eV) in the La-based system, which greatly affects the stability of
superconductivity. These differences are owing to the fact that the
$d_{z^{2}}$ orbital is hybridized with the $p_{x/y}$ orbital in the La-based
system much more strongly than in the Hg-based system. The $d_{z^{2}}$
electrons obtain the itinerancy through the hybridization and therefore the
$d_{z^{2}}$ component of the La-based system is more extended than that of the
Hg-based system. (ii) The site-energy difference between the $d_{x^{2}-y^{2}}$
and $p_{x/y}$ orbitals,
$\Delta_{dp}=\varepsilon_{d_{x^{2}-y^{2}}}-\varepsilon_{p}$, is larger in the
La-based system (2.26 eV) than in the Hg-based system (1.84 eV). This affects
the occupancy of each orbital and thus the strength of the electron
correlation. When $\Delta_{dp}$ is small, the energy band crossing the Fermi
energy contains more weight of the $p_{x/y}$ orbital, which leads to a weaker
electron correlation.

Table 1:  The tight-binding parameters for the La2CuO4 and HgBa2CuO4 systems in eV units estimated on the basis of maximally localized Wannier orbitals from the first-principles LDA and QSGW calculations. The definitions of $t_{i}$ and $\varepsilon_{\alpha}$ are described in Appendix A.1.  | $t_{1}$ | $t_{2}$ | $t_{3}$ | $t_{4}$ | $t_{5}$ | $t_{6}$ | $\varepsilon_{d_{x2-y^{2}}}$ | $\varepsilon_{d_{z^{2}}}$ | $\varepsilon_{p_{x/y}}$
---|---|---|---|---|---|---|---|---|---
La | 1.42 | 0.61 | 0.07 | 0.51 | 0.03 | 0.07 | -0.87 | -0.68 | -3.13
Hg | 1.26 | 0.65 | 0.13 | 0.33 | 0.00 | 0.05 | -1.41 | -1.68 | -3.25
Figure 2:  The energy dispersions and projected density of states onto the
$d_{x^{2}-y^{2}}$, $d_{z^{2}}$, and $p_{x/y}$ orbitals for the noninteracting
tight-binding models for (a) La2CuO4, and (b) HgBa2CuO4 systems. Fermi energy
(defined as zero energy) is set to be the case of 15% hole doping ($x$=0.15).
The high symmetric momenta are indicated as $\Gamma=(0,0)$,
$\text{M}=(\pi,\pi)$, and $\text{X}=(\pi,0)$.

Starting from these energy band structures, we will investigate the ground
state property of the La- and Hg-based systems using the VMC method. We assume
that the tight-binding parameters remain unchanged with hole doping. In the
following, we set $t_{1}$ as a unit of energy. The $d$-orbital on-site Coulomb
interaction $U_{d}/t_{1}$ is varied as a parameter. The other Coulomb
interaction parameters are set as
$(U^{\prime}_{d},J,U_{p},V_{dp},V_{pp})=(0.8,0.1,0.5,0.25,0.2)\,U_{d}$ for
both La- and Hg-based systems with reference to Ref. [57], unless otherwise
noted. Note that the values for the La-based system are larger in eV units
because of the larger $t_{1}$. We first set $V_{dd}=0$ and then discuss the
effect of finite $V_{dd}$ in Sec. III.1.3.

##  III RESULTS

###  III.1 AF and Stripe phases

In this paper, the stripe phase with charge and spin periodicity
$\lambda_{\mathrm{c}}$ and $\lambda_{\mathrm{s}}$ is denoted as a
C$\lambda_{\mathrm{c}}$S$\lambda_{\mathrm{s}}$ phase according to the
convention. For example, the most familiar stripe phase observed around
$x=1/8$ in several La-based systems corresponds to the C4S8 phase. Our main
finding in the charge/spin stripe structure is the “role-sharing?f?f of each
orbital; charge modulation occurs mainly on the $p_{x}$ and $p_{y}$ orbitals,
while spin modulation on the $d_{x^{2}-y^{2}}$ orbital. This is because there
is an efficient energy gain due to the orbital degrees of freedom, which is an
aspect not present in the one-band model. Here we show the energetic
competition between AF and stripe phases, and discuss the changes in stripe
structures along with the doping and material dependence.

Figure 3:  Ground state phase diagrams for the La-based system for the hole
doping rate $x$ with $U_{d}/t_{1}=6$, 7, and 8. AFI, AFM, and PS represent AF
insulator, AF metal, and phase separation, respectively.
C$\lambda_{\mathrm{c}}$S$\lambda_{\mathrm{s}}$ represents the stripe phase
with charge and spin periodicity $\lambda_{\mathrm{c}}$ and
$\lambda_{\mathrm{s}}$, respectively.  Figure 4:  (a) The $x$ dependence of
variational energies per site $E(x)/N_{\mathrm{S}}$ for different phases (AFI,
AFM, and stripe phases with charge periodicity $\lambda_{\mathrm{c}}$) for the
La-based system with $U_{d}/t_{1}=8$.
$f(x)=[-E(0)+8x\\{E(0)-E(1/8)\\}]/N_{\mathrm{S}}$ is added for visibility. The
dashed line represents the Maxwell construction, which indicates the PS for
$0<x<1/8$ (gray shaded area). (b) The excess energy per added hole
$e_{\mathrm{h}}$ for stripe phase as a function of stripe filling $\nu$. (c)
The spin incommensurability $\delta$ as a function of $x$.

####  III.1.1 La2CuO4

Figure 3 depicts the ground-state phase diagram obtained in the La-based
system as a function of the hole doping rate $x$. The AF insulator at undoped
$x=0$ readily becomes unstable with doping. For large $U_{d}/t_{1}$, the AF
metallic phase does not appear. For $0<x<1/8$, the PS occurs between the AFI
and the C4S8 stripe phase. Near $x\sim 1/8$, the C4S8 stripe phase is quite
stable. The structure of the C4S8 phase is consistent with the experimentally
observed $x=1/8$ stripe structure [58]. The stripe phases are widely observed
for $x>1/8$, although the period of the stripes $\lambda_{\mathrm{c}}$
decreases monotonically with further doping. In the obtained C3S3 and C2S4
stripe phases, the ground-state energies increase drastically. Unlike the C4S8
phase, these stripe phases have no clear reason for their stabilization and
also have never been observed experimentally. We believe that these are
artificial states of our VMC calculations, and are in fact liquid-like states
due to the thermal/quantum fluctuation. Let us discuss these stripe structure
in more detail below.

Figure 4(a) shows ground-state energies of AF and stripe phases as a function
of $x$. One can see that AFI and C4S8 phases are fairly stable. The period of
stripe structure $\lambda_{\mathrm{c}}$ monotonically decreases with $x$.
Although the value of $\lambda_{\mathrm{c}}$ is limited to an integer value
due to commensurability and finite size effects, it will vary smoothly in the
limit of infinite system size. For $x<1/8$, the PS between AFI and C4S8 phases
occurs according to the Maxwell construction, and thus the stripe phase for
$\lambda_{\mathrm{c}}>4$ is a metastable state. The presence of the PS has
been reported in previous studies [59, 16, 60, 61]. In real materials, the
effect of lattice distortion is also important for the stability of the stripe
phase. In the La-based systems, such as La2-xSrxCuO4 (LSCO), La2-xBaxCuO4
(LBCO), and La2-x-yNdySrxCuO4 (LNSCO), the CuO2 square lattice undergoes a
structural transition from a high-temperature tetragonal (HTT) to a low-
temperature orthorombic (LTO) phase. The LTO phase is further deformed to a
low-temperature tetragonal (LTT) phase in LBCO and LNSCO. Such structural
distortion couples to the stripe phases and stabilizes them. It is possible
that some metastable stripe structures in the PS become stable through the
structural distortion.

Next, to dissect the stripe structure, we compute the excess energy per added
hole [12], $e_{h}=\frac{E(x)-E(0)}{x},$ and illustrate it in Fig. 4(b) as a
function of the stripe filling $\nu=x\lambda_{\mathrm{c}}$. Here, $E(x)$ is
the VMC total energy at the hole doping rate $x$. The stripe filling $\nu$
denotes the hole filling of the folded band structure in the stripe phase. For
instance, $\nu=0.5$ and $\nu=1$ correspond to a half-filled metallic and a
fully-filled insulating states, respectively. In Fig. 4(b), one can see that
the behavior of $e_{h}$ is drastically changed around
$\lambda_{\mathrm{c}}=4$. For $\lambda_{\mathrm{c}}>4$, $\nu=0.5$ is always
stable, especially, the $\lambda_{\mathrm{c}}=4$ (C4S8) phase observed
experimentally is the most stable. This can be the result of maintaining the
undoped AFI state as much as possible due to the strong AF correlations, as
mentioned in the next paragraph. Note that in our case, the C8S16 phase with
$\nu=1$ proposed in the previous studies does not appear at $x=1/8$. This
phase appears only in peculiar band structures with nearly zero diagonal
(next-nearest-neighbor) hopping [12, 61, 62, 63]. It is thus not so realistic.
On the other hand, for $\lambda_{\mathrm{c}}<4$, the stable value of $\nu$
shifts to smaller values as $\lambda_{\mathrm{c}}$ decreases. Unlike the
stripe structure with $\nu=0.5$, there is no clear reason for such a stripe
structure to be stabilized. In fact, stable stripe phases with
$\lambda_{\mathrm{c}}<4$ have never been observed experimentally for $x>1/8$.
It is natural to assume that such a stripe state is a fluctuating state like a
liquid for $x>1/8$, rather than in a long-range order [64, 65]. More improved
calculations incorporating quantum fluctuations can lead to such liquid states
[26]. It is important that such short-range correlations develop for $x>1/8$
[66], which can have significant effects for the mechanism of
high-$T_{\mathrm{c}}$ superconductivity.

Next, let us discuss the spin degree of freedom. For all the stripe phases
considered here, we find the so-called “spin-charge locking” [58, 67]. Namely,
the relation $2\lambda_{\mathrm{c}}=\lambda_{\mathrm{s}}$
($\lambda_{\mathrm{c}}=\lambda_{\mathrm{s}}$) for even (odd)
$\lambda_{\mathrm{c}}$ is always satisfied. This is because the spin
modulation is a driving force of the stripe ordering and the charge modulation
only follows it [68]. Indeed, we have checked that the stripe phase without
spin modulation, namely, a pure charge stripe phase, is not stabilized within
our calculation. Figure 4(c) shows the spin incommensurability $\delta$ as a
function of $x$. $\delta$ is defined as the difference from the AF wave vector
in $k_{x}$ direction, $\textbf{k}/2\pi=(\frac{1}{2}\pm\delta,\frac{1}{2})$,
where k is a peak position of the spin structure factor of the stripe phase.
In the stripe phases for $\lambda_{\mathrm{c}}\geq 4$, $\delta=x$ is satisfied
[blue shaded area for $1/16\leq x\leq 1/8$ in Fig. 4(c)], which is consistent
with “Yamada relation” for the La-based system confirmed in neutron scattering
experiments [69]. This implies that the stripe phases with
$\lambda_{\mathrm{c}}\geq 4$ are realized. We obtain the relation
$\delta=x=1/\lambda_{\mathrm{s}}~{}(1/2\lambda_{\mathrm{s}})$ for even (odd)
$\lambda_{\mathrm{c}}$ from these stripe filling
$\nu=x\lambda_{\mathrm{c}}=0.5$ and the spin-charge locking. These results
suggest that the origin of the stripe phase in the La-based system is not a
band (nesting) effect but a strong correlation effect. Actually, it has been
reported that $\lambda_{\mathrm{c}}$ expected from the Fermi-surface nesting
has opposite $x$ dependence ($\lambda_{\mathrm{c}}$ increases with increasing
$x$) [70].

Figure 5:  Ground state phase diagrams for the Hg-based system with
$U_{d}/t_{1}=6$, 7, and 8. Notation is the same as in Fig. 3.  Figure 6:  (a)
The $x$ dependence of variational energies per site $E(x)/N_{\mathrm{S}}$ for
different phases (AFI, AFM, and stripe phases with charge periodicity
$\lambda_{\mathrm{c}}$) for the Hg-based system with $U_{d}/t_{1}=6$.
$f(x)=[-E(0)+8x\\{E(0)-E(1/8)\\}]/N_{\mathrm{S}}$ is added for visibility. The
dashed line represents the Maxwell construction, which indicates the PS for
$0.06<x<1/8$ (gray shaded area). (b) The excess energy per added hole
$e_{\mathrm{h}}$ for stripe phase as a function of stripe filling $\nu$. (c)
The spin incommensurability $\delta$ as a function of $x$.

####  III.1.2 HgBa2CuO4

Figure 5 shows the ground state phase diagram of the Hg-based system. Unlike
the La-based system, the AFM phase appears in the low-doping region and the PS
region is narrower. As shown in Fig. 6(a), the energy of several metastable
stripe phases in the PS region is comparable. Although the excess energy per
added hole, $e_{h}$, for $\lambda_{\mathrm{c}}>4$ has the minimum at
$\nu=0.5$, the energy curve is somewhat shallow compared to the La-based
system, as shown in Fig. 6(b). Therefore, the half-filled stripe phases are
not so robust as in the case of La-based system. This implies that the
relation $\delta=x$ for incommensurability [blue shaded area in Fig. 6(c)] is
fragile. In fact, no such relationship has been observed experimentally in the
Hg-based systems. The stripe phases for $x>1/8$ is almost the same as in the
La-based system.

####  III.1.3 Internal structure of stripe phases and the effect of $V_{dd}$

Next, let us consider the internal structure of the stripe phases. Figure 7(a)
shows the obtained C4S8 structure. Reflecting that the cuprates are the
charge-transfer insulators, doped hole carriers are mainly introduced into the
$p$ orbitals, and the occupation number of the $d$ orbitals remains almost
unchanged. Thus, in Fig. 7, only the spin density is depicted for the $d$
orbitals, while only the hole density is depicted for the $p$ orbitals. As can
be seen in Fig. 7(a), the charge-density wave (CDW) of the doped hole carriers
introduced into the $p_{x}$ and $p_{y}$ orbitals appears in-phase. This
structure is actually the same as the $s^{\prime}$-wave bond-order CDW
observed in the La-based systems by the resonant soft X-ray scattering
experiment [71]. For the spin density, the spin-rich Cu sites are surrounded
by the hole-poor O sites, and vice versa. The system efficiently gains the
exchange (kinetic) energy around the hole-poor (hole-rich) sites with this
configuration. In this way, spin-active and charge-active areas alternate in a
stripe pattern as shown in Fig. 7(a). The observation of the $s^{\prime}$-wave
CDW is consistent with a slave-boson mean-field approximation [12] and a
density-matrix-renormalization-group study on a three-band $d$-$p$ model [25].

Figure 7:  Schematic pictures of the obtained hole and spin density
modulations of the C4S8 phase with (a) $s^{\prime}$-wave CDW+SDW and (b)
$d$-wave CDW+SDW. Orange and blue solid circles represent $p_{x}$ and $p_{y}$
holes and each radius is proportional to the hole density. Arrows represent
the $d_{x^{2}-y^{2}}$ spins and each length is proportional to the spin
density.

On the other hand, a stripe structure with an anti-phase CDW of the $p$
orbital in Fig. 7(b) (so-called $d$-wave CDW) has been proposed for non-La-
based systems, such as Bi2Sr2CaCu2O8+δ (Bi2212) [72], Ca2-xNaxCuO2Cl2 (Na-
CCOC) [72], and YBa2Cu3O6+y (YBCO) [73]. In the present calculation, the
$s^{\prime}$-wave CDW is stable, but the $d$-wave CDW appears in cases where
$V_{pp}$ is a little larger, because $V_{pp}$ has a repulsive effect between
neighboring $p_{x}$ and $p_{y}$ holes and favors the anti-phase CDW. We can
also confirm the appearance of the $d$-wave CDW by introducing a finite
$V_{dd}$, although we have so far assumed $V_{dd}=0$. As shown in Figs. 8(a)
and 8(b), the $d$-wave CDW is more likely to appear in the moderately
correlated Hg-based systems than in the strongly correlated La-based systems,
which is consistent with the experimental results.

Figure 8:  Internal structure of the C4S8 phase at $x=56/576\sim 0.097$ for
(a) La-based and (b) Hg-based systems. The order parameter of the SDW in the
C4S8 phase, $M$, for (c) La-based and (d) Hg-based systems. Solid (open)
symbols correspond to the $s^{\prime}$-($d$-)wave CDW+SDW.

In order to clarify the effect of $V_{dd}$ on the stripe, the magnetization in
each state,

$M=\sqrt{\frac{S_{z}(\textbf{q}_{\mathrm{peak}})}{N_{\mathrm{S}}}},$ (10)

is shown in Figs. 8(c) and 8(d), where $S_{z}(\textbf{q}_{\mathrm{peak}})$
represents the peak value of the $z$-component of the $d_{x^{2}-y^{2}}$ spin
structure factor. With the introduction of $V_{dd}$, the average magnetization
decreases monotonically, and the $d$-wave CDW appears for $M\lesssim 0.16$.
Namely, the introduction of $V_{dd}$ works destructively on the SDW, and
eventually the $d$-wave CDW appears. This is because the $V_{dd}$ effectively
weakens $U_{d}$ by increasing the density of doubly-occupied Cu sites [38],
and thus the system becomes less correlated. This situation seems to be a
natural connection to the weak correlation approach. Indeed, it has been
argued in several theoretical proposals that the $d$-wave CDW is induced by
the nesting effects on the Fermi surface [74, 75, 23]. Since our finite-size
calculations have difficulty in capturing the fine structure of the Fermi
surfaces, the approach from weak correlations is complementary.

###  III.2 Superconductivity

Next let us consider $d_{x^{2}-y^{2}}$-wave spin-singlet superconductivity. As
mentioned earlier, we believe that stripe fluctuations are important for the
mechanism of superconductivity, but it is difficult to incorporate this effect
explicitly into the variational functions. Here we calculate the
superconducting correlation function without charge modulation, defined as

$P^{dd}(\textbf{r})=\frac{1}{N_{\text{S}}}\sum_{i}\sum_{\tau,\tau^{\prime}}f^{(dd)}_{\tau\tau^{\prime}}\bigl{<}\Delta^{\dagger}_{\tau}(\textbf{R}_{i})\Delta_{\tau^{\prime}}(\textbf{R}_{i}+\textbf{r})\bigr{>},$
(11)

where $\Delta^{\dagger}_{\tau}(\textbf{R}_{i})$ is a creation operator of
singlet pairs between nearest-neighbor $d_{x^{2}-y^{2}}$ orbitals,

$\Delta^{\dagger}_{\tau}(\textbf{R}_{i})=\frac{1}{\sqrt{2}}(c^{\dagger}_{i1\uparrow}c^{\dagger}_{i+\tau
1\downarrow}+c^{\dagger}_{i+\tau 1\uparrow}c^{\dagger}_{i1\downarrow}),$ (12)

and $\tau$ runs over four nearest-neighbor Cu sites
($\tau=\pm\mathbf{e}_{x},\pm\mathbf{e}_{y}$). $f^{(dd)}_{\tau\tau^{\prime}}$
is a form factor of a superconducting gap function with $d_{x^{2}-y^{2}}$
symmetry, namely, $f^{(dd)}_{\tau\tau^{\prime}}=1$ for
$\tau\parallel\tau^{\prime}$ and $-1$ for $\tau\perp\tau^{\prime}$.
$\langle\cdots\rangle$ denotes
$\langle\Psi|\cdots|\Psi\rangle/\langle\Psi|\Psi\rangle$ for the optimized
variational wave function $|\Psi\rangle$. If $P^{dd}(\textbf{r})$ is saturated
to a finite value for $r=|\mathbf{r}|\rightarrow\infty$, superconducting long-
range order exists. In the following, we take the saturated value of
$P^{dd}(r\rightarrow\infty)$ as a strength of superconductivity $P^{dd}$.

The $x$ dependence of $P^{dd}$ for the La-based system is shown in Fig. 9(a).
At $x=0$, the system is an insulator due to correlation effects, and thus
superconductivity is completely suppressed, $P^{dd}=0$. As $x$ increases,
mobile carriers are introduced into the system and the mobility of the Cooper
pairs increases. On the other hand, the strength of the $d$-$d$ pairing itself
is reduced by doping because the electron correlation is also reduced. The
former effect is predominant in a low-doping, strongly-correlated region and
the latter effect is predominant in a high-doping, weakly-correlated region.
In other words, the peak position of the dome corresponds to a boundary
between strongly- and weakly-correlated regions. Therefore, the peak position
is moved to a smaller value of $x$ with decreasing $U_{d}/t_{1}$ because the
weakly-correlated region is extended.

Figure 9:  The $x$ dependence of superconducting correlation function $P^{dd}$
for (a) La-based and (b) Hg-based systems.

The $x$ dependence of $P^{dd}$ for the Hg-based system is shown in Fig. 9(b).
While the dome-shaped behavior is observed for $U_{d}/t_{1}=7$ and 8, $P^{dd}$
monotonically decreases with $x$ for $U_{d}/t_{1}=6$. This is because for
$U_{d}/t_{1}=6$, the system is not insulating but metallic at $x=0$ when the
paramagnetic state is assumed. When the AF order is taken into account, the
crossover from the Mott insulator to the Slater insulator occurs at
$U_{d}/t_{1}\sim 6.3$ in the present model [38]. Namely, the view of a “doped
Mott insulator” is no longer valid for $U_{d}/t_{1}=6$. The shape of $P^{dd}$
can be a measure of the strength of an electron correlation. From this point
of view, the La-based system is more strongly correlated than the Hg-based
system, because the system is still a Mott insulator at $x=0$ for
$U_{d}/t_{1}=6$, as shown in Fig. 9(a). It results from the larger value of
$\Delta_{dp}$ in the La-based system, which leads to a larger $d$-orbital
occupancy of holes and a stronger electron correlation [38].

We can observe that $P^{dd}$ of the Hg-based system is larger than that of the
La-based system. It is consistent with the higher $T_{\mathrm{c}}$ in the Hg-
based system ($\sim$90 K) compared with the La-based system ($\sim$40 K). The
reason why $P^{dd}$ of the Hg-based system is larger can be understood as
follows: (i) The low $\varepsilon_{d_{z^{2}}}$ in the Hg-based system is
favorable for superconductivity through the less $d_{z^{2}}$-orbital
contribution to the Fermi surface; the $d_{z^{2}}$-orbital contribution to the
Fermi surface is destructive for superconductivity due to the localized
character of the $d_{z^{2}}$ electrons [38]. This effect is significant for a
high-doping ($x\gtrsim 0.15$) region. (ii) The small $\Delta_{dp}(>0)$ in the
Hg-based system is favorable for superconductivity in a low-doping region
($x\lesssim 0.15$) through the weaker electron correlation due to the
hybridization of $p$ orbitals; the electron correlation is destructive for
superconductivity in a low-doping region because the mobility of the Cooper
pairs is suppressed. $P^{dd}$ itself is a physical quantity in the ground
state but is closely related to $T_{\mathrm{c}}$. We consider that both the
band structure and electron correlation effects are necessary to explain the
material dependence of $T_{\mathrm{c}}$.

We also note that the dip structure around $x=1/8$ becomes slightly more
visible as $U_{d}/t_{1}$ increases, in both La-based and Hg-based systems. In
fact, the dip structure around $x=1/8$ is widely observed in cuprate
superconductors. This structure is more pronounced the more stable the stripe
states are, as in LBCO and LNSCO with the LTT distortion mentioned in Sec.
III.1.1. The dip structure we found is tightly related to this fact. This
implies that although our variational wave functions for the superconductivity
do not explicitly include the charge modulation, stripe fluctuations are
partially included through the Jastrow factors.

###  III.3 Heavily-overdoped region

Here, we discuss the hole density of the $d_{z^{2}}$ orbital $n_{z^{2}}$ and
novel magnetism in the heavily-overdoped region. First of all, in Figs. 10(a)
and 10(b), we illustrate $n_{z^{2}}$ as a function of $x$ for several
$J/U_{d}$ at $U_{d}/t_{1}=8$ in the La-based and Hg-based systems, along with
the superconducting correlation function $P_{dd}$. In the La-based system,
$n_{z^{2}}$ increases almost linearly with $x$, and is strongly enhanced by
the Hund’s coupling $J/U_{d}$. Such enhancement of $n_{z^{2}}$ in the heavily-
overdoped region has been also reported by the high-resolution Compton
scattering experiment in LSCO [76]. Note that $n_{z^{2}}$ is suppressed in a
finite $P_{dd}$ region. The negative correlation between $P_{dd}$ and
$n_{z^{2}}$ suggests that the $d_{z^{2}}$ orbital works destructively for
superconductivity, which is consistent with the previous studies [36, 38].
Such superconducting-elastic effect could explain the intriguing pressure
effect where $c$-axis compression reduces superconducting $T_{\mathrm{c}}$
[77], because the $c$-axis compression leads to the increase of $n_{z^{2}}$
through lowering the apical oxygen height. On the other hand, in the Hg-based
system, such an increase of $n_{z^{2}}$ and suppression of $P_{dd}$ are not so
pronounced even at relatively large $J/U_{d}=0.2$. This is because the
$\varepsilon_{d_{z^{2}}}$ is much lower in the Hg-based system and the
$d_{z^{2}}$ orbital is almost inactive.

Figure 10:  The $x$ dependence of superconducting correlation function
$P^{dd}$ (dotted lines) and the average hole density of the $d_{z^{2}}$
orbital $n_{z^{2}}$ (solid lines with symbols) for (a) La-based and (b) Hg-
based systems. The $x$ dependence of local moment $m_{\mathrm{local}}$ for (c)
La-based and (d) Hg-based systems. Solid (open) symbols correspond to the
superconducting (paramagnetic) phase. The colors indicate several $J/U_{d}$
values.

Second, in Figs. 10(c) and 10(d), we show the local moment
$m_{\mathrm{local}}$,

$m_{\mathrm{local}}=\sqrt{\left<S^{z}_{i}S^{z}_{i}\right>}$ (13)

where $S^{z}_{i}=\sum_{\alpha}S^{z}_{i\alpha}$ is the $z$ component of the
total spin angular momentum operator at site $i$ and $\alpha$ runs over four
orbitals in the unit cell. As shown in Figs. 10(c) and 10(d),
$m^{2}_{\mathrm{local}}$ exhibits quite similar behavior to $n_{z^{2}}$. In
the La-based system, it suggests that the spin of doped holes in the
$d_{z^{2}}$ orbital is ferromagnetically aligned with that of the
$d_{x^{2}-y^{2}}$ orbital due to the Hund’s coupling, which leads to the
increase of $m_{\mathrm{local}}$. Namely, new magnetic “seeds” appear in the
heavily-overdoped region originated from the $d_{z^{2}}$ orbital degree of
freedom, which is absent in the one-band Hubbard and even in the three-band
$d$-$p$ models. On the other hand, in the Hg-based system, such a multiorbital
effect is not observed because the $d_{z^{2}}$ orbital is almost inactive. The
relative difference in the position of the $d_{z^{2}}$ orbital causes a clear
difference in the doping dependence of the La and Hg systems.

Finally, let us comment a novel magnetism in the heavily-overdoped region. It
has been discussed in various ways. Magnetic susceptibility measurements
suggest the emergence of local paramagnetic moments in LSCO [78] or
ferromagnetic spin fluctuations in (Bi,Pb)2Sr2CuO6+δ (Bi-2201) [79, 80]. A
zero-field $\mu$SR measurement in LSCO also suggests the development of
ferromagnetic spin fluctuations [81]. Theoretical studies show the possibility
of ferromagnetism using a quantum scaling theory [82], first-principles
calculations for supercells of LBCO [83], determinant quantum Monte Carlo
[84], dynamical cluster approximation [85], and fluctuation exchange
approximation [86]. We can consider two mechanisms for the ferromagnetism in
the heavily-overdoped region. One is the band effect. The peak position of the
spin susceptibility at $\textbf{q}=(\pi,\pi)$ in the low-doping region moves
toward $\textbf{q}=(0,0)$ with hole doping. Indeed, the Lifshitz transition
where the topology of the Fermi surface changes from hole-like to electron-
like occurs around $x\sim 0.3$ in the one-band Hubbard model with
$t^{\prime}/t=-0.3$, leading to the enhancement of $\chi_{0}(0,0)$. The other
is the multiorbital effect, i.e., the double-exchange mechanism originated
from the Hund’s coupling between the $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$
orbitals. Our result supports this mechanism. We expect that these two
mechanisms will not contradict but rather cooperatively work toward
ferromagnetism.

##  IV DISCUSSION AND SUMMARY

In this paper, we present a unified explanation of the doping and material
dependence of the ground state in cuprate superconductors using the VMC
analysis of the four-band $d$-$p$ model derived from the first-principles
calculations. First, we found the following results from a comparison of total
energies of the AFI, AFM, and stripe phases. In the Hg-based system, a slight
hole-doping causes the appearance of AFM phase, but in the La-based systems,
the AFM state is less likely to appear due to strong correlations. In the La-
based system, the experimentally observed $s^{\prime}$-wave bond stripe phase
is quite stable, and all stable stripe structures are in the spin-charge
locking state. Overall features of spin incommensurability observed in the
neutron scattering experiments are explained. In the Hg-based system,
introducing a small but finite $V_{dd}$ or slightly increasing $V_{pp}$ can
lead to a transition to the $d$-wave bond stripe structure and even a spin-
charge unlocked state observed by the nesting effect. PS appears in the low-
doping region of both La- and Hg-based systems. All of these behaviors are
consistent with experimental observations, and our results generally capture
their characteristic features well.

Next, we computed the superconducting correlation functions $P^{dd}$, and
clarified that the dome-shaped $P^{dd}$, indicative of the $T_{\mathrm{c}}$
dome as a function of the hole doping rate, is originated from the dual nature
of the Coulomb repulsion, that is, the enhancement of the pairing interaction
and the depairing effect by the decrease of mobile electron pairs.
$\varepsilon_{d_{z^{2}}}$ and $\Delta_{dp}$ are material-dependent key
parameters for high $T_{\mathrm{c}}$. Our result suggests that deep
$\varepsilon_{d_{z^{2}}}$ and small $\Delta_{dp}$ raise the $T_{\mathrm{c}}$.
This is the reason why $T_{\mathrm{c}}$ of the Hg-based system is higher than
that of the La-based system.

Finally, we investigated novel magnetism in heavily-overdoped region. For the
hole doping rate $x>1/8$, both the $d_{z^{2}}$ hole density $n_{z^{2}}$ and
the local moment at Cu sites increases almost linearly with $x$ in the La-
based system. This means that the $e_{g}$ orbital spins are aligned parallel
via the Hund’s coupling. These local spins can be ferromagnetically aligned by
the double-exchange mechanism, which is consistent with previous proposals. On
the other hand, such a multiorbital effect is not observed in the Hg-based
system because the $d_{z^{2}}$ orbital is located much lower and thus almost
inactive.

These findings, which are beyond one-band description, could provide a major
step toward a full explanation of unconventional normal state and
high-$T_{\mathrm{c}}$ in cuprate supercondutors.

###### Acknowledgements.

The authors thank K. Kuroki and T. Adachi for useful discussions. The
computation has been done using the facilities of the Supercomputer Center,
Institute for Solid State Physics, University of Tokyo and the supercomputer
system HOKUSAI in RIKEN. This work was supported by a Grant-in-Aid for
Scientific Research on Innovative Areas “Quantum Liquid Crystals” (KAKENHI
Grant No. JP19H05825) from JSPS of Japan, and also supported by JSPS KAKENHI
(Grant Nos. JP22K03520, JP22K03512, JP21H04446, JP20K03847, JP19K23433,
JP19H01842, and JP18H01183).

*

##  Appendix A Construction of the trial wave function

Here, we describe the details of the trial wave function together with the
noninteracting tight-binding model obtained on the basis of the first-
principles calculations in Sec. II.3. The construction of the trial wave
function is the most important part for the VMC method. Depending on the trial
state, both real and k space representations are used.

###  A.1 Noninteracting energy band

First, we describe how to construct the noninteracting tight-binding energy
band discussed in Sec. II.3. The noninteracting energy band is obtained by
diagonalizing the following one-body Hamiltonian:

$\displaystyle H_{\text{kin}}$
$\displaystyle=\sum_{\textbf{k},\sigma}\left(c_{\textbf{k}1\sigma}^{\dagger},c_{\textbf{k}2\sigma}^{\dagger},c_{\textbf{k}3\sigma}^{\dagger},c_{\textbf{k}4\sigma}^{\dagger}\right)\begin{pmatrix}t_{11}&t^{*}_{21}&t^{*}_{31}&t^{*}_{41}\\\
t_{21}&t_{22}&t^{*}_{32}&t^{*}_{42}\\\ t_{31}&t_{32}&t_{33}&t^{*}_{43}\\\
t_{41}&t_{42}&t_{43}&t_{44}\end{pmatrix}\begin{pmatrix}c_{\textbf{k}1\sigma}\\\
c_{\textbf{k}2\sigma}\\\ c_{\textbf{k}3\sigma}\\\
c_{\textbf{k}4\sigma}\end{pmatrix}$ (14)
$\displaystyle=\sum_{\textbf{k},\sigma}\sum_{m}E_{m}(\textbf{k})a^{\dagger}_{\textbf{k}m\sigma}a_{\textbf{k}m\sigma}$
(15)

with the hopping matrix elements given as

$\displaystyle t_{11}$ $\displaystyle=\varepsilon_{d_{x^{2}-y^{2}}},$ (16)
$\displaystyle t_{21}$ $\displaystyle=0,$ (17) $\displaystyle t_{22}$
$\displaystyle=\varepsilon_{d_{z^{2}}}-2t_{5}(\cos k_{x}+\cos k_{y}),$ (18)
$\displaystyle t_{31}$ $\displaystyle=2\text{i}t_{1}\sin\frac{1}{2}k_{x},$
(19) $\displaystyle t_{32}$
$\displaystyle=-2\text{i}t_{4}\sin\frac{1}{2}k_{x},$ (20) $\displaystyle
t_{33}$ $\displaystyle=\varepsilon_{p_{x}}+2t_{3}\cos
k_{x}+2t_{6}[\cos(k_{x}+k_{y})+\cos(k_{x}-k_{y})],$ (21) $\displaystyle
t_{41}$ $\displaystyle=-2\text{i}t_{1}\sin\frac{1}{2}k_{y},$ (22)
$\displaystyle t_{42}$ $\displaystyle=-2\text{i}t_{4}\sin\frac{1}{2}k_{y},$
(23) $\displaystyle t_{43}$
$\displaystyle=2t_{2}\left[\cos\left(\frac{1}{2}k_{x}+\frac{1}{2}k_{y}\right)-\cos\left(\frac{1}{2}k_{x}-\frac{1}{2}k_{y}\right)\right],$
(24) $\displaystyle t_{44}$ $\displaystyle=\varepsilon_{p_{y}}+2t_{3}\cos
k_{y}+2t_{6}[\cos(k_{x}+k_{y})+\cos(k_{x}-k_{y})],$ (25)

where $c^{\dagger}_{\textbf{k}\alpha\sigma}$ ($c_{\textbf{k}\alpha\sigma}$) is
a creation (annihilation) operator of an electron with momentum k, spin
$\sigma\,(=\uparrow,\downarrow)$, and orbital $\alpha\,(=1,2,3,4)$
corresponding to ($d_{x^{2}-y^{2}}$, $d_{z^{2}}$, $p_{x}$, $p_{y}$),
respectively. The bonds between nearest-neighbor Cu sites are set as unit
vectors ($|\textbf{e}_{x}|=|\textbf{e}_{y}|=1$) and the bonds between nearest-
neighbor Cu and O sites are one-half of them (see Fig. 1). The hopping
integrals $t_{i}$ ($i=1-6$) and the site energy of each orbital
$\varepsilon_{\alpha}$ are determined by fitting the band structures that are
obtained by the LDA or QSGW calculation. The specific values for the La and
Hg-based systems are listed in Table 1. Equation (15) is obtained by
diagonalizing the Hamiltonian matrix in Eq. (14) and is the same with Eq. (3)
in Sec. II.1. $E_{m}(\textbf{k})$ is the noninteracting energy band
characterized by the wave vector k and the energy band index $m\,(=1,2,3,4)$
with $a^{\dagger}_{\textbf{k}m\sigma}$ ($a_{\textbf{k}m\sigma}$) being a
creation (annihilation) operator of the corresponding energy band with spin
$\sigma$.

###  A.2 Trial wave function

####  A.2.1 Superconductivity

To construct the trial wave function for superconductivity, we employ the
Bogoliubov de-Gennes (BdG) type Hamiltonian in real space [87], i.e.,

$H_{\text{BdG}}=\sum_{i,j}\sum_{\alpha,\beta}\left(c^{\dagger}_{i\alpha\uparrow},c_{i\alpha\downarrow}\right)\begin{pmatrix}T^{\alpha\beta}_{ij\uparrow}&\Delta^{\alpha\beta}_{ij}\\\
\Delta^{\alpha\beta*}_{ji}&-T^{\alpha\beta}_{ji\downarrow}\end{pmatrix}\begin{pmatrix}c_{j\beta\uparrow}\\\
c^{\dagger}_{j\beta\downarrow}\end{pmatrix}.$ (26)

Here, $T^{\alpha\beta}_{ij\sigma}$ is obtained from the Fourier transform of
the matrix in Eq. (14) with renormalized hopping integrals $\tilde{t}_{i}$ and
also includes the chemical potential term. The chemical potential $\mu$ is set
to the corresponding Fermi energy. $\Delta^{\alpha\beta}_{ij}$ corresponds to
an anomalous part that describes the superconducting pairing in real space.
Therefore, the variational parameters to be optimized in $\left|\Phi\right>$
are $\tilde{t}_{i}$ ($i=2-6$) and $\\{\Delta^{\alpha\beta}_{ij}\\}$ with
$\tilde{t}_{1}=t_{1}$ being fixed as a unit of energy. In this study, we
consider the pairing between nearest-neighbor orbitals, $d$-$d$, $d$-$p_{x}$,
$d$-$p_{y}$, $p_{x}$-$p_{x}$, and $p_{y}$-$p_{y}$, where $d$ denotes the
$d_{x^{2}-y^{2}}$ orbital. In the paramagnetic phase, we simply set
$\Delta^{\alpha\beta}_{ij}=0$.

We can also construct the trial wave function with the band (k space)
representation, where
$\Delta^{mn}(\textbf{k})\propto\bigl{<}a_{\textbf{k}m\uparrow}a_{-\textbf{k}n\downarrow}\bigr{>}$
is the variational parameter. However, we find that the trial wave function
with the real space representation always gives lower (i.e., better) energy
than that with the band representation, especially, for large Coulomb
interaction parameters. This is because the Coulomb interaction depends on the
orbital, not on the band, and thus the trial wave function with the real space
(orbital) representation gives the better result.

####  A.2.2 Uniform spin AF and stripe phases

As mentioned in Sec. II.2, various long-range orderings of charge and spin can
be described by introducing $\\{\rho^{\alpha}_{i}\\}$ and
$\\{m^{\alpha}_{i}\\}$. A uniform spin AF phase along the $z$\- and
$x$-direction can be introduced as

$m^{\alpha}_{(z)}\sum_{i}e^{i\textbf{Q}\cdot\textbf{r}_{i}}\left(c^{\dagger}_{i\alpha\uparrow}c_{i\alpha\uparrow}-c^{\dagger}_{i\alpha\downarrow}c_{i\alpha\downarrow}\right)$
(27)

and

$m^{\alpha}_{(x)}\sum_{i}e^{i\textbf{Q}\cdot\textbf{r}_{i}}\left(c^{\dagger}_{i\alpha\uparrow}c_{i\alpha\downarrow}+c^{\dagger}_{i\alpha\downarrow}c_{i\alpha\uparrow}\right)$
(28)

for each orbital $\alpha\,(=1,2,3,4)$, where $\textbf{Q}=(\pi,\pi)$, and
$m^{\alpha}_{(z)}$ and $m^{\alpha}_{(x)}$ are treated as variational
parameters. They are energetically degenerate within a one-body description.
However, they will give different variational energies with the Gutzwiller and
Jastrow factors in Eq. (6) that break the SU(2) symmetry. We have confirmed
that $m^{\alpha}_{(x)}$ always gives a lower variational energy and thus a
better trial state. This is because the Gutzwiller and Jastrow factors
generate spin fluctuations in the direction orthogonal to that of
$m^{\alpha}_{(x)}$.

For a stripe phase with charge and spin periodicities
$\lambda^{\alpha}_{\text{c}}=2\pi/q^{\alpha}_{\text{c}}$ and
$\lambda^{\alpha}_{\text{s}}=2\pi/q^{\alpha}_{\text{s}}$, respectively, the
following potentials with spatial modulation in the $x$ direction should be
introduced at site $i$ for each orbital $\alpha\,(=1,2,3,4)$:

$\rho^{\alpha}_{i}=\rho^{\alpha}\cos[q^{\alpha}_{\text{c}}(x_{i}-x^{\alpha}_{\text{c}})]$
(29)

and

$m^{\alpha}_{i}=(-1)^{x_{i}+y_{i}}m^{\alpha}\sin[q^{\alpha}_{\text{s}}(x_{i}-x^{\alpha}_{\text{s}})],$
(30)

where $\rho^{\alpha}$ and $m^{\alpha}$ are the amplitude of charge and spin
orderings, respectively. $x^{\alpha}_{\text{c}}$ and $x^{\alpha}_{\text{s}}$
control the relative phases of charge and spin orderings, respectively.
$\lambda^{\alpha}_{\text{c}}$, $\lambda^{\alpha}_{\text{s }}$,
$\rho^{\alpha}$, $m^{\alpha}$, $x^{\alpha}_{\text{c}}$ and
$x^{\alpha}_{\text{s}}$ are all variational parameters to be optimized. For
example, the stripe phase observed around $x=1/8$ in several cuprate
superconductors corresponds to $\lambda_{\mathrm{c}}=4$ and
$\lambda_{\mathrm{s}}=8$, although the orbital dependence and relative phase
are still under debate.

## References

  * [1] J. G. Bednorz and K. A. Muller, Possible high$T_{\text{c}}$ superconductivity in the Ba-La-Cu-O system, Z. Phys. B 64, 189 (1986).
  * [2] Y. Sato, S. Kasahara, H. Murayama, Y. Kasahara, E.-G. Moon, T. Nishizaki, T. Loew, J. Porras, B. Keimer, T. Shibauchi, and Y. Matsuda, Thermodynamic evidence for a nematic phase transition at the onset of the pseudogap in YBa2Cu3Oy, Nat. Phys. 13, 1074 (2017).
  * [3] S. Murakami, N. Nagaosa, and S.-C. Zhang, Dissipationless Quantum Spin Current at Room Temperature, Science 301, 1348 (2003).
  * [4] Y. K. Kato, R. C. Myers, A. C. Gossard, and D. D. Awschalom, Observation of the Spin Hall Effect in Semiconductors, Science 306, 1910 (2004).
  * [5] K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka, H. Kontani, and H. Aoki, Unconventional Pairing Originating from the Disconnected Fermi Surfaces of Superconducting LaFeAsO1-xFx, Phys. Rev. Lett. 101, 087004 (2008).
  * [6] D. F. Agterberg, P. M. R. Brydon, and C. Timm, Bogoliubov Fermi Surfaces in Superconductors with Broken Time-Reversal Symmetry, Phys. Rev. Lett. 118, 127001 (2017).
  * [7] A. Damascelli, Z. Hussain, and Z.-X. Shen, Angle-resolved photoemission studies of the cuprate superconductors, Rev. Mod. Phys. 75, 473 (2003).
  * [8] V. J. Emery, Theory of high-$T_{c}$ superconductivity in oxides, Phys. Rev. Lett. 58, 2794 (1987).
  * [9] T. Asahata, A. Oguri, and S. Maekawa, Superconducting State in the Three-Band Hubbard Model: A Variational Monte Carlo Study, J. Phys. Soc. Jpn. 65, 365 (1996).
  * [10] T. Takimoto and T. Moriya, Theory of Spin Fluctuation-Induced Superconductivity Based on a $d$-$p$ Model, J. Phys. Soc. Jpn. 66, 2459 (1997); T. Takimoto and T. Moriya, Theory of Spin Fluctuation-Induced Superconductivity Based on a $d$-$p$ Model. II. Superconducting State, ibid. 67, 3570 (1998).
  * [11] T. Yanagisawa, S. Koike, and K. Yamaji, Ground state of the three-band Hubbard model, Phys. Rev. B 64, 184509 (2001).
  * [12] J. Lorenzana and G. Seibold, Metallic Mean-Field Stripes, Incommensurability, and Chemical Potential in Cuprates, Phys. Rev. Lett. 89, 136401 (2002).
  * [13] S. Shinkai, H. Ikeda, and K. Yamada, Study of Superconducting Transition Temperature in d-p Model on Basis of Perturbation Theory, J. Phys. Soc. Jpn. 75, 104712 (2006).
  * [14] P. R. C. Kent, T. Saha-Dasgupta, O. Jepsen, O. K. Andersen, A. Macridin, T. A. Maier, M. Jarrell, and T. C. Schulthess, Combined density functional and dynamical cluster quantum Monte Carlo calculations of the three-band Hubbard model for hole-doped cuprate superconductors, Phys. Rev. B 78, 035132 (2008).
  * [15] R. Thomale and M. Greiter, Numerical analysis of three-band models for CuO planes as candidates for a spontaneous T-violating orbital current phase, Phys. Rev. B 77, 094511 (2008).
  * [16] E. Arrigoni, M. Aichhorn, M. Daghofer, and W. Hanke, Phase diagram and single-particle spectrum of CuO2 high-$T_{c}$ layers: variational cluster approach to the three-band Hubbard model, New J. Phys. 11, 055066 (2009).
  * [17] C. Weber, A. Läuchli, F. Mila, and T. Giamarchi, Orbital Currents in Extended Hubbard Models of High-$T_{c}$ Cuprate Superconductors, Phys. Rev. Lett. 102, 017005 (2009).
  * [18] C. Weber, K. Haule, and G. Kotliar, Strength of correlations in electron- and hole-doped cuprates, Nat. Phys. 6, 574 (2010).
  * [19] C. Weber, K. Haule, and G. Kotliar, Apical oxygens and correlation strength in electron- and hole-doped copper oxides, Phys. Rev. B 82, 125107 (2010).
  * [20] M. H. Fischer and E.-A. Kim, Mean-field analysis of intra-unit-cell order in the Emery model of the CuO2 plane, Phys. Rev. B 84, 144502 (2011).
  * [21] C. Weber, C. Yee, K. Haule, and G. Kotliar, Scaling of the transition temperature of hole-doped cuprate superconductors with the charge-transfer energy, Europhys. Lett. 100, 37001 (2012).
  * [22] S. Bulut, W. A. Atkinson, and A. P. Kampf, Spatially modulated electronic nematicity in the three-band model of cuprate superconductors, Phys. Rev. B 88, 155132 (2013).
  * [23] Y. Yamakawa and H. Kontani, Spin-Fluctuation-Driven Nematic Charge-Density Wave in Cuprate Superconductors: Impact of Aslamazov-Larkin Vertex Corrections, Phys. Rev. Lett. 114, 257001 (2015).
  * [24] D. Ogura and K. Kuroki, Asymmetry of superconductivity in hole- and electron-doped cuprates: Explanation within two-particle self-consistent analysis for the three-band model, Phys. Rev. B 92, 144511 (2015).
  * [25] S. R. White and D. J. Scalapino, Doping asymmetry and striping in a three-orbital CuO2 Hubbard model, Phys. Rev. B 92, 205112 (2015).
  * [26] E. W. Huang, C. B. Mendl, S. Liu, S. Johnston, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Numerical evidence of fluctuating stripes in the normal state of high-$T_{\mathrm{c}}$ cuprate superconductors, Science 358, 1161 (2017).
  * [27] M. Tsuchiizu, K. Kawaguchi, Y. Yamakawa, and H. Kontani, Multistage electronic nematic transitions in cuprate superconductors: A functional-renormalization-group analysis, Phys. Rev. B 97, 165131 (2018).
  * [28] P. P. Orth, B. Jeevanesan, R. M. Fernandes, and J. Schmalian, Enhanced nematic fluctuations near an antiferromagnetic Mott insulator and possible application to high-$T_{c}$ cuprates, npj Quantum Mater. 4, 4 (2019).
  * [29] M. Zegrodnik, A. Biborski, M. Fidrysiak, and J. Spałek, Superconductivity in the three-band model of cuprates: Variational wave function study and relation to the single-band case, Phys. Rev. B 99, 104511 (2019).
  * [30] S. S. Dash and D. Sénéchal, Pseudogap transition within the superconducting phase in the three-band Hubbard model, Phys. Rev. B 100, 214509 (2019).
  * [31] A. Moreo and E. Dagotto, Minimal-size real-space $d$-wave pairing operator in CuO2 planes, Phys. Rev. B 100, 214502 (2019).
  * [32] A. Biborski, M. Zegrodnik, and J. Spałek, Superconducting properties of the hole-doped three-band $d-p$ model studied with minimal-size real-space $d$-wave pairing operators, Phys. Rev. B 101, 214504 (2020).
  * [33] Z.-H. Cui, C. Sun, U. Ray, B.-X. Zheng, Q. Sun, and G. K.-L. Chan, Ground-state phase diagram of the three-band Hubbard model from density matrix embedding theory, Phys. Rev. Research 2, 043259 (2020).
  * [34] A. Chiciak, E. Vitali, and S. Zhang, Magnetic and charge orders in the ground state of the Emery model: Accurate numerical results, Phys. Rev. B 102, 214512 (2020).
  * [35] P. Mai, G. Balduzzi, S. Johnston, and T. Maier, Orbital structure of the effective pairing interaction in the high-temperature superconducting cuprates, npj Quantum Mater. 6, 26 (2021).
  * [36] H. Sakakibara, H. Usui, K. Kuroki, R. Arita, and H. Aoki, Two-Orbital Model Explains the Higher Transition Temperature of the Single-Layer Hg-Cuprate Superconductor Compared to That of the La-Cuprate Superconductor, Phys. Rev. Lett. 105, 057003 (2010).
  * [37] C. E. Matt, D. Sutter, A. M. Cook, Y. Sassa, M. Månsson, O. Tjernberg, L. Das, M. Horio, D. Destraz, C. G. Fatuzzo, K. Hauser, M. Shi, M. Kobayashi, V. N. Strocov, T. Schmitt, P. Dudin, M. Hoesch, S. Pyon, T. Takayama, H. Takagi, O. J. Lipscombe, S. M. Hayden, T. Kurosawa, N. Momono, M. Oda, T. Neupert, and J. Chang, Direct observation of orbital hybridisation in a cuprate superconductor, Nat. Commun. 9, 972 (2018).
  * [38] H. Watanabe, T. Shirakawa, K. Seki, H. Sakakibara, T. Kotani, H. Ikeda, and S. Yunoki, Unified description of cuprate superconductors using a four-band $d$-$p$ model, Phys. Rev. Research 3, 033157 (2021).
  * [39] Y. Ohta, T. Tohyama, and S. Maekawa, Apex oxygen and critical temperature in copper oxide superconductors: Universal correlation with the stability of local singlets, Phys. Rev. B 43, 2968 (1991).
  * [40] J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Prog. Theor. Phys. 30, 275 (1963).
  * [41] P. Hansmann, N. Parragh, A. Toschi, G. Sangiovanni, and K. Held, Importance of $d-p$ Coulomb interaction for high Tc cuprates and other oxides, New J. Phys. 16, 033009 (2014).
  * [42] W. L. McMillan, Ground State of Liquid He4, Phys. Rev. 138, A442 (1965).
  * [43] D. Ceperley, G. V. Chester, and M. H. Kalos, Monte Carlo simulation of a many-fermion study, Phys. Rev. B 16, 3081 (1977).
  * [44] H. Yokoyama and H. Shiba, Variational Monte-Carlo Studies of Hubbard Model. I, J. Phys. Soc. Jpn. 56, 1490 (1987).
  * [45] J. Bünemann, W. Weber, and F. Gebhard, Multiband Gutzwiller wave functions for general on-site interactions, Phys. Rev. B 57, 6896 (1998).
  * [46] H. Watanabe, K. Seki, and S. Yunoki, Charge-density wave induced by combined electron-electron and electron-phonon interactions in 1$T$-TiSe2: A variational Monte Carlo study, Phys. Rev. B 91, 205135 (2015).
  * [47] S. Sorella, Generalized Lanczos algorithm for variational quantum Monte Carlo, Phys. Rev. B 64, 024512 (2001); S. Yunoki and S. Sorella, Two spin liquid phases in the spatially anisotropic triangular Heisenberg model, ibid. 74, 014408 (2006).
  * [48] N. Marzari and D. Vanderbilt, Maximally localized generalized Wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
  * [49] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized Wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001).
  * [50] A first-principles electronic-structure suite based on the plane-wave + muffin-tin orbital (PMT) method [51, 52], ecalj package, is freely available from https://github.com/tkotani/ecalj.
  * [51] T. Kotani and M. van Schilfgaarde, Fusion of the LAPW and LMTO methods: The augmented plane wave plus muffin-tin orbital method, Phys. Rev. B 81, 125117 (2010).
  * [52] T. Kotani, H. Kino, and H. Akai, Formulation of the Augmented Plane-Wave and Muffin-Tin Orbital Method, J. Phys. Soc. Jpn. 84, 034702 (2015).
  * [53] S. V. Faleev, M. van Schilfgaarde, and T. Kotani, All-Electron Self-Consistent $GW$ Approximation: Application to Si, MnO, and NiO, Phys. Rev. Lett. 93, 126406 (2004).
  * [54] M. van Schilfgaarde, T. Kotani, and S. V. Faleev, Quasiparticle Self-Consistent $GW$ Theory, Phys. Rev. Lett. 96, 226402 (2006).
  * [55] T. Kotani, M. van Schilfgaarde, and S. V. Faleev, Quasiparticle self-consistent $GW$ method: A basis for the independent-particle approximation, Phys. Rev. B 76, 165106 (2007).
  * [56] S. W. Jang, T. Kotani, H. Kino, K. Kuroki, and M. J. Han, Quasiparticle self-consistent $GW$ study of cuprates: electronic structure, model parameters and the two-band theory for Tc, Sci. Rep. 5, 12050 (2015).
  * [57] M. Hirayama, T. Misawa, T. Ohgoe, Y. Yamaji, and M. Imada, Effective Hamiltonian for cuprate superconductors derived from multiscale ab initio scheme with level renormalization, Phys. Rev. B 99, 245155 (2019).
  * [58] J. M. Tranquada, B. J. Sternlleb, J. D. Axe, Y. Nakamura, and S. Uchida, Evidence for stripe correlations of spins and holes in copper oxide superconductors, Nature 375, 561 (1995).
  * [59] V. J. Emery and S. A. Kivelson, Frustrated electronic phase separation and high-temperature superconductors, Physica C 209, 597 (1993).
  * [60] T. Misawa and M. Imada, Origin of high-$T_{c}$ superconductivity in doped Hubbard models and their extensions: Roles of uniform charge fluctuations, Phys. Rev. B 90, 115137 (2014).
  * [61] K, Ido, T, Ohgoe, and M. Imada, Competition among various charge-inhomogeneous states and $d$-wave superconducting state in Hubbard models on square lattices, Phys. Rev. B 97, 045138 (2018).
  * [62] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Stripe order in the underdoped region of the two-dimensional Hubbard model, Science 358, 1155 (2017).
  * [63] H.-C. Jiang and T. P. Devereaux, Superconductivity in the doped Hubbard model and its interplay with next-nearest hopping $t^{\prime}$, Science 365, 1424 (2019).
  * [64] J.-J. Wen, H. Huang, S.-J. Lee, H. Jang, J. Knight, Y. S. Lee, M. Fujita, K. M. Suzuki, S. Asano, S. A. Kivelson, C.-C. Kao, and J.-S. Lee, Observation of two types of charge-density-wave orders in superconducting La2-xSrxCuO4, Nat. Commun. 10, 3269 (2019).
  * [65] K. von Arx, Q. Wang, S. Mustafi, D. G. Mazzone, M. Horio, D. John Mukkattukavil, E. Pomjakushina, S. Pyon, T. Takayama, H. Takagi, T. Kurosawa, N. Momono, M. Oda, N. B. Brookes, D. Betto, W. Zhang, T. C. Asmara, Y. Tseng, T. Schmitt, Y. Sassa, and J. Chang, Fate of charge order in overdoped La-based cuprates, arXiv: 2206.06695 (2022).
  * [66] G. Seibold, R. Arpaia, Y. Y. Peng, R. Fumagalli, L. Braicovich, C. Di Castro, M. Grilli, G. Ghiringhelli, and S. Caprara, Strange metal behaviour from charge density fluctuations in cuprates, Commun. Phys. 4, 7 (2021).
  * [67] E. Blackburn, J. Chang, M. Hücker, A. T. Holmes, N. B. Christensen, Ruixing Liang, D. A. Bonn, W. N. Hardy, U. Ruütt, O. Gutowski, M. v. Zimmermann, E. M. Forgan, and S. M. Hayden, X-Ray Diffraction Observations of a Charge-Density-Wave Order in Superconducting Ortho-II YBa2Cu3O6.54 Single Crystals in Zero Magnetic Field, Phys. Rev. Lett. 110, 137004 (2013).
  * [68] L. F. Tocchio, A. Montorsi, and F. Becca, Metallic and insulating stripes and their relation with superconductivity in the doped Hubbard model, SciPost Phys. 7, 021 (2019).
  * [69] K. Yamada, C. H. Lee, K. Kurahashi, J. Wada, S. Wakimoto, S. Ueki, H. Kimura, Y. Endoh, S. Hosoya, G. Shirane, R. J. Birgeneau, M. Greven, M. A. Kastner, and Y. J. Kim, Doping dependence of the spatially modulated dynamical spin correlations and the superconducting-transition temperature in La2-xSrxCuO4, Phys. Rev. B 57, 6165 (1998).
  * [70] H. Miao, R. Fumagalli, M. Rossi, J. Lorenzana, G. Seibold, F. Yakhou-Harris, K. Kummer, N. B. Brookes, G. D. Gu, L. Braicovich, G. Ghiringhelli, and M. P. M. Dean, Formation of Incommensurate Charge Density Waves in Cuprates, Phys. Rev. X 9, 031042 (2019).
  * [71] A. J. Achkar, F. He, R. Sutarto, Christopher McMahon, M. Zwiebler, M. Hücker, G. D. Gu, Ruixing Liang, D. A. Bonn, W. N. Hardy, J. Geck, and D. G. Hawthorn, Orbital symmetry of charge-density-wave order in La1.875Ba0.125CuO4 and YBa2Cu3O6.67, Nature Mater. 15, 616 (2016).
  * [72] K. Fujita, M. H. Hamidian, S. D. Edkins, C. K. Kim, Y. Kohsaka, M. Azuma, M. Takano, H. Takagi, H. Eisaki, S. Uchida, A. Allais, M. J. Lawler, E.-A. Kim, S. Sachdev, and J. C. S. Davis, Direct phase-sensitive identification of a $d$-form factor density wave in underdoped cuprates, Proc. Natl. Acad. Sci. USA 111, E3026 (2014).
  * [73] R. Comin, R. Sutarto, F. He, E. H. da Silva Neto, L. Chauviere, A. Fraño, R. Liang, W. N. Hardy, D. A. Bonn, Y. Yoshida, H. Eisaki, A. J. Achkar, D. G. Hawthorn, B. Keimer, G. A. Sawatzky, and A. Damascelli, Symmetry of charge order in cuprates, Nature Mater. 14, 796 (2015).
  * [74] S. Sachdev and R. La Placa, Bond Order in Two-Dimensional Metals with Antiferromagnetic Exchange Interactions, Phys. Rev. Lett. 111, 027202 (2013).
  * [75] K. B. Efetov, H. Meier, and C. Pépin, Pseudogap state near a quantum critical point, Nature Phys. 9, 442 (2013).
  * [76] Y. Sakurai, M. Itou, B. Barbiellini, P. E. Mijnarends, R. S. Markiewicz, S. Kaprzyk, J.-M. Gillet, S. Wakimoto, M. Fujita, S. Basak, Yung Jui Wang, W. Al-Sawai, H. Lin, A. Bansil, K. Yamada, Imaging Doped Holes in a Cuprate Superconductor with High-Resolution Compton Scattering, Science 332, 698 (2011).
  * [77] F. Hardy, N. J. Hillier, C. Meingast, D. Colson, Y. Li, N. Barišić, G. Yu, X. Zhao, M. Greven, and J. S. Schilling, Enhancement of the Critical Temperature of HgBa2CuO4+δ by Applying Uniaxial and Hydrostatic Pressure: Implications for a Universal Trend in Cuprate Superconductors, Phys. Rev. Lett. 105, 167002 (2010).
  * [78] S. Wakimoto, R. J. Birgeneau, A. Kagedan, Hyunkyung Kim, I. Swainson, K. Yamada, and H. Zhang, Magnetic properties of the overdoped superconductor La2-xSrxCuO4 with and without Zn impurities, Phys. Rev. B 72, 064521 (2005).
  * [79] K. Kurashima, T. Adachi, K. M. Suzuki, Y. Fukunaga, T. Kawamata, T. Noji, H. Miyasaka, I. Watanabe, M. Miyazaki, A. Koda, R. Kadono, and Y. Koike, Development of Ferromagnetic Fluctuations in Heavily Overdoped (Bi,Pb)2Sr2CuO6+δ Copper Oxides, Phys. Rev. Lett. 121, 057002 (2018).
  * [80] Y. Komiyama, S. Onishi, M. Harada, H. Kuwahara, H. Kuroe, K. Kurashima, T. Kawamata, Y. Koike, I. Watanabe, and T. Adachi, Magnetic Impurity Effects on Ferromagnetic Fluctuations in Heavily Overdoped (Bi,Pb)2Sr2Cu1-yFeyO6+δ Cuprates, J. Phys. Soc. Jpn, 90, 084701 (2021).
  * [81] J. E. Sonier, C. V. Kaiser, V. Pacradouni, S. A. Sabok-Sayr, C. Cochrane, D. E. MacLaughlin, S. Komiya, and N. E. Hussey, Direct search for a ferromagnetic phase in a heavily overdoped nonsuperconducting copper oxide, Proc. Natl. Acad. Sci. U.S.A. 107, 17131 (2010).
  * [82] A. Kopp, A. Ghosal, and S. Chakravarty, Competing ferromagnetism in high-temperature copper oxide superconductors, Proc. Natl. Acad. Sci. USA 104, 6123 (2007).
  * [83] B. Barbiellini and T. Jarlborg, Importance of Local Band Effects for Ferromagnetism in Hole-Doped La2CuO4 Cuprate Superconductors, Phys. Rev. Lett. 101, 157002 (2008).
  * [84] C. J. Jia, E. A. Nowadnick, K. Wohlfeld, Y. F. Kung, C.-C. Chen, S. Johnston, T. Tohyama, B. Moritz, and T. P. Devereaux, Persistent spin excitations in doped antiferromagnets revealed by resonant inelastic light scattering, Nat. Commun. 5, 3314 (2014).
  * [85] T. A. Maier and D. J. Scalapino, Disappearance of Superconductivity in the Overdoped Cuprates, J. Supercond. Novel Magn. 33, 15 (2020).
  * [86] S. Teranishi, K. Nishiguchi, S. Yunoki, and K. Kusakabe, Effect of On-site Coulomb Repulsion on Ferromagnetic Fluctuations in Heavily Overdoped Cuprates, J. Phys. Soc. Jpn. 90, 094707 (2021).
  * [87] A. Himeda, T. Kato, and M. Ogata, Stripe States with Spatially Oscillating $d$-Wave Superconductivity in the Two-Dimensional $t$-$t^{\prime}$-$J$ Model, Phys. Rev. Lett. 88, 117001 (2002).


# Learning Globally Optimized Language Structure via Adversarial Training

Xuwang Yin
<EMAIL_ADDRESS>
###### Abstract

Recent work has explored integrating autoregressive language models with
energy-based models (EBMs) to enhance text generation capabilities. However,
learning effective EBMs for text is challenged by the discrete nature of
language. This work proposes an adversarial training strategy to address
limitations in prior efforts. Specifically, an iterative adversarial attack
algorithm is presented to generate negative samples for training the EBM by
perturbing text from the autoregressive model. This aims to enable the EBM to
suppress spurious modes outside the support of the data distribution.
Experiments on an arithmetic sequence generation task demonstrate that the
proposed adversarial training approach can substantially enhance the quality
of generated sequences compared to prior methods. The results highlight the
promise of adversarial techniques to improve discrete EBM training. Key
contributions include: (1) an adversarial attack strategy tailored to text to
generate negative samples, circumventing MCMC limitations; (2) an adversarial
training algorithm for EBMs leveraging these attacks; (3) empirical validation
of performance improvements on a sequence generation task.

##  1 Introduction

The recent advancement in natural language generation (NLG) has been
predominantly driven by autoregressive models, particularly Locally Normalized
Language Models (LNLMs). These models, including variations of Recurrent
Neural Networks (RNNs) and the Transformer architecture, have shown remarkable
capabilities in generating coherent and contextually relevant text. However,
they often suffer from limitations inherent to their local scope — optimizing
the likelihood of each token solely based on its preceding context. This local
perspective can result in sequences that, although locally plausible, are
globally inconsistent or suboptimal, leading to limitations in capturing long-
range dependencies and semantic coherence.

To address these shortcomings, researchers have explored the integration of
LNLMs with Energy-Based Models (EBMs), which offer a more global view of
sequence generation by leveraging an energy function parameterized across the
entire sequence (Wang and Ou, 2018; Parshakova et al., 2019; Deng et al.,
2019).

However, the application of EBMs to text-based tasks is challenged by the
discrete nature of text, which hinders the use of gradient-based Markov Chain
Monte Carlo (MCMC) methods for learning EBMs. To sidestep this issue, the
above work proposes a joint model $P_{\theta}(x)\propto
P_{LM}(x)\exp\left(-E_{\theta}(x)\right)$ with which the energy function
$E_{\theta}(x)$ can be learned using _noise contrastive estimation (NCE)_ by
training a binary classifier to discriminate text generated by $P_{LM}$ and
text samples from the target data distribution. The quality of generated text
from $P_{LM}$ can then be improved by resampling according to the learned
energy function.

While theoretically $P_{\theta}(x)$ can be trained to approximate the target
data distribution $P_{\text{data}}$, given that
$\text{Supp}(P_{LM})=\text{Supp}(P_{\text{data}})$, the practical
effectiveness of this approach is constrained by $P_{LM}$’s ability to
accurately capture the support of $P_{\text{data}}$. The study by Bakhtin et
al. (2019) reveals that it is easy to find text of low probability under
$P_{LM}$ but high probability under $E_{\theta}(x)$ by perturbing the text
generated by $P_{LM}$. This observation suggests that the learned energy model
has a multitude of _spurious modes_ in the residual space of $P_{LM}$
(although their impact on the aforementioned resampling strategy may be
limited, given that $P_{LM}$ does not operate in this residual space).

In this work, we propose an _adversarial training strategy_ to address
spurious modes and assess its efficacy in enhancing the sequence generation
abilities of $P_{LM}$. Unlike the aforementioned approach that trains a binary
classifier to distinguish between samples from the data distribution
$P_{\text{data}}$ and generated text from $P_{LM}$, our strategy involves
training the energy model to contrast samples of $P_{\text{data}}$ with
_adversarially perturbed_ versions of $P_{LM}$ samples. These perturbed
samples are generated through gradient-free adversarial attacks targeted at
the energy model within the discrete input space. This methodology, while
analogous to traditional MCMC-based maximum likelihood training protocols for
EBMs, deviates by substituting MCMC sampling with adversarial attacks for
generating negative samples, thereby circumventing the challenges associated
with gradient-based sampling in a discrete space. While there is no
theoretical guarantee for the convergence of adversarial attacks to the model
distribution, empirical evidence validates their efficacy in isolating
spurious modes and generating statistically representative samples from the
model distribution. By training the energy model to distinguish text samples
from the target distribution from these perturbed samples, we enable it to
suppress the spurious modes and more accurately allocate probability mass to
the support of the target data distribution. Consequently, this allows us to
produce higher-quality samples by first drawing from $P_{LM}$ and subsequently
applying adversarial perturbations to steer these base instances toward the
support of the target distribution.

We assess the efficacy of our proposed method using an arithmetic sequence
generation task, where it demonstrates a substantial performance enhancement
compared to the baseline approach presented in Deng et al. (2019).

##  2 Methodology

###  2.1 Notations and Definitions

Let $\mathcal{X}$ represent the input space and $\mathbf{x}\in\mathcal{X}$
denote a specific instance of a text sequence. Importantly, $\mathbf{x}$ can
be represented as a sequence of words $(x_{1},x_{2},\ldots,x_{T})$, where $T$
is the length of the sequence. We define $P_{LM}(\mathbf{x})$ as the
probability of $\mathbf{x}$ under the Locally Normalized Language Model
(LNLM). Further, $P_{\text{data}}(\mathbf{x})$ denotes the empirical
distribution of the target data. The ultimate objective of an Energy-Based
Model (EBM) is to learn an energy function $E_{\theta}(\mathbf{x})$,
parameterized by $\theta$, which aims to assign lower energy values to the
support of $P_{\text{data}}$ and higher energy values to regions outside the
support.

In the context of EBMs, probability distributions are formulated by linking
the energy function $E_{\theta}(x)$ to probabilities via a Gibbs distribution:

$p_{\theta}(x)=\frac{\exp\left(-E_{\theta}(x)\right)}{Z(\theta)},$ (1)

where the normalizing constant $Z(\theta)$, also known as the partition
function, is an integral over the unnormalized probabilities across all
states:

$Z(\theta)=\int\exp\left(-E_{\theta}(x)\right)\,dx.$ (2)

In many complex models, the partition function $Z(\theta)$ is computationally
intractable, making the direct application of maximum likelihood estimation
(MLE) for the model parameters $\theta$ impractical. Traditional MLE for EBMs
relies on computing the gradient of the log-likelihood function. Let
$p_{\text{data}}$ represent the distribution of the observed data; then the
gradient of the log-likelihood can be expressed as:

$\nabla_{\theta}\mathbb{E}_{x\sim p_{\text{data}}}\left[\log
p_{\theta}(x)\right]=\mathbb{E}_{x\sim
p_{\text{data}}}\left[\nabla_{\theta}(-E_{\theta}(x))\right]-\mathbb{E}_{x\sim
p_{\theta}(x)}\left[\nabla_{\theta}(-E_{\theta}(x))\right].$ (3)

Intuitively, maximizing the log-likelihood using this gradient serves to
reduce the energy function $E_{\theta}(x)$ on samples from $p_{\text{data}}$
while elevating it on samples drawn from $p_{\theta}$. Samples of $p_{\theta}$
are also known as _negative samples_ , and are typically generated with Markov
Chain Monte Carlo (MCMC) sampling. However, as mentioned earlier, the discrete
nature of text complicates the direct use of MCMC methods for text generation.

###  2.2 Generating Negative Samples With Adversarial Attacks

Given the challenges associated with sampling from $p_{\theta}$ in a discrete
textual space, we propose to generate negative samples by first drawing
samples from $P_{LM}$ and then performing adversarial attacks on these base
samples.

Let $x_{i}$ denote the $i^{th}$ token in the text sequence
$\mathbf{x}=(x_{1},x_{2},\ldots,x_{T})$, which is sampled from $P_{LM}$. The
corresponding token embedding for $x_{i}$ is represented as $\mathbf{e}_{i}$,
where $\mathbf{e}_{i}$ belongs to $\mathbb{R}^{d}$ and $d$ signifies the
dimensionality of the embedding space. To perform an adversarial attack, we
first identify the index of the token that most significantly influences the
energy function, as given by:

$i^{*}=\arg\max_{i=1,...,T}\|\nabla_{\mathbf{e}_{i}}(-E_{\theta}(x))\|_{1}$
(4)

Subsequently, we form a set of candidate sequences, denoted by $\mathcal{C}$,
by replacing $x_{i^{*}}$ with tokens from its synonym set
$\mathbb{L}(x_{i^{*}})$ within the vocabulary $\mathcal{V}$. Specifically,
each candidate sequence $\mathbf{x^{\prime}}$ is crafted by substituting
$x_{i^{*}}$ in the original sequence $\mathbf{x}$ with an alternative token
$v$ from $\mathbb{L}(x_{i^{*}})$. Mathematically, each candidate sequence
$\mathbf{x^{\prime}}$ can be expressed as:

$\mathbf{x^{\prime}}=(x_{1},...,x_{i^{*}-1},v,x_{i^{*}+1},...,x_{T})$ (5)

The adversarial example is then selected as the candidate sequence that
maximizes the negative energy function:

$\mathbf{x^{\prime}}^{*}=\arg\max_{\mathbf{x^{\prime}}\in\mathcal{C}}(-E_{\theta}(\mathbf{x^{\prime}}))$
(6)

In line with the iterative nature of MCMC methods used in EBM training, we
propose an iterative multi-step adversarial attack. Starting with an initial
text sequence, each iteration identifies the most influential token based on
the energy function’s gradient. A new adversarial sample is crafted by
replacing this token with a synonym to maximize the negative energy function.
This new sample serves as the starting point for the subsequent iteration. The
process continues iteratively, chaining together multiple attacks to
incrementally exploit the model’s spurious modes.

It should be noted that the above algorithm is specifically tailored for the
arithmetic sequence generation task, as outlined in Section 3.1, and shares
notable similarities with substitution-based text attacks prevalent in
existing literature. For real-world language modeling tasks, it may be
advantageous to employ more sophisticated attacks, such as those outlined in
Morris et al. (2020)

###  2.3 Training the Energy Function $E_{\theta}$

We outline the training algorithm employed to optimize the energy function
$E_{\theta}$ in Algorithm 1. The primary goal is to enable $E_{\theta}$ to
approximate the target distribution $p_{\text{data}}$ by training the model to
discriminate between samples drawn from $p_{\text{data}}$ and negative samples
generated with the attack algorithm as outlined in Section 2.2.

Algorithm 1 Training Sequence Energy-Based Model with Iterative Multi-Step
Adversarial Examples

1: Input: Number of epochs $N$, batch size $B$, learning rate $\alpha$, number
of attack steps $S$

2: Initialize: EBM parameters $\theta$

3: for epoch $=1,2,\ldots,N$ do

4: Sample mini-batch $\mathbf{X}_{\text{data}}$ from $p_{\text{data}}$

5: Sample mini-batch $\mathbf{X}_{\text{LM}}$ from $P_{\text{LM}}$

6: # Generate multi-step adversarial examples for each sample in mini-batch

7: for each $\mathbf{x}_{\text{LM}}$ in $\mathbf{X}_{\text{LM}}$ do

8: for step $=1,2,\ldots,S$ do

9: Compute
$i^{*}=\arg\max_{i=1,...,T}\|\nabla_{\mathbf{e}_{i}}(-E_{\theta}(\mathbf{x}_{\text{LM}}))\|_{1}$

10: Generate candidate sequences $\mathcal{C}$ by replacing $x_{i^{*}}$ with
synonyms

11: Compute
$\mathbf{x^{\prime}}^{*}=\arg\max_{\mathbf{x^{\prime}}\in\mathcal{C}}(-E_{\theta}(\mathbf{x^{\prime}}))$

12: Set $\mathbf{x}_{\text{LM}}\leftarrow\mathbf{x^{\prime}}^{*}$ for the next
step

13: end for

14: Store the final $\mathbf{x^{\prime}}^{*}$ in set $\mathbf{X^{\prime}}^{*}$

15: end for

16: # Update model parameters using the entire mini-batch

17: Compute gradient
$\nabla_{\theta}\sum_{\mathbf{x}_{\text{data}}\in\mathbf{X}_{\text{data}}}E_{\theta}(\mathbf{x}_{\text{data}})$
and
$\nabla_{\theta}\sum_{\mathbf{x^{\prime}}^{*}\in\mathbf{X^{\prime}}^{*}}E_{\theta}(\mathbf{x^{\prime}}^{*})$

18: Update $\theta$ using gradient descent:
$\theta\leftarrow\theta-\alpha\left(\nabla_{\theta}\sum_{\mathbf{x}_{\text{data}}\in\mathbf{X}_{\text{data}}}E_{\theta}(\mathbf{x}_{\text{data}})-\nabla_{\theta}\sum_{\mathbf{x^{\prime}}^{*}\in\mathbf{X^{\prime}}^{*}}E_{\theta}(\mathbf{x^{\prime}}^{*})\right)$

19: end for

###  2.4 Improving Generation With the Energy Function $E_{\theta}$

After training both $P_{LM}$ and $E_{\theta}$, we enhance the quality of
generated samples by initially drawing from $P_{LM}$ and subsequently applying
adversarial attacks to these samples via the iterative attack algorithm
described in Section 2.2. Intuitively, if $E_{\theta}$ is a faithful
approximation of the target distribution, it will assign lower energy values
specifically to the support of the target distribution. Therefore, when the
iterative adversarial attacks are employed to guide the perturbed samples
toward these low-energy regions, the samples are expected to align more
closely with the target distribution, thereby elevating their overall quality.

##  3 Experiments

###  3.1 Task

In this work, our primary aim is to empirically validate our proposed
adversarial training approach using a synthetic arithmetic sequence generation
task, identical to the sequence generation task investigated in Deng (2023).
The arithmetic sequences in this dataset are constructed based on elementary
addition operations, as demonstrated by the sequence format:

$a_{10}\ a_{9}\ a_{8}\ a_{7}\ a_{6}\ a_{5}\ a_{4}\ a_{3}\ a_{2}\ a_{1}+b_{10}\
b_{9}\ b_{8}\ b_{7}\ b_{6}\ b_{5}\ b_{4}\ b_{3}\ b_{2}\ b_{1}=c_{11}\ c_{10}\
c_{9}\ c_{8}\ c_{7}\ c_{6}\ c_{5}\ c_{4}\ c_{3}\ c_{2}\ c_{1}$

Here, $a_{i},b_{i}$, and $c_{i}$ denote individual digits ranging between 0
and 9. The $c$ sequence is a result of the element-wise addition of $a$ and
$b$, following conventional rules for arithmetic addition, including carrying
over, if needed:

$c_{i}=a_{i}+b_{i}\quad\text{for}\quad i=1,\ldots,10$ $c_{11}=\text{carry-over
from the last addition, defaults to 0}$

An example sequence looks like:

$4\ 6\ 8\ 0\ 1\ 5\ 6\ 8\ 2\ 3+2\ 3\ 9\ 0\ 9\ 2\ 2\ 6\ 4\ 5=0\ 7\ 3\ 3\ 2\ 2\
6\ 5\ 1\ 6\ 1$

Here, the sequence $a$ is $4\ 6\ 8\ 0\ 1\ 5\ 6\ 8\ 2\ 3$ and sequence $b$ is
$2\ 3\ 9\ 0\ 9\ 2\ 2\ 6\ 4\ 5$. We note that the language model $P_{LM}$ is
trained to generate the whole sequence $a+b=c$, as opposed to only $c$.

By employing this synthetic task, we can precisely assess the quality of
generated sequences through their compliance with the arithmetic rules of
addition. To quantify performance, we employ an average correctness metric
computed over multiple generated sequences. The metric is particularly
tailored to the characteristics of the arithmetic sequence generation task.
The correctness for a single generated sequence is calculated as:

$\text{correctness}=\frac{\text{number of correctly predicted digits in
}c}{\text{total number of digits in }c}$

Here, a digit in sequence $c$ is considered “correctly predicted” if it
exactly coincides with the corresponding digit in the ground truth sequence.
Subsequently, the average correctness over $N$ sequences is given by:

$\text{Average correctness}=\frac{1}{N}\sum_{i=1}^{N}\text{correctness}_{i}$

where $\text{correctness}_{i}$ is the correctness of the $i^{th}$ sequence.
Thus, a perfect model would attain an average correctness of 1, indicating
flawless sequence generation across all test instances.

By using this synthetic dataset, we aim to empirically demonstrate the
efficacy of our proposed adversarial training approach in mitigating spurious
modes, thereby enhancing the capabilities of $P_{LM}$ in sequence generation
tasks.

###  3.2 Training

We set up our training protocol based on Deng (2023). The full training
details for $P_{LM}$ and $E_{\theta}$ (pretraining, after attaching a language
model head) can be found in Table 1. In terms of the training objective, both
models are trained by maximizing the conditional likelihood of the next token
$x_{i+1}$ given the sequence of previous tokens $x_{1},x_{2},\ldots,x_{i}$:

$\mathcal{L}(\theta)=\sum_{i=1}^{T}\log P(x_{i+1}\mid
x_{1},x_{2},\ldots,x_{i};\theta).$ (7)

Following pretraining, $E_{\theta}$ undergoes a fine-tuning process aimed at
approximating the target distribution $P_{\text{data}}$. In the baseline
approach of Residual-EBM (Deng et al., 2019), $E_{\theta}$ is fine-tuned using
the conventional binary cross-entropy loss function. This fine-tuning process
utilizes 1000 samples drawn from $P_{\text{data}}$ and 1000 synthetic samples
from $P_{LM}$, spanning 5 epochs with early stopping. In contrast, our
proposed methodology fine-tunes $E_{\theta}$ according to Algorithm 1,
employing the same number of $P_{\text{data}}$ and $P_{LM}$ samples but
extending the training to 55 epochs, also with early stopping. We use 5 attack
steps in Algorithm 1.

Owing to the inherent stochasticity involved in the training procedure, we
note fluctuations in model performance when utilizing Algorithm 1. To account
for this variability, we conduct the experiment ten times and report both
individual and aggregate performance metrics.

Table 1: Training setups for $P_{LM}$ and $E_{\theta}$ (pretraining)

| $P_{\text{LM}}$ |  $E_{\theta}$ (pretraining)
---|---|---
Architecture |  | One-layer LSTM
---
Embedding Size = 256
Hidden Layer Size = 256
| Two-layer Transformer
---
Attention Heads = 4
Embedding Size = 128
Hidden Layer Size = 128
Optimization |  | Optimizer = AdamW
---
Learning Rate = 0.001
Weight Decay = 0.05
Epochs = 500
| Optimizer = AdamW
---
Learning Rate = $6\times 10^{-4}$
Epochs = 500
Training Samples | 1,000 | 100,000

###  3.3 Results

Fig. 1 shows the performances of the proposed approach to learning EBMs with
adversarial training and the baseline approach of Residual-EBM (Deng et al.,
2019). The Residual-EBM methodology employs importance sampling to produce
novel sequence samples. To create a single sequence instance, the method
initially draws multiple instances from $P_{LM}$. Subsequently, it reweighs
these instances using the learned energy function $E_{\theta}$. Finally, it
performs resampling based on the computed instance weights. Fig. 1(b)
illustrates the average correctness of generated sequences as a function of
the number of instances sampled in this procedure. It can be seen that the
importance sampling can only marginally improve the quality of generated
sequences. In contrast, our approach significantly improves the quality of
generated sequences as the number of attack steps increases.

(a) Average correctness of generated sequences as the number of attack steps
varies

(b) Average correctness of the generated sequence as the number of samples (in
resampling) varies

Figure 1: Performances of the proposed approach to learning EMBs with
adversarial training and the baseline approach of Residual-EBM

####  3.3.1 Ablation Study on Effect of Adversarial Training

In this ablation study we investigate the impact of adversarial training on
the energy function. Specifically, we examine two distinct training paradigms:

  1. 1.

The Proposed Approach: EBM with Adversarial Training — In this setting, the
energy function $E_{\theta}$ is trained using samples from $P_{\text{data}}$
as positive samples and _adversarially perturbed_ samples from $P_{LM}$ as
negative samples. Further details can be found in Algorithm 1.

  2. 2.

Baseline: EBM without Adversarial Training — Here, the energy function
$E_{\theta}$ is trained using samples from $P_{\text{data}}$ as positive
samples and samples from $P_{LM}$ as negative samples. Importantly, we do not
apply adversarial perturbations to samples from $P_{LM}$, making it
functionally equivalent to setting the number of adversarial attack steps $S$
to zero in Algorithm 1.

After training $E_{\theta}$, we generate test samples by first drawing from
$P_{LM}$ and then applying adversarial perturbations to the sampled data. The
quality of the test samples is then evaluated using the average correctness
metric as specified in Section 3.1.

Fig. 2 illustrates the average correctness of the generated test samples as a
function of the number of attack steps. In the absence of adversarial
training, the quality of the generated samples experiences a marginal
improvement after a few attack steps. In contrast, when attacking the EBM that
has undergone adversarial training, the quality of the generated samples sees
a substantial improvement.

Figure 2: Average correctness of perturbed sequences as the number of attack
steps are varied

##  4 Conclusion

This work demonstrates the potential for adversarial training techniques to
enhance the training of energy-based models for discrete sequence generation
tasks. By employing adversarial attacks to generate negative samples, the
proposed methodology enables more effective suppression of spurious modes
compared to prior binary classification-based training protocols.

However, the results also highlight key limitations and challenges associated
with adversarial training that warrant further investigation. Most notably,
the adversarial perturbation process introduces significant computational
overhead compared to standard training paradigms. Each iteration of the attack
algorithm requires evaluating the gradient of the energy function across all
tokens in the sequence, followed by synonym substitutions and re-evaluation of
the energy function. As the sequence length and size of the synonym set grow,
these computations can become prohibitively expensive. Developing efficient
and scalable adversarial training protocols is an important area for future
work towards enabling the application of discrete EBMs to complex, real-world
NLG tasks.

In conclusion, this preliminary study demonstrates promising results but also
highlights algorithmic and computational challenges associated with
adversarial discrete EBM training that necessitate continued research and
innovation. With further progress in this emerging field, adversarial
techniques may prove instrumental in enhancing the generative capacities of
discrete sequential models.

## References

  * Bakhtin et al. [2019] A. Bakhtin, S. Gross, M. Ott, Y. Deng, M. Ranzato, and A. Szlam.  Real or fake? learning to discriminate machine from human generated text.  _arXiv preprint arXiv:1906.03351_ , 2019.
  * Deng [2023] Y. Deng.  Residual-ebm: Code for residual energy-based models.  https://github.com/da03/Residual-EBM/tree/main, 2023.  GitHub repository.
  * Deng et al. [2019] Y. Deng, A. Bakhtin, M. Ott, A. Szlam, and M. Ranzato.  Residual energy-based models for text generation.  In _International Conference on Learning Representations_ , 2019.  URL https://arxiv.org/pdf/2004.11714.pdf.
  * Morris et al. [2020] J. X. Morris, E. Lifland, J. Y. Yoo, J. Grigsby, D. Jin, and Y. Qi.  Textattack: A framework for adversarial attacks, data augmentation, and adversarial training in nlp, 2020.
  * Parshakova et al. [2019] T. Parshakova, J.-M. Andreoli, and M. Dymetman.  Global autoregressive models for data-efficient sequence learning.  _arXiv preprint arXiv:1909.07063_ , 2019.  URL https://arxiv.org/pdf/1909.07063.pdf.
  * Wang and Ou [2018] B. Wang and Z. Ou.  Learning neural trans-dimensional random field language models with noise-contrastive estimation.  In _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 6134–6138. IEEE, 2018.