11institutetext: Institute of Theoretical Physics
Faculty of Mathematics and Physics, Charles University in Prague
V Holešovičkách 2, 180 00 Praha 8, Czech Republic
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Averaging in LRS class II spacetimes

Petr Kašpar† Otakar Svítek§

(Received: date / Accepted: date)

###### Abstract

We generalize Buchert’s averaged equations [Gen. Rel. Grav. 32, 105 (2000);
Gen. Rel. Grav. 33, 1381 (2001)] to LRS class II dust model in the sense that
all Einstein equations are averaged, not only the trace part. We derive the
relevant averaged equations and we investigate backreaction on expansion and
shear scalars in an approximate LTB model. Finally we propose a way to close
the system of averaged equations.

###### Keywords:

LRS family; Cosmology; Averaging

##  1 Introduction

Our universe is considered to be homogeneous and isotropic on the large scale
leading to the FLRW model. However, if we move to smaller scales, we can
observe a strongly inhomogeneous distribution of structures. If we want to
deal with inhomogeneity rigorously and at the same time keep a consistent
connection with the FLRW geometry we may consider an averaging formalism to
smooth out the metric tensor and at the same time average Einstein equations
as well. The problem is that Einstein equations are nonlinear and if we
average them straightforwardly we do not obtain an averaged metric tensor as a
solution of averaged equations. Instead, we should consider an additional term
— the so-called correlation term, which can change the evolution of a smooth
metric tensor and lead to the so-called backreaction. This term arises due to
the nonlinearity of Einstein equations. It does not need to satisfy the usual
energy conditions so it can possibly act as dark energy.

While building a rigorous averaging scheme we face the problem that the
average value of a tensor field is not well defined. There are several
different approaches to define averages of tensors. One of the most promising
ones is the scheme by Zalaletdinov Zal1 , Zal2  where not only Einstein
equations but also Cartan structure equations (and their integrability
conditions) are averaged. A theorem about isometric embedding of a 2-sphere
into Euclidian space is applied in the averaging method developed by Korzyński
Korz . In Brannlund  Weitzenböck connection for parallel transport is used to
define the average value of a tensor field.

One of the most popular approaches to averaging is the one investigated by
Buchert Buchert1 , Buchert2 , where only the scalar part of Einstein equations
is averaged. Wiltshire used this approach to give an alternative explanation
of cosmic acceleration Wiltshire . This theory was also applied to the
cosmological perturbation theory Li1 , Li2 , Behrend , Clarkson . For
observational issues see e.g. Larena . In this paper, we will generalize
Buchert’s equations to the locally rotationally symmetric (LRS) class II dust
family of spacetimes. The LRS family was classified in dust , Stewart  and
recently in LRS . It contains e.g. LRS Bianchi cosmologies, Kantowski-Sachs
model or LTB model and its generalizations. We will use the fact that this
family is described by scalars to average the complete set of Einstein
equations, including constraints. Although the averaged constraints are shown
to be preserved during evolution the averaged system of equations is not
closed and additional information has to be supplemented.

In the past there have been many attempts to apply Buchert’s approach to the
LTB model. Papers Sussman  and Chuang  are comprehensive studies of Buchert’s
formalism applied to generic LTB models. For a treatment of Buchert’s
equations inside LTB spacetime see e.g. Paranjape1 , Rasanen1 , Bolejko2 ,
Sussman2 , Sussman3 , Mattsson1 , Mattsson2  and Mattsson3 . For an
application of Buchert’s formalism to the structure formation see e.g.
Rasanen2 , Rasanen3  and Paranjape2 .

Naturally, one can study inhomogeneities perturbatively on a homogeneous
background and many important results are based on this approach. However, we
should be cautious about relying solely on a linear perturbative analysis when
dealing with a nonlinear theory. The effects of the correlation term indicate
what kind of effects one might be missing when using a simple approach. In
this sense, rigorous averaging of exact inhomogeneous spacetimes leading to
standard cosmological models provides a possibility to qualitatively estimate
these effects.

The paper is organized as follows. In Section 2 we review the LRS family and
its characterizations, then we briefly mention Buchert’s equations. Next, we
average equations describing dust LRS class II family. After a short review of
LTB metric in Section 5 we investigate the backreaction in the so-called onion
model. We proceed by attempting to close the averaged equations and we finish
with conclusion.

##  2 LRS family

Locally rotationally symmetric (LRS) dust spacetimes are defined by the
following features dust : In an open neighborhood of each point $p$, there
exists a nondiscrete subgroup of the Lorentz group which leaves invariant the
Riemann tensor and its covariant derivatives up to the third order. Therefore,
in LRS spacetimes there exists a preferred direction $e^{\mu}$ (the axis of
symmetry) at every point. The subgroup can be one or three-dimensional. In the
latter case, we can rotate the axis of symmetry and spacetimes are everywhere
isotropic - these are the FLRW models.

We will use the covariant 1+3 splitting of spacetime with the timelike vector
field $u^{\mu}$ normalized by the condition $u_{\rho}u^{\rho}=-1$ and the
projection tensor $h_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}$. In this section we
will follow the article of van Elst and Ellis LRS .

The preferred spacelike vector field $e^{\mu}$ satisfies the following
conditions:

$e_{\rho}u^{\rho}=0,\,\,\,\,\,\,e_{\rho}e^{\rho}=1.$ (1)

Because of the defining property of the LRS spacetime, all covariantly defined
spacelike vectors orthogonal to $u^{\mu}$ (acceleration $\dot{u}^{\mu}$,
vorticity $\omega^{\mu}$, projected gradient of density
$h_{\,\,\mu}^{\sigma}\nabla_{\sigma}\rho$, pressure
$h_{\,\,\mu}^{\sigma}\nabla_{\sigma}p$ and expansion
$h_{\,\,\mu}^{\sigma}\nabla_{\sigma}\theta$) must be proportional to $e^{\mu}$
\- if this condition does not hold, spacelike vectors will not be invariant
under the rotation about $e^{\mu}$.

$\dot{u}^{\mu}=\dot{u}e^{\mu},\,\,\,\,\,\,\,\,\,\ \omega^{\mu}=\omega
e^{\mu},$ (2)

$h_{\,\,\mu}^{\sigma}\nabla_{\sigma}\rho=\rho^{\prime}e_{\mu},\,\,\,\,\,\
h_{\,\,\mu}^{\sigma}\nabla_{\sigma}p=p^{\prime}e_{\mu},\,\,\,\,\,\
h_{\,\,\mu}^{\sigma}\nabla_{\sigma}\theta=\theta^{\prime}e_{\mu}.$ (3)

Dots denote covariant derivative along the flow vector $u^{\mu}$ and primes
denote covariant derivative along the vector $e^{\mu}$. We define the
magnitude of the spatial rotation $k$ and the magnitude of the spatial
divergence $a$ by

$k:=\left|{\eta^{\alpha\beta\gamma\delta}\left({\nabla_{\beta}e_{\gamma}}\right)u_{\delta}}\right|,$
(4)

$a:=h_{\,\,\beta}^{\alpha}\left({\nabla_{\alpha}e^{\beta}}\right),$ (5)

where $\eta^{\alpha\beta\gamma\delta}$ is the totally antisymmetric Levi-
Civita pseudotensor ($\eta_{0123}=-\sqrt{-g}$). A similar rule works also for
the spacelike tracefree symmetric tensors orthogonal to $u^{\mu}$. Following
LRS  we introduce a tensor field $e_{\mu\nu}$ defined by $e^{\mu}$

$e_{\mu\nu}:={1\over 2}\left({3e_{\mu}e_{\nu}-h_{\mu\nu}}\right).$ (6)

Then we have the relations for the shear tensor and the electric and magnetic
parts of the Weyl tensor

$\sigma_{\mu\nu}={2\over{\sqrt{3}}}\sigma
e_{\mu\nu},\,\,\,\,\,\,\,\,E_{\mu\nu}={2\over{\sqrt{3}}}Ee_{\mu\nu},\,\,\,\,\,\,\,\,H_{\mu\nu}={2\over{\sqrt{3}}}He_{\mu\nu}.$
(7)

Here we can see that the LRS spacetimes are characterized only by a finite set
of scalar functions.

For simplicity we will restrict our attention to the LRS class II dust models
defined by the relation $k=\omega=0$. One can also show that the magnetic part
of the Weyl tensor is equal to zero, $H=0$. This family of spacetimes includes
the LTB metric and its generalizations based on foliation by spacelike
2-surfaces with negative or zero curvature scalar. The relevant evolution
equations are

$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle-{1\over
3}\theta^{2}-2\sigma^{2}-4\pi\rho,$ (8) $\displaystyle\dot{\sigma}$
$\displaystyle=$ $\displaystyle-{1\over{\sqrt{3}}}\sigma^{2}-{2\over
3}\theta\sigma-E,$ (9) $\displaystyle\dot{E}$ $\displaystyle=$
$\displaystyle-4\pi\rho\sigma+\sqrt{3}E\sigma-\theta E,$ (10)
$\displaystyle\dot{\rho}$ $\displaystyle=$ $\displaystyle-\rho\theta,$ (11)
$\displaystyle\dot{a}$ $\displaystyle=$ $\displaystyle-{1\over
3}a\theta+{1\over{\sqrt{3}}}a\sigma,$ (12)

and the constraints

$\displaystyle\sigma^{\prime}$ $\displaystyle=$
$\displaystyle{1\over{\sqrt{3}}}\theta^{\prime}-{2\over 3}a\sigma,$ (13)
$\displaystyle E^{\prime}$ $\displaystyle=$ $\displaystyle-{3\over
2}aE+{{4\pi}\over{\sqrt{3}}}\rho^{\prime},$ (14) $\displaystyle a^{\prime}$
$\displaystyle=$ $\displaystyle{2\over
9}\theta^{2}+{2\over{3\sqrt{3}}}\theta\sigma-{4\over
3}\sigma^{2}-{2\over{\sqrt{3}}}E-{1\over 2}a^{2}-{{16\pi}\over 3}\rho.$ (15)

If we take the time derivative of the constraints, we can prove that they do
not change with time.

##  3 Buchert’s equations

In this section we will review an averaging method developed by Buchert
Buchert1 . We will consider only the dust case - for generalization to perfect
fluid see Buchert2 . This approach uses 1+3 splitting of spacetime, which is
well defined by irrotational dust 4-velocity. However, averaging is well
defined only for scalars, therefore only scalar part of Einstein equations is
averaged. Given a scalar field $A$, the average value over three-dimensional
spacelike domain ${\cal D}$ is defined by

$\left\langle A\right\rangle_{\cal D}={1\over{V_{\cal D}}}\int\limits_{\cal
D}{d^{3}XJA=}{1\over{V_{\cal D}}}\int\limits_{\cal D}{d^{3}X\sqrt{\det
g_{ij}}A},$ (16)

where $J:=\sqrt{detg_{ij}}$, $g_{ij}$ is the metric of the spacelike
hypersurface, $X^{i}$ are the comoving coordinates and $V_{\cal D}$ is the
proper volume of the three-dimensional domain ${\cal D}$. From this definition
we can see that time derivative and averaging do not commute. We have a
commutation relation

$\displaystyle\left\langle A\right\rangle_{\cal D}^{\cdot}$ $\displaystyle=$
$\displaystyle{d\over{dt}}\left({{1\over{V_{\cal D}}}\int\limits_{\cal
D}{d^{3}XJA}}\right)=-{{\dot{V}_{\cal D}}\over{V_{\cal D}}}\left\langle
A\right\rangle_{\cal D}+{1\over{V_{\cal D}}}\int\limits_{\cal
D}{d^{3}X\left({\dot{J}A+J\dot{A}}\right)}$ (17) $\displaystyle=$
$\displaystyle-\left\langle\theta\right\rangle_{\cal D}\left\langle
A\right\rangle_{\cal D}+\left\langle{A\theta}\right\rangle_{\cal
D}+\left\langle{\dot{A}}\right\rangle_{\cal D},$

where the expansion rate $\theta$ is related to the velocity of the fluid
$u^{\mu}$ according to the definition $\theta=u^{\mu}_{\,\,\,;\mu}$. Next, in
analogy with the FLRW spacetime we introduce the dimensionless scale factor
$a_{\cal D}$ and the effective Hubble parameter $H_{\cal D}$

${a_{\cal D}}=\left(\frac{V_{\cal D}}{V_{{\cal D}_{i}}}\right)^{\frac{1}{3}},$
(18) $\left\langle\theta\right\rangle_{\cal D}=\frac{\dot{V_{\cal D}}}{V_{\cal
D}}=3\frac{\dot{a_{\cal D}}}{a_{\cal D}}=:3H_{\cal D}.$ (19)

$V_{{\cal D}_{i}}$ is the volume of the initial domain which is geodetically
evolved into $V_{\cal D}$. Now we have a formalism for averaging scalars. To
obtain a scalar equation from Einstein equations, we have to contract them
with available tensors \- i.e. $g^{\mu\nu},u^{\mu}$ and $\nabla^{\mu}$. After
contracting we get the Raychaudhuri equation, the Hamiltonian constraint and
the continuity equation. Now we perform averaging and use the commutation rule
(17)

$3\frac{\ddot{a_{\cal D}}}{a_{\cal D}}+4\pi
G\left\langle\rho\right\rangle_{\cal D}={\cal Q}_{\cal D},$ (20)
$\left(\frac{\dot{a}_{\cal D}}{a_{\cal D}}\right)^{2}-\frac{8\pi
G}{3}\left\langle\rho\right\rangle_{\cal D}+\frac{\left\langle{\cal
R}\right\rangle_{\cal D}}{6}=-\frac{{\cal Q}_{\cal D}}{6},$ (21)
$\partial_{t}\left\langle\rho\right\rangle_{\cal D}+3\frac{\dot{a_{\cal
D}}}{a_{\cal D}}\left\langle\rho\right\rangle_{\cal D}=0.$ (22)

${\left\langle{\cal R}\right\rangle_{\cal D}}$ denotes the average value of
the spatial Ricci scalar, ${\left\langle{\rho}\right\rangle_{\cal D}}$ means
the average density of the fluid and ${\cal Q}_{\cal D}$ that shows possible
backreaction (due to inhomogeneity and anisotropy) is defined by

${\cal Q}_{\cal D}:={2\over
3}\left\langle{\left({\theta-\left\langle\theta\right\rangle_{\cal
D}}\right)^{2}}\right\rangle_{\cal
D}-2\left\langle{\sigma^{2}}\right\rangle_{\cal D}.$ (23)

The scalar $\sigma^{2}={1\over 2}\sigma_{\mu\nu}\sigma^{\mu\nu}$ is
constructed from the shear tensor.

##  4 Averaging LRS class II dust spacetime

Now, we will generalize the above approach to LRS class II dust solutions.
Originally, Buchert considered spacetimes with a dust Buchert1  or a perfect
fluid Buchert2  source. He did not assume any symmetries or simplifications
and his equations can be applied to a large class of metrics. Here we will
restrict to spacetimes with the special LRS symmetry. For this family we will
generalize Buchert’s equations in the sense that all Einstein equations are
averaged consistently.

Given a preferred spacelike direction $e_{\mu}$, all the equations describing
the LRS metric are scalar. It means that we can perform averaging (which is
covariantly defined for scalars). We will define averaging over the spacelike
domain ${\cal D}$ according to (16). In order to obtain averaged equations we
need to derive commutation relations for the time and spatial derivatives
(with respect to the preferred direction). For the LRS class II spacetime it
is possible to express derivative along the preferred direction $\vec{e}$ by
the formula $\vec{e}=\sqrt{g_{11}}\partial_{1}$, where $g_{11}$ is the metric
function of a particular solution inside the LRS class II spacetime. The most
studied solution inside this class is the LTB model, for which the square root
of the metric function reads $\sqrt{g_{11}}=\frac{R^{\prime}}{\sqrt{1+2B}}$.
We will show the basic facts about the LTB metric in the next section. The
commutation rule between coordinate derivative and averaging reads

$\partial_{1}\left\langle
A\right\rangle-\left\langle\partial_{1}A\right\rangle=\left\langle
A\frac{\partial_{1}J}{J}\right\rangle-\left\langle
A\right\rangle\left\langle\frac{\partial_{1}J}{J}\right\rangle.$ (24)

For simplicity we will restrict to the class II LRS spacetime with the
condition $p=0\Leftrightarrow\dot{\rho}=-\rho\theta$ (dust models) which
includes LTB spacetimes and their generalizations. For simpler notation we
shall omit symbol $\cal D$ at the averaging bracket in the rest of the paper
(however we retain the symbols not defined as direct averages - $a_{\cal
D},{\cal Q}_{\cal D}$). If we average the equations (8) \- (15) we obtain

$\displaystyle\left\langle\theta\right\rangle^{\cdot}$ $\displaystyle=$
$\displaystyle-{1\over
3}\left\langle\theta\right\rangle^{2}-4\pi\left\langle\rho\right\rangle+\underline{{2\over
3}\left({\left\langle{\theta^{2}}\right\rangle-\left\langle\theta\right\rangle^{2}}\right)-2\left\langle{\sigma^{2}}\right\rangle}\
,$ (25) $\displaystyle\left\langle\sigma\right\rangle^{\cdot}$
$\displaystyle=$
$\displaystyle-{1\over{\sqrt{3}}}\left\langle\sigma\right\rangle^{2}-{2\over
3}\left\langle\theta\right\rangle\left\langle\sigma\right\rangle-\left\langle
E\right\rangle+\underline{{1\over{\sqrt{3}}}\left({\left\langle\sigma\right\rangle^{2}-\left\langle{\sigma^{2}}\right\rangle}\right)}$
(26) $\displaystyle+\underline{{1\over
3}\left({\left\langle{\theta\sigma}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma\right\rangle}\right)}\
,$ $\displaystyle\left\langle E\right\rangle^{\cdot}$ $\displaystyle=$
$\displaystyle-4\pi\left\langle\rho\right\rangle\left\langle\sigma\right\rangle+\sqrt{3}\left\langle
E\right\rangle\left\langle\sigma\right\rangle-\left\langle\theta\right\rangle\left\langle
E\right\rangle$ (27)
$\displaystyle\underline{-4\pi\left({\left\langle{\rho\sigma}\right\rangle-\left\langle\rho\right\rangle\left\langle\sigma\right\rangle}\right)+\sqrt{3}\left({\left\langle{E\sigma}\right\rangle-\left\langle
E\right\rangle\left\langle\sigma\right\rangle}\right)}\ ,$
$\displaystyle\left\langle\rho\right\rangle^{\cdot}$ $\displaystyle=$
$\displaystyle-\left\langle\rho\right\rangle\left\langle\theta\right\rangle\
,$ (28) $\displaystyle\left\langle a\right\rangle^{\cdot}$ $\displaystyle=$
$\displaystyle-{1\over 3}\left\langle
a\right\rangle\left\langle\theta\right\rangle+{1\over{\sqrt{3}}}\left\langle
a\right\rangle\left\langle\sigma\right\rangle+\underline{{2\over
3}\left({\left\langle{a\theta}\right\rangle-\left\langle
a\right\rangle\left\langle\theta\right\rangle}\right)}$ (29)
$\displaystyle+\underline{{1\over{\sqrt{3}}}\left({\left\langle{a\sigma}\right\rangle-\left\langle
a\right\rangle\left\langle\sigma\right\rangle}\right)}\ ,$
$\displaystyle\left\langle\sigma\right\rangle^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{3}}\left\langle\theta\right\rangle^{\prime}-\frac{2}{3}\left\langle
a\right\rangle\left\langle\sigma\right\rangle+\underline{\sqrt{g_{11}}\left(\left\langle\sigma\frac{\partial_{1}J}{J}\right\rangle-\left\langle\sigma\right\rangle\left\langle\frac{\partial_{1}J}{J}\right\rangle\right)}$
$\displaystyle-$
$\displaystyle\underline{\sqrt{\frac{g_{11}}{3}}\left(\left\langle\theta\frac{\partial_{1}J}{J}\right\rangle+\left\langle\theta\right\rangle\left\langle\frac{\partial_{1}J}{J}\right\rangle\right)-\frac{2}{3}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}a\sigma\right\rangle-\left\langle
a\right\rangle\left\langle\sigma\right\rangle\right)}\ ,$
$\displaystyle\left\langle E\right\rangle^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left\langle a\right\rangle\left\langle
E\right\rangle+\frac{4\pi}{\sqrt{3}}\left\langle\rho\right\rangle^{\prime}+\underline{\sqrt{g_{11}}\left(\left\langle
E\frac{\partial_{1}J}{J}\right\rangle-\left\langle
E\right\rangle\left\langle\frac{\partial_{1}J}{J}\right\rangle\right)}$
$\displaystyle-$
$\displaystyle\underline{4\pi\sqrt{\frac{g_{11}}{3}}\left(\left\langle\rho\frac{\partial_{1}J}{J}\right\rangle-\left\langle\rho\right\rangle\left\langle\frac{\partial_{1}J}{J}\right\rangle\right)-\frac{3}{2}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}aE\right\rangle-\left\langle
a\right\rangle\left\langle E\right\rangle\right)},$ $\displaystyle\left\langle
a\right\rangle^{\prime}$ $\displaystyle=$
$\displaystyle\frac{2}{9}\left\langle\theta\right\rangle^{2}+\frac{2}{3\sqrt{3}}\left\langle\theta\right\rangle\left\langle\sigma\right\rangle-\frac{4}{3}\left\langle\sigma\right\rangle^{2}-\frac{2}{\sqrt{3}}\left\langle
E\right\rangle-\frac{1}{2}\left\langle
a\right\rangle^{2}-\frac{16\pi}{3}\left\langle\rho\right\rangle$
$\displaystyle+$ $\displaystyle\underline{\sqrt{g_{11}}\left(\left\langle
a\frac{\partial_{1}J}{J}\right\rangle-\left\langle
a\right\rangle\left\langle\frac{\partial_{1}J}{J}\right\rangle\right)+\frac{2}{9}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}\theta^{2}\right\rangle-\left\langle\theta\right\rangle^{2}\right)}$
$\displaystyle+$
$\displaystyle\underline{\frac{2}{3\sqrt{3}}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}\theta\sigma\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma\right\rangle\right)-\frac{4}{3}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}\sigma^{2}\right\rangle-\left\langle\sigma\right\rangle^{2}\right)}$
$\displaystyle-$
$\displaystyle\underline{\frac{2}{\sqrt{3}}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}E\right\rangle-\left\langle
E\right\rangle\right)-\frac{1}{2}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}a^{2}\right\rangle-\left\langle
a\right\rangle^{2}\right)}$ $\displaystyle-$
$\displaystyle\underline{\frac{16\pi}{3}\left(\sqrt{g_{11}}\left\langle\frac{1}{\sqrt{g_{11}}}\rho\right\rangle-\left\langle\rho\right\rangle\right)}\
.$

The underlined parts of the equations denote additional terms due to
averaging. We can recognize the well-known Buchert’s equation (25) with the
kinematical backreaction term and the mass conservation equation (28).

If we want to restrict the above equations to the LTB model (the line element
is given in the next section) we can substitute the following expressions for
the magnitude of the spatial divergence $a$, its average value $\left\langle
a\right\rangle_{\cal D}$ and the square root of the spatial part of the metric
$J=\sqrt{detg_{ij}}$

$a=h^{\nu}_{\mu}e^{\mu}_{\,\,\,\,;\nu}=\frac{2\sqrt{1+2B}}{R},\,\,\,\,\left\langle
a\right\rangle_{\cal D}=\frac{4\pi R^{2}}{V_{\cal
D}},\,\,\,\,J=\frac{R^{\prime}R^{2}\sin{\theta}}{\sqrt{1+2B}}.$ (33)

One can show that the averaged constraint equations (4) \- (4) are preserved
in time. The key role in the calculation is played by the equation LRS

$\left({f^{\prime}}\right)^{\cdot}=\left({\dot{f}}\right)^{\prime}-{2\over{\sqrt{3}}}\sigma
f^{\prime}-{1\over 3}\theta f^{\prime},$ (34)

and its averaged version. Now, we can take time derivative of the constraints
(4) \- (4). Using commutation rules and equations (8) \- (15), a slow but
straightforward computation will show that the constraints do not evolve in
time. The explicit computation for equation (4) is shown in the Appendix.

All of Einstein equations are averaged now. It means that we can investigate
not only backreaction on the expansion rate but also on shear scalar or
electric part of the Weyl scalar. The problem is that the equations are not
closed. We need additional relations to close the system because for example
$\left\langle\theta\right\rangle$ is independent of
$\left\langle\theta^{2}\right\rangle$. In the next chapters we will give some
suggestions for closing the system of equations.

So far we have not seen an analog of the Hamiltonian constraint. For example
Sussman in Sussman  used three-dimensional curvature $\cal R$ instead of the
function $a$. The relation between three-dimensional curvature $\cal R$ and
the magnitude of the spatial divergence $a$ reads LRS

${\cal R}=-\left(2a^{\prime}+\frac{2}{3}a^{2}-2K\right),$ (35)

where $K$ is the Gaussian curvature of the 2-D spacelike group orbits
orthogonal to $e^{\mu}$ and $u^{\mu}$. The Hamiltonian constraint has the form

${\cal R}=16\pi\rho-\frac{2}{3}\theta^{2}+2\sigma^{2}.$ (36)

The averaged Hamiltonian constraint is Buchert’s equation (21)

$\displaystyle\left\langle{\cal R}\right\rangle$ $\displaystyle=$
$\displaystyle
16\pi\left\langle\rho\right\rangle-\frac{2}{3}\left\langle\theta\right\rangle^{2}+2\left\langle\sigma\right\rangle^{2}$
(37)
$\displaystyle\underline{-\frac{2}{3}\left(\left\langle\theta^{2}\right\rangle-\left\langle\theta\right\rangle^{2}\right)+2\left(\left\langle\sigma^{2}\right\rangle-\left\langle\sigma\right\rangle^{2}\right)}.$

##  5 LTB metric

The most important representative of the dust LRS class II family is the LTB
spacetime. In this section we will briefly review its properties.

The Lemaître-Tolman-Bondi (LTB) metric Lemaitre , Tolman , Bondi  is a
spherically symmetric exact solution of Einstein equations. It corresponds to
an inhomogeneous dust with the stress energy tensor

$T_{\mu\nu}=\rho u_{\mu}u_{\nu},$ (38)

where $u_{\mu}$ is 4-velocity of the dust with density $\rho$. For a recent
review of LTB metric, see e.g. Bolejko , Hellaby . The line element reads

$ds^{2}=-dt^{2}+\frac{(R^{\prime})^{2}}{1+2B(r)}dr^{2}+R^{2}(t,r)[d\theta^{2}+\mbox{sin}^{2}(\theta)d\phi^{2}],$
(39)

where $B(r)$ is an arbitrary function and the prime denotes partial derivative
with respect to $r$. Function $R(t,r)$ is a solution of Einstein equation

$R_{,t}^{2}=2B+\frac{2M}{R}+\frac{\Lambda}{3}R^{2}.$ (40)

$M=M(r)$ is another arbitrary function of integration. The energy density
$\rho$ is determined by the equation

$4\pi\rho=\frac{M^{\prime}}{R^{\prime}R^{2}}.$ (41)

The function $B(r)$ is related to the quasi-local spatial curvature Sussman
and $M(r)$ is the gravitational mass contained within a comoving spherical
shell at a given $r$. Equation (40) can be integrated to yield

$\int\limits_{0}^{R}{{{d\tilde{R}}\over{\sqrt{2B+{{2M}\over{\tilde{R}}}+{1\over
3}\Lambda\tilde{R}^{2}}}}}=t-t_{B}(r).$ (42)

Here $t_{B}(r)$ is a third free function of $r$ (called the bang time
function). In the LTB model, in general, the Big Bang is not simultaneous as
in the FRW case, but it depends on the radial coordinate $r$. The given
formulas are invariant under a transformation $\tilde{r}=g(r)$. We can use
this freedom to choose one of the functions $B(r),M(r)$ and $t_{B}(r)$. For
$\Lambda=0$ the above equation can be solved explicitly. The evolution can be
elliptic ($B<0$), parabolic ($B=0$) or hyperbolic ($B>0$).

##  6 Backreaction inside the LTB onion model

As an example of backreaction computation we consider an approximate LTB model
(the so-called onion model) investigated in Biswas  by Biswas, Mansouri and
Notari, who computed corrections to the luminosity distance–redshift relation.

Figure 1: Density profile at the time $t=20$.

It represents a spacetime with radial shells of overdense and underdense
regions. The function $B(r)$ is nonzero ($B(r)>0$), so the evolution of the
LTB model is hyperbolic. The metric function $R(t,r)$ reads

$R(t,r):=\left(\frac{6}{\pi}\right)^{1/3}t^{2/3}r\left(1+\left(\frac{81}{4000\pi^{2}}\right)^{1/3}\left(\frac{1}{2\pi}\right)t^{2/3}\frac{1}{r}\sin^{2}{\pi
r}\right).$ (43)

The function $B(r)$ is given as follows

$B(r)=\frac{r}{2\pi}\sin^{2}{\pi r}.$ (44)

The density profile at the time $t=20$ can be seen in Figure 1. The density is
computed using the formula (41). The coordinates were chosen so that the
function $M(r)$ is given by $M(r)=4/3\pi r^{3}$.

First, we investigate the backreaction term in Buchert’s equation (25). We
numerically integrate the underlined part of equation (25) depending on the
averaging scale $l$. As one can see from Figure 2, the backreaction normalized
by $\left\langle\theta\right\rangle^{\cdot}$ is negative. It has a peak for
the averaging scale $l\approx 0.8\,.$ The value of backreaction normalized by
$\left\langle\theta\right\rangle^{\cdot}$ is of the order of $10^{-4}$. The
backreaction term (without normalization) is positive and it leads to an
increase of expansion.

Figure 2: Backreaction term ${2\over
3}\left({\left\langle{\theta^{2}}\right\rangle-\left\langle\theta\right\rangle^{2}}\right)-2\left\langle{\sigma^{2}}\right\rangle$
in the evolution equation for expansion depending on the averaging scale $l$
and normalized by $\left\langle\theta\right\rangle^{\cdot}$.

We can investigate also the backreaction terms in other equations which do not
appear in the Buchert framework and which can supplement his equations. For
example here we will show the result for backreaction in the averaged
evolution equation for shear (26) (specifically the whole underlined part of
the equation is considered). All results depend on the averaging scale $l$. As
we can see from Figure 3 \- for small scales, the contribution of all
backreaction terms in the evolution equation for shear normalized by
$\left\langle\sigma\right\rangle^{\cdot}$ is negative with a peak around
$l\approx 0.9\,.$ For larger scales the contribution is smaller and positive.
The turning point is for $l\approx 1.2\,.$ To be more precise, there exist
regions where the backreaction changes the sign twice for a very small
increase of $l$. If we compare the backreaction with the time derivative of
the shear scalar, we can see that their ratio is of the order of
$10^{-4}-10^{-3}$. It means that the backreaction plays a more important role
in the averaged equation for shear than in the averaged equation for
expansion.

Note that we investigated only an approximate LTB model. Due to nonlinearity
it is not clear if the backreaction behavior shown above will be similar for
exact solutions (even when they are close to the onion model in some specific
sense). We used this non-exact model because it has suitable properties for
investigation of averaging and backreaction.

Figure 3: Backreaction term
${1\over{\sqrt{3}}}\left({\left\langle\sigma\right\rangle^{2}-\left\langle{\sigma^{2}}\right\rangle}\right)+{1\over
3}\left({\left\langle{\theta\sigma}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma\right\rangle}\right)$
in the evolution equation for shear depending on the averaging scale $l$ and
normalized by $\left\langle\sigma\right\rangle^{\cdot}$.

##  7 Averaged LRS dust class II equations

One of the most important equations in cosmology is the evolution equation for
the expansion scalar. In the averaged equation (25) we have independent
variables $\left\langle{\theta^{2}}\right\rangle$ and
$\left\langle{\sigma^{2}}\right\rangle$. To obtain evolution equation for
$\left\langle{\theta^{2}}\right\rangle$ we multiply (8) by $2\theta$. Then we
perform averaging and we obtain the equation

$\displaystyle\left\langle\theta^{2}\right\rangle^{\cdot}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}\left\langle\theta\right\rangle^{3}-4\left\langle\theta\right\rangle\left\langle\sigma\right\rangle^{2}-8\pi\left\langle\rho\right\rangle\left\langle\theta\right\rangle+\underline{\frac{1}{3}\left(\left\langle\theta^{3}\right\rangle-\left\langle\theta\right\rangle^{3}\right)}$
$\displaystyle\underline{+\left(\left\langle\theta\right\rangle^{3}-\left\langle\theta\right\rangle\left\langle\theta^{2}\right\rangle\right)-4\left(\left\langle\theta\sigma^{2}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma\right\rangle^{2}\right)-4\pi\left(\left\langle\rho\theta\right\rangle-\left\langle\rho\right\rangle\left\langle\theta\right\rangle\right)}.$

In a similar way we derive an evolution equation for
$\left\langle{\sigma^{2}}\right\rangle$

$\displaystyle\left\langle\sigma^{2}\right\rangle^{\cdot}$ $\displaystyle=$
$\displaystyle-\frac{2}{\sqrt{3}}\left\langle\sigma\right\rangle^{3}-\frac{4}{3}\left\langle\theta\right\rangle\left\langle\sigma\right\rangle^{2}-2\left\langle
E\right\rangle\left\langle\sigma\right\rangle+\underline{\frac{1}{3}\left(\left\langle\theta\right\rangle\left\langle\sigma\right\rangle^{2}-\left\langle\theta\sigma^{2}\right\rangle\right)}$
(46)
$\displaystyle\underline{-\frac{4}{3}\left(\left\langle\theta\right\rangle\left\langle\sigma^{2}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma\right\rangle^{2}\right)-\frac{2}{\sqrt{3}}\left(\left\langle\sigma^{3}\right\rangle-\left\langle\sigma\right\rangle^{3}\right)}$
$\displaystyle-\underline{2\left(\left\langle
E\sigma\right\rangle-\left\langle
E\right\rangle\left\langle\sigma\right\rangle\right)}.$

Now, we also need equations for e.g. $\left\langle{\theta^{3}}\right\rangle$,
$\left\langle{\theta\sigma^{2}}\right\rangle$ or
$\left\langle{\sigma^{3}}\right\rangle$ (and of course an evolution equation
for $\left\langle\rho\right\rangle$, $\left\langle E\right\rangle$ and
$\left\langle\sigma\right\rangle$ given in Section 4). We could obtain these
evolution equations by the same procedure. Thus we have an infinite number of
equations for the correlation terms. Here we need to adopt an ansatz. For
example we can consider a reasonable assumption that for a given order the
correlation terms are negligibly small and we can truncate the hierarchy to
obtain a finite set of equations. We can also assume that some terms are
proportional to each other. In this approach the inhomogeneities are modeled
by different relations for correlation functions.

The question is what kind of spacetime may correspond to the given set of
averaged equations. We have started with scalar equations characterizing a LRS
class II dust spacetime (containing the LTB metric and its generalizations).
By averaging we can not leave this class, instead we may end up in a special
subclass of LRS class II dust models. We performed averaging of the evolution
equation for expansion (25) and of the evolution equations for different
products of expansion, shear, density and electric part of the Weyl tensor.
From this construction we can see that the averaged equations contain an
averaged LTB model, but generally not e.g. the homogeneous LRS Bianchi
cosmologies.

In the above described approach we have evolution equations for averages of
different powers and products of the expansion, shear, density and electric
part of the Weyl tensor. The problem is that if we derive an evolution
equation for the averaged nonlinear terms, then more complicated terms appear
in the relevant correlation terms (as is evident in equations (7) and (46)).
If we want to close the system of equations we need to effectively eliminate
these higher-order terms. In general, we may express the ”unwanted terms” as a
suitable function of the lower-order averaged terms whose evolution equation
is known. This kind of ansatz does not need to make all higher-order
correlation terms necessarily vanish, it only serves to close the system of
averaged equations through the selected relations - these may be for example
expressed as products of the averaged terms of lower order.

##  8 Conclusion

We generalized Buchert’s equations for the LRS class II dust model. We used
the property that this family is characterized only by scalars and we employed
a similar technique for averaging. However, the system of averaged equations
is not closed. Buchert considered the so-called scaling solutions Buchert3  to
close the system of equations. In our work, we first investigated the
influence of backreaction on the expansion and shear scalars for an
approximate LTB model which describes fluctuating radial inhomogeneities. Then
we proposed an infinite system of equations which supplement the averaged
equations for expansion. In this approach inhomogeneities are modeled by the
form of the correlation terms. Finally, we discussed how to close the system
of averaged equations.

## Acknowledgments

We would like to thank R. Sussman for a useful discussion. We would also like
to express our gratitude to the referee for correcting and improving our
paper. P.K. was supported by grants GAUK 398911 and SVV-267301. O.S.
acknowledges the support of grant GAČR 14-37086G.

## Appendix

In this appendix we will show the computations demonstrating that the averaged
constraint equations (4) \- (4) are preserved in time. We start with an
unaveraged constraint (13) and perform straightforward averaging without using
the commutation rules which gives us

$\left\langle\sigma^{\prime}\right\rangle={1\over{\sqrt{3}}}\left\langle\theta^{\prime}\right\rangle-{2\over
3}\left\langle a\sigma\right\rangle.$ (A.1)

Now we take time derivative of (A.1) and use the commutation rule (17). We
obtain the following expression

$\displaystyle\left\langle(\sigma^{\prime})^{\cdot}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma^{\prime}\right\rangle+\left\langle\theta\sigma^{\prime}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{3}}\left[\left\langle(\theta^{\prime})^{\cdot}\right\rangle-\left\langle\theta\right\rangle\left\langle\theta^{\prime}\right\rangle+\left\langle\theta\theta^{\prime}\right\rangle\right]$
(A.2)
$\displaystyle-\frac{2}{3}\left[\left\langle(a\sigma)^{\cdot}\right\rangle-\left\langle\theta\right\rangle\left\langle
a\sigma\right\rangle+\left\langle\theta a\sigma\right\rangle\right].$

Now we need to commute prime and dot derivatives. This is done applying the
commutation rule (34)

$\displaystyle\left\langle(\dot{\sigma})^{\prime}\right\rangle-\frac{2}{\sqrt{3}}\left\langle\sigma\sigma^{\prime}\right\rangle-\frac{1}{3}\left\langle\theta\sigma^{\prime}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma^{\prime}\right\rangle+\left\langle\theta\sigma^{\prime}\right\rangle=$
(A.3)
$\displaystyle\frac{1}{\sqrt{3}}\left[\left\langle(\dot{\theta})^{\prime}\right\rangle-\frac{2}{\sqrt{3}}\left\langle\sigma\theta^{\prime}\right\rangle-\frac{1}{3}\left\langle\theta\theta^{\prime}\right\rangle-\left\langle\theta\right\rangle\left\langle\theta^{\prime}\right\rangle+\left\langle\theta\theta^{\prime}\right\rangle\right]-\frac{2}{3}\left[\left\langle(a\sigma)^{\cdot}\right\rangle-\left\langle\theta\right\rangle\left\langle
a\sigma\right\rangle+\left\langle\theta a\sigma\right\rangle\right].$

Next, we apply the unaveraged evolution equations for $\theta,\sigma$ and $a$
((8), (9) and (12)) and we obtain the following expression:

$\displaystyle-\frac{2}{\sqrt{3}}\left\langle\sigma\sigma^{\prime}\right\rangle-\frac{2}{3}\left\langle\theta^{\prime}\sigma\right\rangle-\frac{2}{3}\left\langle\theta\sigma^{\prime}\right\rangle-\left\langle
E^{\prime}\right\rangle-\frac{2}{\sqrt{3}}\left\langle\sigma\sigma^{\prime}\right\rangle-\frac{1}{3}\left\langle\theta\sigma^{\prime}\right\rangle-\left\langle\theta\right\rangle\left\langle\sigma^{\prime}\right\rangle+\left\langle\theta\sigma^{\prime}\right\rangle=$
$\displaystyle\frac{1}{\sqrt{3}}\left[-\frac{2}{3}\left\langle\theta\theta^{\prime}\right\rangle-4\left\langle\sigma\sigma^{\prime}\right\rangle-4\pi\left\langle\rho^{\prime}\right\rangle-\frac{2}{\sqrt{3}}\left\langle\sigma\theta^{\prime}\right\rangle-\frac{1}{3}\left\langle\theta\theta^{\prime}\right\rangle-\left\langle\theta\right\rangle\left\langle\theta^{\prime}\right\rangle+\left\langle\theta\theta^{\prime}\right\rangle\right]$
(A.4) $\displaystyle-\frac{2}{3}\left[-\frac{1}{3}\left\langle
a\theta\sigma\right\rangle+\frac{1}{\sqrt{3}}\left\langle
a\sigma^{2}\right\rangle-\frac{1}{\sqrt{3}}\left\langle
a\sigma^{2}\right\rangle-\frac{2}{3}\left\langle
a\theta\sigma\right\rangle-\left\langle
aE\right\rangle-\left\langle\theta\right\rangle\left\langle
a\sigma\right\rangle+\left\langle\theta a\sigma\right\rangle\right].$

From the above expression we can see that several terms cancel each other.
Moreover, using constraint equations (13) and (14) to further simplify the
above equation (Appendix) it is straightforward to see that the left hand side
is equal to the right hand side. This means that the constraint equation (4)
does not change in time. In the same way it can be shown that the constraint
equations (4) and (4) are preserved in time, too.

## References

  * (1) Zalaletdinov, R.M.,: Averaging out the Einstein equations, Gen. Rel. Grav. 24, 1015 (1992).
  * (2) Zalaletdinov, R.M.,: Towards a theory of macroscopic gravity, Gen. Rel. Grav. 25, 673 (1993).
  * (3) Korzyński, M.: Covariant coarse graining of inhomogeneous dust flow in general relativity, Class. Quantum Grav. 27, 105015 (2010).
  * (4) Brannlund, J., van den Hoogen, R., Coley, A.: Averaging geometrical objects on a differentiable manifold, Int.J.Mod.Phys. D19 1915-1923 (2010)
  * (5) Buchert, T.: On Average Properties of Inhomogeneous Fluids in General Relativity: Dust Cosmologies, Gen. Rel. Grav. 32, 105 (2000).
  * (6) Buchert, T.: On Average Properties of Inhomogeneous Fluids in General Relativity: Perfect Fluid Cosmologies , Gen. Rel. Grav. 33, 1381 (2001).
  * (7) Wiltshire, D.L.: Cosmic clocks, cosmic variance and cosmic averages , New J. Phys., 9, 377 (2007).
  * (8) Li, N., Schwarz, D.: Onset of cosmological backreaction, Phys. Rev. D 76, 083011 (2007).
  * (9) Li, N., Schwarz, D.: Scale dependence of cosmological backreaction, Phys. Rev. D, 78, 083531 (2008).
  * (10) Behrend, J., Brown, I. A., Robbers, G.: Cosmological backreaction from perturbations, JCAP 01, 013 (2008).
  * (11) Clarkson, Ch., Ananda, K., Larena, J.: Influence of structure formation on the cosmic expansion, Phys. Rev. D 80, 083525 (2009).
  * (12) Larena, J., Alimi, J. M., Buchert, T., Kunz, M., Corasaniti, P. S.: Testing backreaction effects with observations, Phys. Rev. D 79, 083011 (2009).
  * (13) Ellis, G.F.R.: Dynamics of Pressure-Free Matter in General Relativity, J. Math. Phys. 8, 1171 (1967).
  * (14) Stewart, J.M., Ellis, G.F.R.: Solutions of Einstein’s Equations for a Fluid Which Exhibit Local Rotational Symmetry , J. Math. Phys. 9, 1072 (1968).
  * (15) van Elst, H., Ellis, G.F.R.: The covariant approach to LRS perfect fluid spacetime geometries, Class. Quantum Grav. 13 1099 (1996).
  * (16) Sussman, R.: Back-reaction and effective acceleration in generic LTB dust models, Class. Quant. Grav. 28, 235002 (2011).
  * (17) Chuang, C. H., Gu, J. A., Hwang, W. Y. P.: Inhomogeneity-induced cosmic acceleration in a dust universe, Class. Quantum Grav. 25, 175001 (2008).
  * (18) Paranjape, A., Singh, T. P.: The possibility of cosmic acceleration via spatial averaging in Lemaître-Tolman-Bondi models, Class. Quantum Grav. 23, 6955 (2006).
  * (19) Räsänen S.: Backreaction in the Lemaître-Tolman-Bondi model, JCAP 0411, 010 (2004).
  * (20) Bolejko, K., Andersson, L.: Apparent and average accelerations of the Universe, JCAP 10, 003 (2008).
  * (21) Sussman, R.: On spatial volume averaging in Lemaître-Tolman-Bondi dust models. Part I: back reaction, spacial curvature and binding energy, arXiv:0807.1145 [gr-qc], (2008).
  * (22) Sussman, R.: Quasi-local variables and scalar averaging in LTB dust models, AIP Conf. Proc. 1241, 1146-55 (2010).
  * (23) Mattsson, T., Mattsson, M.: Exploiting scale dependence in cosmological averaging, JCAP02, 004 (2008).
  * (24) Mattsson, T., Mattsson, M.: On the role of shear in cosmological averaging, JCAP 1010, 021 (2010).
  * (25) Mattsson, T., Mattsson, M.: On the role of shear in cosmological averaging II: large voids, non-empty voids and a network of different voids, JCAP 1105, 003 (2011).
  * (26) Räsänen S.: Cosmological acceleration from structure formation , Int.J.Mod.Phys. D15, 2141-2146 (2006).
  * (27) Räsänen S.: Evaluating backreaction with the peak model of structure formation, JCAP 04, 026 (2008).
  * (28) Paranjape, A., Singh, T. P.: Structure formation, backreaction and weak gravitational fields, JCAP 03, 023 (2008).
  * (29) Lemaître, G.: L’Universe en Expansion, Ann. Soc. Sci. Bruxelles A53, 51-85 (1933).
  * (30) Tolman, R.C.: Effect of Inhomogeneity on Cosmological Models, Proc. Nat. Acad. Sci. U.S.A. 20, 169-76 (1934).
  * (31) Bondi, H.: Spherically Symmetric Models in General Relativity, Mon. Not. Roy. Astron. Soc. 107, 410 (1947).
  * (32) Bolejko, K., Krasiński, A., Hellaby, C., Célérier, M.N.: Structures in the Universe by exact methods: formation, evolution, interactions, Cambridge: Cambridge University Press, (2009).
  * (33) Hellaby, C.: Modelling Inhomogeneity in the Universe, 5th International School on Field Theory and Gravitation, Cuiab a, Brazil, 20-24 April 2009, Proc. Sci. PoS(ISFTG) 005 (2009).
  * (34) Biswas, T., Mansouri, R., Notari, A.: Non-linear structure formation and ’apparent’ acceleration: an investigation, JCAP 12, 0712 (2007).
  * (35) Buchert, T., Larena, J., Alimi, J.M.: Correspondence between kinematical backreaction and scalar field cosmologies - the ‘morphon field’, Class. Quant. Grav. 23, 6379 (2006). 320, 1 (1997).


# AFM insulating state and normal state of cuprates from first principle

Zheting Jin  Department of Applied Physics, Yale University, New Haven,
Connecticut 06511, USA  Sohrab Ismail-Beigi  Department of Applied Physics,
Yale University, New Haven, Connecticut 06511, USA  Department of Physics,
Yale University, New Haven, Connecticut 06511, USA  Department of Mechanical
Engineering and Materials Science, Yale University, New Haven, Connecticut
06511, USA

###### Abstract

Materials-realistic microscopic theoretical descriptions of copper-based
superconductors are challenging due to their complex crystal structures
combined with strong electron interactions. Here, we demonstrate how ab initio
calculations can accurately describe key structural, electronic, and magnetic
properties of the normal state of the prototypical cuprate Bi2Sr2CaCu2O8+x
(Bi-2212). We emphasize the importance of accounting for energy-lowering
structural distortions, which then allows us to: (a) accurately describe the
insulating antiferromagnetic (AFM) ground state of the undoped parent compound
(in contrast to the metallic state predicted by previous ab initio studies);
(b) identify numerous low-energy competing spin and charge stripe orders in
the hole-overdoped material that are nearly degenerate in energy with the
G-AFM ordered state, indicating strong spin fluctuations; (c) predict the
lowest-energy hole-doped crystal structure which exhibits long-range
structural distortions and oxygen dopant positions matching high-resolution
scanning transmission electron microscopy (STEM) measurements; and (d)
describe electronic bands near the Fermi energy with flat antinodal
dispersions and Fermi surfaces that in agreement with angle-resolved
photoemission spectroscopy (ARPES) measurements and provide a clear
explanation for the structural origin of the so-called “shadow bands”.

##  I Introduction

The cuprate superconductors continue to be a fascinating and actively
researched class of materials. In addition to their superconducting phases,
their normal state has attracted broad research interest due to a wide range
of unusual properties. Understanding the physical origin of the AFM insulating
phase Kastner et al. (1998); Ho et al. (2001), the pseudogap Warren Jr et al.
(1989); Homes et al. (1993); Loeser et al. (1996); Ding et al. (1996), the
flat antinodal dispersion Andersen et al. (1994); Liechtenstein et al. (1996);
He et al. (2021); Li et al. (2021), the strange metallicity Gurvitch and Fiory
(1987); Martin et al. (1990); Cooper et al. (2009); Greene et al. (2020);
Ayres et al. (2021), and the presence of quantum critical fluctuations Dai et
al. (2001); Le Tacon et al. (2011); Dean et al. (2013); Keimer et al. (2015)
in the normal state can provide important insights into the underlying
mechanisms that can give rise to the superconductivity. Despite plenty of
proposed possible mechanisms, including competing orders Kivelson et al.
(1998); Lee et al. (2006); Chang et al. (2012) and preformed pairs Emery and
Kivelson (1995); Corson et al. (1999); Li et al. (2010); Dubroka et al.
(2011), the physical origin of this complicated normal state is still unclear.

A comprehensive description of the normal state is challenging due to the
strong electronic interactions combined with the structural complexity of
typical doped cuprates. While the electronic correlations can be captured by
effective model Hamiltonian approaches using accurate methods such as density
matrix renormalization group (DMRG) White (1992); Schollwöck (2005) or quantum
Monte Carlo (QMC) Ceperley and Alder (1986); Foulkes et al. (2001), the
effects of the complicated lattice distortions introduced by dopants and
impurities necessitate a realistic and detailed understanding of the
materials. These structural distortions can result in significant changes to
the materials such as additional symmetry breaking Gao et al. (1988); Mans et
al. (2006) and modified superconducting temperatures Andersen et al. (2007).

Density functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965)
offers a potent foundational method for investigating the ground-state
properties of materials from first principles. For cuprates, DFT has played a
pivotal role in constructing effective model Hamiltonians Hybertsen et al.
(1989); Martin (1996); Markiewicz et al. (2005); Korshunov et al. (2005);
Delannoy et al. (2009); Hanke et al. (2010). To ensure the accuracy of these
models, it is essential that the DFT calculations capture correctly both the
structural properties and the predominant electronic properties of the ground
state. Indeed, recent DFT studies on transition metal oxides have highlighted
the significance of allowing energy-lowering structural distortions to achieve
high-quality predictions of materials properties Varignon et al. (2019); Wang
et al. (2021); Zunger (2022).

Bismuth strontium calcium copper oxide Bi2Sr2CaCu2O8+x (BSCCO or Bi-2212)
Maeda et al. (1988) is one of the most intensively studied cuprates and is the
focus of in this work. However, prior DFT studies faced challenges due to its
intricate structural distortions and superlattice modulations Petricek et al.
(1990); Levin et al. (1994); Slezak et al. (2008); He et al. (2008). We
demonstrate that by employing modern exchange-correlation functionals in DFT
and providing an accurate description of the crystalline structure including
energy-lowering lattice distortions, we can directly describe the
antiferromagnetic insulating ground state of BSCCO, its correct crystal
structure upon doping, the presence of competing magnetic and charge stripe
orders, as well as crucial details of photoemission spectra such as “shadow
bands”. This means that we can, from first principles, correctly ascribe
certain experimental observations to specific structural motifs (e.g., shadow
bands) while simultaneously helping build microscopically well-justified model
Hamiltonians that can be used to compute the effects of strong electron
correlations.

As a first example, we note that to date, the ground state properties of
Bi-2212 have not been described well by DFT, which hampers further theoretical
studies of the effects of doping and other perturbations on the system at a
microscopic level. For instance, an important problem involves the undoped
parent material, which serves as the starting point for any systematic
analysis. Prior works were unable to reproduce the AFM insulating behavior of
undoped Bi-2212 Wang et al. (2005); Lin et al. (2006); Foyevtsova et al.
(2010), even when considering the crystal modulation correction He et al.
(2008); Fan and Zeng (2011). Here, we provide a realistic description of
Bi-2212 by capturing the dominant lattice and electronic properties of the
normal state using DFT. For the undoped system, we demonstrate insulating
bismuth-derived bands as a direct result of strong Bi-O bonding in the
optimized and distorted BiO layers which explains the absence of bismuth-
derived Fermi surface in all experiments regardless of doping level. In
contrast, previous DFT studies Lin et al. (2006); He et al. (2006); Foyevtsova
et al. (2010) had to resort to additional “manual” hole doping to eliminate
the finite bismuth density of states at the Fermi level that plagued the
predictions for undoped Bi-2212.

Second, armed with a proper description of the undoped material, our doped
Bi-2212 results predict correct structural modulations and realistic oxygen
dopant positions that match STEM measurements Song et al. (2019). Third, using
this optimized crystal structure, we discover many stripe-ordered electronic
states that are almost degenerate with the G-type (checkerboard) AFM-ordered
state. This indicates the presence of strong spin and charge fluctuations and
competing orders in the normal state and might be a possible physical origin
of the pseudogap phase in cuprates Zhang et al. (2020).

Fourth, we compute the electronic band structure near the Fermi energy and
reproduce the flat band dispersions around the antinodal regions of the
Brillouin zone, as well as details of the Fermi surface including the shadow
bands observed via ARPES He et al. (2021). Our analysis identifies two
distinct structural origins for different types of shadow bands. By capturing
the dominant properties of the crystal and electronic structure in Bi-2212
under a first principle framework, our work can serve as a useful platform for
future theoretical modeling and possible engineering of these superconductors.

##  II Undoped system

Figure 1:  Crystal and electronic structure of high-symmetry undoped Bi-2212.
(a) Side view of the crystal structure along the $b$-axis. Sr, Ca, Cu, Bi, and
O atoms are marked by green, gray, blue, purple, and red balls, respectively.
The black square in the crystal structure marks the G-AFM unit cell. (b) Top
view of the BiO layer from the $c$-axis. (c) Projected band structures (left)
and density of states (DOS) (right) from DFT calculation, where the Fermi
energy is set to be the reference energy. The inset of the band structures
shows the first Brillouin zones (BZ) of the 60-atom supercell. The unit of DOS
is the number of states per unit cell (UC) per eV. Red circles and solid lines
show Bi $p$-orbitals; blue squares and dash lines show Cu $d$-orbitals; black
dash-dotted line shows O $p$-orbitals in BiO layers.

The undoped Bi-2212 possesses a bilayer crystal structure as depicted in Fig.
1(a)(b), where each bilayer comprises two CuO layers separated by one Ca layer
and sandwiched between SrO and BiO layers. The smallest unit cell contains 30
atoms and crystallizes in the tetragonal I4/mmm space group, with only one Cu
atom in each CuO layer Liang et al. (1988); Jain et al. (2013); mat . However,
this small unit cell leads to a false non-magnetic ground state due to the
artificial assumption of translational invariance of the Cu local moments
Yelpo et al. (2021). Here, we focus on a larger supercell (60 atoms/cell) with
two Cu atoms in each CuO layer, which allows for the spontaneous symmetry
breaking of the magnetic local moments. First, we perform a relaxation of the
database crystal structure mat ; Jain et al. (2013); Subramanian et al. (1988)
to identify the nearest local minimum with the same symmetry, referred to as
the “high-symmetry” structure. As a result, we find a G-AFM ordered ground
state with local Cu magnetic moments of $\pm 0.45\mu_{B}$, while the non-
magnetic state is about 0.25 eV/Cu higher in energy. These local moments agree
with the experimental measurements, typically falling within the range of
0.4-0.6$\mu_{B}$ in cuprates without chlorine Schrieffer and Brooks (2007).
Fig. 1(c) shows the projected band structure and density of states of this AFM
ground state. The AFM order opens a gap of about 0.6 eV for Cu $d$-orbitals
between $M^{\prime}$ and $X^{\prime}$. However, this calculation still results
in a false metallic ground state due to the widely dispersive Bi $p$-bands
located within the AFM gap.

Figure 2:  Crystal and electronic structure of two stable low-symmetry
crystals of undoped Bi-2212. (a)(b) The top views of the “zigzag” and
“orthorhombic” distortion patterns of the BiO layer, where the orthorhombic
distortion pattern in (b) is the most energetically favorable structure. Large
purple and small red balls represent Bi and O atoms. The black squares in the
crystal structure illustrate the G-AFM unit cells. (c)(d) Projected band
structures (left) and the density of states (right) of the crystal structures
in (a) and (b), separately. Red circles and solid lines show Bi $p$-orbitals;
blue squares and dash lines show Cu $d$-orbitals; black dash-dotted line shows
O $p$-orbitals in BiO layers.

To improve the theoretical description, we find that allowing crystal
structural distortions in the BiO layers can further lower the total energy of
the calculation. Using the conjugate gradient algorithm for structural
relaxation, we optimize the structures by initially lowering the symmetry in
the BiO layers manually. Fig. 2(a)(b) depicts two typical stable BiO layer
patterns whose ground states are both G-AFM ordered on the CuO2 planes. Fig.
2(a) shows a zigzag pattern, where all the oxygen atoms on BiO layers are
moved along the diagonal direction of the in-plane unit cell by the same
amount. With all BiO layers displaying this zigzag distortion pattern, the AFM
ground state energy is about 0.29 eV/Cu lower than the one with the high-
symmetry structure in Fig. 1(b), and the local Cu magnetic moments are $\pm
0.48\mu_{B}$, slightly larger than the high-symmetry case. Fig. 2(b) shows an
orthorhombic pattern, characterized by the lowest symmetry among all three
structures presented in Fig. 1 and 2. This pattern further breaks the mirror
symmetry of the crystal, resulting in an orthorhombic lattice consistent with
observations in experiments Zeljkovic et al. (2012). This orthorhombic pattern
exhibits the lowest AFM ground state energy, 0.5 eV/Cu lower than the high-
symmetry case. It also features the largest local Cu magnetic moments of $\pm
0.53\mu_{B}$ among all three Bi-O structural motifs mentioned above.
Consequently, this structure displays the largest AFM gap of about 1 eV,
consistent with the gap size observed in a recent scanning tunneling
microscopy (STM) experiment Wang et al. (2023).

The electronic structures of the distorted crystals in Figs. 2(a) and (b) are
presented in Figs. 2(c) and (d), respectively. Compared to the electronic
structure of the high-symmetry case in Fig. 1(c), most of the changes are
observed within the Bi bands due to the BiO distortion patterns. Although the
zigzag distortion pattern helps reduce the size of the Bi electron pocket,
it’s only the lowest-energy orthorhombic distortion pattern that elevates the
entire Bi bands above the Fermi level, resulting in an insulating AFM ground
state.

A straightforward microscopic picture helps elucidate how the distortions
contribute to elevating the Bi bands to higher energy levels. The densities of
states (DOS) plots reveal that the in-plane BiO system possesses filled low-
energy bonding states dominated by oxygen and antibonding states dominated by
bismuth, and the bonding/anti-bonding gaps are centered at about an energy 1
eV below the Fermi energy for all three structures. However, the size of the
bonding/anti-bonding gap varies with the distortion pattern. The DOS shows
that this gap is smallest for the high-symmetry structure in Fig. 1(c)
($\sim$1 eV) and largest for the orthorhombic structure in Fig. 2(d) ($\sim$2
eV). The gap size difference arises from different Bi-O hybridization
strengths among the three structures. The coupling is strong enough only in
the orthorhombic structure to lift the anti-bonding state above the Fermi
level, while the couplings in high-symmetry and zigzag structures are too
weak, resulting in metallic states.

This microscopic picture is consistent with a structural analysis: the Bi-O
bond lengths in the high-symmetry, zigzag, and orthorhombic structures are
2.65, 2.28, and 2.17Å, respectively. In general, a shorter bond length implies
a larger hopping element between two atoms. We substantiate this picture
quantitatively by computing the tight-binding Kohn-Sham Hamiltonian on the
maximally localized Wannier basis Marzari and Vanderbilt (1997) extracted from
our DFT calculations using Wannier90 Pizzi et al. (2020). The projected
Wannier orbitals encompass the $p$-orbitals of Bi and O, as well as the
$d$-orbitals of Cu, which are sufficient to describe the bands near the Fermi
level as shown in supplementary Fig. S LABEL:supp-fig:wannier90_vs_vasp Jin
and Ismail-Beigi . Our calculations reveal that the Bi-O hopping (tunneling)
matrix element for nearest-neighbor in-plane Bi-O pairs is about 1.0, 1.8, and
2.3 eV for the high-symmetry, zigzag, and orthorhombic structures,
respectively, which is consistent with the behavior of the gaps. Notably, the
hopping strength in the orthorhombic structure is the only one larger than 2
eV, a value sufficient to open a gap and raise the Bi electron pockets
(originating from the anti-bonding BiO bands) above the Fermi level, which
leads to the insulating ground state.

##  III Oxygen-doped system

Figure 3:  Crystal structure of Bi2Sr2CaCu2O8.25, where the hole doping level
is $x=0.25$. Sr, Ca, Cu, Bi, O, and dopant O atoms are marked by green, gray,
blue, purple, red, and black balls, respectively. Red arrows further highlight
the O dopants. The black solid square in the crystal structure marks the
244-atom unit cell of the doped crystal. The black dashed rectangle highlights
the CuO layers in one of the bilayers in the unit cell.

Commencing with undoped Bi-2212, hole doping is incorporated through
interstitial oxygen dopants within the material. These additional hole dopants
give rise to a diverse range of physical properties including
superconductivity, the pseudogap phenomenon, and the presence of shadow bands
within the Fermi surface. The subsections below present our findings on
various aspects: 1. The crystal structure including the long-range
superlattice modulation and associated placement of the oxygen dopants. 2. The
low-energy magnetic states including spin- and charge-stripes. 3. Detailed
ARPES spectra around the Fermi level, which allows us to understand the role
of structural distortions in creating the shadow bands.

###  III.1 Crystal structure

For the oxygen-doped Bi2Sr2CaCu$2$O8+x system with $x=0.25$, we account for
the superlattice modulation Gao et al. (1988); Petricek et al. (1990); Levin
et al. (1994); Slezak et al. (2008); He et al. (2008) and explore various
possible positions for the oxygen dopants. These structural properties are
essential for describing the BiO layers properly and have a significant impact
on the superconducting gap Andersen et al. (2007); Massee et al. (2020).

Fig. 3 shows our optimized lowest-energy crystal structure for
Bi2Sr2CaCu2O8.25, where four oxygen dopants are added into a stoichiometric
undoped Bi-2212 unit cell of 240 atoms which is a $4\times 1\times 1$
enlargement of the 60-atom unit cell of Sec. II. Several metastable crystals
are listed in the supplementary material Jin and Ismail-Beigi , showing an
energy cost of at least 0.9 eV/dopant to move oxygen dopants between two BiO
layers, and 168 meV/dopant to move dopants between BiO and SrO layers. The
optimized super modulation periodicity is 8 atom sites as enforced by the
$x=0.25$ doping level. The most stable positions for the oxygen dopants are
located at the necking region between BiO and SrO layers, all in good
agreement with a recent high-resolution STEM study Song et al. (2019) on a
similar doping level.

Figure 4:  Competing G-AFM and stripe order phases in a CuO bilayer. (a,b)
Illustrations of G-AFM (a) and a typical stripe order phase (b) for the
bilayer in the black dashed square of Fig. 3. The black arrows represent the
local moments on Cu atoms, with exaggerated length and thickness highlighting
their relative magnitudes. Yellow dashed lines show the magnetic domain walls
in the stripe order phase. (c) The magnitude of local moments in $\mu_{B}$ on
the Cu sites $i$ along $a$-axis defined as
$|m_{i}|\equiv|n_{i\uparrow}-n_{i\downarrow}|\mu_{B}$, where $n_{i\sigma}$ is
the occupancy of the $d_{x^{2}-y^{2}}$ Wannier orbital at site $i$ with spin
$\sigma$. (d) The $d_{x^{2}-y^{2}}$ occupancy on Cu atomic sites defined as
$|n_{i}|\equiv|n_{i\uparrow}+n_{i\downarrow}|$. Red dashed lines and blue
solid lines represent the GAFM and stripe-order states, respectively. Circles
and diamonds show the results of the upper and lower layers.

###  III.2 Magnetic and charge ordering

With the lattice structure of our doped BSCCO confirmed, we turn our attention
to the spin structure within this system. We initiate with short-period
magnetic orderings of the Cu magnetic moments. Not surprisingly, the most
energetically favorable magnetic order among them is the G-AFM order with
antiparallel nearest-neighbor spins on Cu atoms as illustrated in Fig. 4(a).
Other short-period meta-stable magnetic orders exhibit aligned nearest-
neighbor spins, either in an intralayer ($\sim 40$meV/Cu higher in energy) or
interlayer ($\sim 2$meV/Cu higher in energy) fashion Jin and Ismail-Beigi .
Separately, the energy cost associated with changing the inter-bilayer spin
alignment from AFM to FM is negligible Jin and Ismail-Beigi . Hence, in the
subsequent discussions, we concentrate the discussion on the magnetic
structure within a single bilayer.

Below, we will delve into longer-period spin and charge orders, making it
crucial to quantify the local Cu moment and electron count in a precise
manner. A widely adopted and chemically intuitive approach involves assessing
the oxidation state of Cu. For example, an undoped Cu ion would have the
formal oxidation state of Cu2+, and each added hole on that Cu would raise its
oxidation state by one unit. Within band theory, this corresponds directly to
band occupancy: we count the occupancy of bands that are labeled to “belong”
to particular atomic orbitals. In our case, the band structures reveal that
the bands at the Fermi energy primarily exhibit Cu $d_{x^{2}-y^{2}}$
character, as demonstrated in supplement Fig.S LABEL:supp-fig:supp_244_wannier
and LABEL:supp-fig:supp_244_wannier_GAFM Jin and Ismail-Beigi . Consequently,
we construct a Wannier basis using only orbitals of Cu $d_{x^{2}-y^{2}}$
character via the maximally localized Wannier function method. As shown in
supplemental material Sec. LABEL:supp-sec:wannier Jin and Ismail-Beigi , this
Wannierized Hamiltonian reproduces the DFT band structure of Cu
$d_{x^{2}-y^{2}}$ derived bands about the Fermi level. In essence, we have
established a localized basis that replicates the ab initio band structure
near the Fermi energy, constructed with one orbital per Cu of
$d_{x^{2}-y^{2}}$ character. Subsequently, we employ this tight-binding
representation to compute the band structure, band occupancies, and local
occupancies of the Wannier orbitals. For a thermal Fermi-Dirac smearing of 100
K, Figure 4(c) presents the local moment magnitudes for the G-AFM state. The
moments are around 0.47$\mu_{B}$ with small modulations of about $\pm
0.04\mu_{B}$ due to the superlattice modulation. A similar modulation of $\pm
0.03e$ also manifests in the $d_{x^{2}-y^{2}}$ occupancies, as shown in Fig.
4(d).

Beyond short-period ordering, we have discovered numerous longer-ranged
stripe-ordered states that exhibit nearly degenerate energies with the G-AFM
state. While DFT-based stripe-ordered states were previously reported in LSCO
Anisimov et al. (2004); Pesant and Côté (2011) and YBCO Zhang et al. (2020),
they have not been observed in BSCCO up to now. For Bi-2212, we present a
typical bond-centered stripe-ordered state in Fig. 4(b), where the nearest
neighbor spins crossing the dashed domain walls align in parallel, in contrast
to the antiparallel alignment in the G-AFM state. Remarkably, the total energy
of this stripe-order state is only 1.9 meV/Cu higher than the G-AFM state. In
the supplementary materials Jin and Ismail-Beigi , we tabulate eight distinct
stripe order patterns with energy costs below 3 meV/Cu. The existence of all
these low-energy orders suggests the presence of strong spin fluctuations in
the normal state, which can play an important role in superconducting pairing
Moriya (2006); Dahm et al. (2009). Fig. 4(c) and (d) show the local moment and
$d_{x^{2}-y^{2}}$ occupancy of this stripe-order state. Notably, the sites
closer to the domain walls exhibit lower occupancies and smaller local moments
than the other sites. These striped spin and charge orders present a
modulation of $\pm 0.2\mu_{B}$ for local moments and $\pm 0.15$ for
occupancies, significantly larger than the modulation caused by structural
supermodulation in the G-AFM state. Consequently, the formation of the stripe
order has an electronic origin, but the precise location of the domain
boundaries can be influenced by the superlattice modulation effect.

Notice that the modulations of local moments and the electron occupancies are
almost the same numerically. This is because we have one orbital per site, and
within band theory, the mechanism to form a local moment is to have an
exchange splitting with an occupied low-energy spin-majority orbital and empty
high-energy minority-spin orbital: the doped holes go into the spin-majority
orbital and reduce both local occupancy and magnetic moment simultaneously.
Hence, the redistribution of doped holes strictly follows the change of spin
structures.

In addition to the Wannier basis described above, we have also analyzed the
electronic structure of the stripe orders with the standard atomic projections
output by the Vienna ab initio simulation package (VASP) software Kresse and
Furthmüller (1996a, b) (Fig.S LABEL:supp-fig:supp_stripe_order_VASP of the
supplementary material Jin and Ismail-Beigi ). The results are difficult to
interpret due to the non-orthogonality of the standard VASP projections as
explained in the supplement: the Cu magnetic moments show spatial modulations
similar to those in Fig. 4 but the total occupancies hardly vary from site to
site; additionally, there is a large oxygen contribution to the bands crossing
the Fermi level. In short, in contrast to our $d_{x^{2}-y^{2}}$ Wannier basis
and the intuitive picture it provides, the Cu VASP projections are
insufficient to explain the stripe order and we do not discuss them further.

Figure 5:  Electronic structures of the 25% hole-doped Bi-2212 system with a
non-magnetic state. (a) The unfolded band structure of Cu $d_{x^{2}-y^{2}}$
orbitals for the 244-atom unit cell. The opacity represents the spectral
weight. (b) Schematic Fermi surfaces are black curves. The black, yellow,
blue, and red rectangles represent the first BZ of the 1, 2, 4, and 8 Cu per
layer unit cell, respectively. (c-e) Unfolded Fermi surfaces for the high-
symmetry crystal of Fig. 1(b), the Bi-O distorted crystal of Fig. 2(b), and
the hole-doped crystal of Fig. 3, respectively. The Fermi level in (c) and (d)
is shifted by $-0.1$eV to allow a fair comparison to the hole-doped Fermi
surface in (e). Yellow dashed ($q_{1}$) and red solid ($q_{2}$) arrows show
two different coupling wave vectors. (f) Measured ARPES Fermi surface at
$T=104$ K, adapted from experiments He et al. (2021). Solid and dashed black
curves highlight two distinct kinds of Fermi surface curves.

###  III.3 ARPES spectra

The topology of the Fermi surface and the associated low-energy electronic
spectrum, usually measured by ARPES Damascelli et al. (2003), provides
important insights into the electronic properties of solids. In particular,
many materials exhibit a so-called “shadow band” (SB) Fermi surface,
resembling a weak-intensity copy of their main band (MB) Fermi surface with
certain shifted vectors in momentum space. Depending on the system, these SB
Fermi surfaces can, in principle, originate from any type of symmetry breaking
such as electronic Brouet et al. (2008), magnetic Rotenberg et al. (2005), or
structural Mans et al. (2006) origins. The physical origin of the SB Fermi
surface is crucial for understanding the physical properties of the material,
but it is difficult to distinguish between the different possible origins from
ARPES alone.

The Fermi surface of Bi-2212, as revealed by intensive ARPES studies, includes
weaker intensity SB in addition to the main bands (see Fig. 5(f), adapted from
a previous ARPES experiment He et al. (2021)). These SB can be described by
two types of symmetry-breaking vectors. One of them always aligns with the
superstructural modulation direction, while the other is along
$\pm(\pi,\pm\pi)$, coincident with the AFM ordering vectors. This has led to a
continued debate regarding the magnetic Kampf and Schrieffer (1990); Langer et
al. (1995); Chubukov (1995) versus the structural origin Koitzsch et al.
(2004); Mans et al. (2006); Nakayama et al. (2006); Yu et al. (2020) of the
$\pm(\pi,\pm\pi)$ folding vector. Theoretical interpretation of these SB Fermi
surfaces from first principles has been lacking. Below, we meticulously
predict the ARPES Fermi surface using our first-principle calculations and
delve into the underlying physical mechanisms behind the emergence of SB. We
will demonstrate that the emergence of these two distinct types of SB is
attributed to two distinct structural symmetry-breaking mechanisms.

Due to the strong spin fluctuations indicated by the competing strip orders
and G-AFM states, the normal state of the hole-doped system cannot be
described by a single magnetic-ordered configuration (i.e., a single Slater
determinant). In principle, an accurate account of quantum spin fluctuations
is needed for a comprehensive theoretical description including the magnetic
susceptibility Stepanov et al. (2018) and pseudogap Zhu et al. (2023)
properties in cuprates. In practice, however, computed spectra in a non-
magnetic state of the cuprates are usually (and surprisingly) comparable to
ARPES measurements Damascelli et al. (2003); Sobota et al. (2021).

In this context, we solely consider the electronic band structure of the non-
magnetic state, while details about the G-AFM and stripe-ordered states are
provided in the supplementary materials Jin and Ismail-Beigi . As an added
benefit, by removing the spin degree of freedom allows us to focus exclusively
on symmetry breaking through structural perturbations. As demonstrated below,
the non-magnetic state successfully accounts for many normal state
spectroscopic properties, including the SB. To facilitate comparison to
experimental ARPES spectra of the Fermi surface, we employ a standard “band-
unfolding” method Ku et al. (2010); Brouet et al. (2012) to project the band
structure onto the primitive unit cell Brillouin zone. This approach is known
to reproduce qualitative spectrum intensities observed in various ARPES
experiments on different materials Lin et al. (2011); Medeiros et al. (2014);
Tomić et al. (2014); Zhu et al. (2018). The band structures and unfoldings
discussed below are computed using the Cu $d_{x^{2}-y^{2}}$ Wannier basis
described in Sec. III.2.

The panels of Fig. 5 display the unfolded band structures. Fig. 5(a) shows the
case of the 244-atom 25% hole-doped system for the non-magnetic phase. Around
the antinodal region near M/X’ in Fig. 5(a), we find two branches of flat
bands below the Fermi level, which are split by the interlayer coupling in the
bilayers, consistent with the ARPES experiments He et al. (2021, 2018).

Fig. 5(e) exhibits the unfolded Fermi surface of this 25% hole-doped system.
The curves with the highest intensity contain two easily visible curves
corresponding to the solid and dashed black curves of Fig. 5(f) from the ARPES
measurement. Aside from these main curves, the ARPES Fermi surface exhibits a
complicated set of shadow bands with lower intensities which are also
reproduced in Fig. 5(e). We will now conduct an analysis to demonstrate that
the different sets of shadow bands on the Fermi surface have two distinct
physical origins.

Our analysis begins with the examination of the high-symmetry undoped crystal
shown in Fig. 1(b). Although the crystal contains 2 Cu atoms in each CuO
layer, the non-magnetic state of this crystal can be described in a smaller
primitive cell with one Cu atom in each CuO layer, due to full translational
symmetry. Consequently, the unfolded Fermi surface of this crystal shown in
Fig. 5(c) only contains the main Fermi curves without any shadow bands, as
expected. (The Fermi levels of the undoped crystals are intentionally lowered
by $-0.1$ eV to simulate hole doping.) Starting from this clean Fermi surface,
we successively add complexity to the crystal to see the emergence of the
shadow bands.

The first type of shadow band exhibits an approximately circular shape as
shown in Fig. 5(d). This circular shape arises from structural distortions in
the BiO layers of the undoped crystal, as depicted in Fig. 2(b). Previous
experiments have observed these shadow bands and suggested their likely
structural origin Koitzsch et al. (2004); Mans et al. (2006); Nakayama et al.
(2006) without providing a specific microscopic picture. Here, we
microscopically observe that the distortions in BiO layers break the
translational symmetry and introduce an inter-band coupling at wave vectors
$q_{1}=\pm(\pi,\pm\pi)$. This coupling leads to folding from the main bright
curves to the shadow bands as indicated by the yellow dashed arrows in Fig.
5(d) and Fig. 5(e).

On top of the circle-like shadow band, there is another type of shadow band
involving $q_{2}=\pm(\pi/4,-\pi/4)$ folding vectors shown in Fig. 5(e). This
type of shadow band, often referred to as a “superstructure” in previous
studies Liu et al. (2019), has long been attributed to a post-emission
modulation effect from BiO layer buckling Yu et al. (2020). Prior experiments
have also shown that this type of shadow band fades away when the crystal
modulation is gradually reduced by Pb doping Liu et al. (2019). Our
calculation presents a consistent microscopic picture: the oxygen dopants and
associated superlattice modulations further break the translational symmetry
and introduce an additional coupling at $q_{2}$. This coupling creates a
further folding (illustrated by the red arrow) from the main bright curves to
another set of shadow bands. It is important to note that this coupling only
occurs in one crystalline direction due to the superlattice modulations being
solely along the $a$-axis as shown in Fig. 3.

##  IV Outlook

In conclusion, within a first principle framework, we have provided a
microscopic understanding of the AFM insulating phase of the undoped Bi-2212
system. Additionally, we have uncovered competing stripe orders in the hole-
overdoped system, which offers a paradigmatic approach (specifically for
BSCCO) to describe stripe orders using DFT without performing complex many-
body calculations. Spectroscopically, our non-magnetic DFT band structure
calculations remarkably reproduce the observed normal state spectral
properties. Furthermore, we have elucidated the structural origin of the ARPES
$\pm(\pi,\pm\pi)$ and $\pm(\pi/4,-\pi/4)$ shadow bands in the hole-overdoped
system. Our work underscores the importance of considering the crystal degrees
of freedom, including structural distortions, modulations, and realistic
oxygen dopant positions, for an accurate description of various material
properties.

We believe that our study establishes a robust first-principle foundation for
more “surgical” structural engineering of cuprates, particularly with regard
to manipulating broken translational and rotational symmetries. Moreover, as
our DFT ground state captures the dominant low-energy properties in the normal
state of Bi-2212, the Wannierized Hamiltonians extracted from our DFT
calculations can provide appropriate theoretical model Hamiltonians that can
serve as platforms for future applications and research. An intriguing
direction for future work involves a deeper exploration of why naive non-
magnetic band structure calculations describe the observed ARPES spectra so
well when the actual materials exhibit numerous nearly degenerate competing,
but clearly magnetic, low-energy configurations.

##  V Computational Details

We use the Vienna ab initio simulation package (VASP) with the projector-
augmented wave method Kresse and Joubert (1999). A relatively high plane-wave
cutoff energy of 500 eV is used. The generalized-gradient-approximation (GGA)
with the semilocal Perdew–Burke–Ernzerhof (PBE) functional Anisimov et al.
(1991); Perdew et al. (1996) is used in all of our calculations. All
calculations are done with collinear spins. The recent finding of non-
collinear spin texture in cuprates Gotlieb et al. (2018) is also an
interesting topic to study but is beyond the scope of this work. To avoid the
known failure of the local-spin-density approximation (LSDA) and the GGA in
reproducing the copper magnetic moment in cuprates Mattheiss and Hamann
(1988); Hybertsen and Mattheiss (1988); Lin et al. (2006); Foyevtsova et al.
(2010); Song et al. (2019); He et al. (2008); Fan and Zeng (2011) due to self-
interaction errors (SIE) in the approximate exchange-correlation functionals
Perdew and Zunger (1981), we add $U=4$ eV for the Cu $3d$ manifold in all our
PBE+$U$ calculations following previous theoretical works Yelpo et al. (2021);
Wang et al. (2006); Deng et al. (2019). Although we treat $U$ as an adjustable
parameter, changing its value only shifts the energy of the unoccupied high-
energy Cu-derived bands Jin and Ismail-Beigi  which does not affect our main
findings around the Fermi energy. On a related note, a recent study Nokelainen
et al. (2020) based on the strongly-constrained-and-appropriately-normed
(SCAN) meta-GGA exchange-correlation functional Sun et al. (2015) predicted
the correct AFM magnetic ground state for undoped Bi-2212 but failed to
reproduce its insulating nature. This is not surprising because the SCAN
functional still displays underestimated band gaps and fails to reproduce the
insulating ground state in a wide range of transition metal oxides Gautam and
Carter (2018) including CuO Long et al. (2020). In the supplement Jin and
Ismail-Beigi , we show that the SCAN+U method Gautam and Carter (2018); Long
et al. (2020) provides very similar results to the PBE+U results we show
below. Details on optimized structures are found in the supplementary
materials Jin and Ismail-Beigi .

To accelerate the structural relaxation of the large 244-atom oxygen-doped
unit cell, we use the Spanish Initiative for Electronic Simulations with
Thousands of Atoms (SIESTA) package Soler et al. (2002) to approximately relax
the structure before doing the final relaxations using VASP. A DZP basis with
EnergyShift of 100 meV and SplitNorm of 0.25 is used in all our SIESTA
calculations.

##  VI Acknowledgements

We thank Byungmin Sohn and Yu He for helpful discussions and comments on the
manuscript. This work was supported by grant NSF DMR 2237469 and NSF ACCESS
supercomputing resources via allocation TG-MCA08X007.

## References

  * Kastner et al. (1998) M. Kastner, R. Birgeneau, G. Shirane, and Y. Endoh, Reviews of Modern Physics 70, 897 (1998).
  * Ho et al. (2001) C.-M. Ho, V. Muthukumar, M. Ogata, and P. Anderson, Physical Review Letters 86, 1626 (2001).
  * Warren Jr et al. (1989) W. Warren Jr, R. Walstedt, G. Brennert, R. Cava, R. Tycko, R. Bell, and G. Dabbagh, Physical review letters 62, 1193 (1989).
  * Homes et al. (1993) C. Homes, T. Timusk, R. Liang, D. Bonn, and W. Hardy, Physical review letters 71, 1645 (1993).
  * Loeser et al. (1996) A. Loeser, Z.-X. Shen, D. Dessau, D. Marshall, C. Park, P. Fournier, and A. Kapitulnik, Science 273, 325 (1996).
  * Ding et al. (1996) H. Ding, T. Yokoya, J. C. Campuzano, T. Takahashi, M. Randeria, M. Norman, T. Mochiku, K. Kadowaki, and J. Giapintzakis, Nature 382, 51 (1996).
  * Andersen et al. (1994) O. Andersen, O. Jepsen, A. Liechtenstein, and I. Mazin, Physical Review B 49, 4145 (1994).
  * Liechtenstein et al. (1996) A. Liechtenstein, O. Gunnarsson, O. Andersen, and R. Martin, Physical Review B 54, 12505 (1996).
  * He et al. (2021) Y. He, S.-D. Chen, Z.-X. Li, D. Zhao, D. Song, Y. Yoshida, H. Eisaki, T. Wu, X.-H. Chen, D.-H. Lu, et al., Physical Review X 11, 031068 (2021).
  * Li et al. (2021) Z.-X. Li, S. A. Kivelson, and D.-H. Lee, npj Quantum Materials 6, 36 (2021).
  * Gurvitch and Fiory (1987) M. Gurvitch and A. Fiory, Physical review letters 59, 1337 (1987).
  * Martin et al. (1990) S. Martin, A. T. Fiory, R. Fleming, L. Schneemeyer, and J. V. Waszczak, Physical Review B 41, 846 (1990).
  * Cooper et al. (2009) R. Cooper, Y. Wang, B. Vignolle, O. Lipscombe, S. Hayden, Y. Tanabe, T. Adachi, Y. Koike, M. Nohara, H. Takagi, et al., Science 323, 603 (2009).
  * Greene et al. (2020) R. L. Greene, P. R. Mandal, N. R. Poniatowski, and T. Sarkar, Annual Review of Condensed Matter Physics 11, 213 (2020).
  * Ayres et al. (2021) J. Ayres, M. Berben, M. Čulo, Y.-T. Hsu, E. van Heumen, Y. Huang, J. Zaanen, T. Kondo, T. Takeuchi, J. Cooper, et al., Nature 595, 661 (2021).
  * Dai et al. (2001) P. Dai, H. A. Mook, R. D. Hunt, and F. Doğan, Physical Review B 63, 054525 (2001).
  * Le Tacon et al. (2011) M. Le Tacon, G. Ghiringhelli, J. Chaloupka, M. M. Sala, V. Hinkov, M. Haverkort, M. Minola, M. Bakr, K. Zhou, S. Blanco-Canosa, et al., Nature Physics 7, 725 (2011).
  * Dean et al. (2013) M. Dean, G. Dellea, R. S. Springell, F. Yakhou-Harris, K. Kummer, N. Brookes, X. Liu, Y. Sun, J. Strle, T. Schmitt, et al., Nature materials 12, 1019 (2013).
  * Keimer et al. (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Kivelson et al. (1998) S. A. Kivelson, E. Fradkin, and V. J. Emery, Nature 393, 550 (1998).
  * Lee et al. (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Reviews of modern physics 78, 17 (2006).
  * Chang et al. (2012) J. Chang, E. Blackburn, A. Holmes, N. B. Christensen, J. Larsen, J. Mesot, R. Liang, D. Bonn, W. Hardy, A. Watenphul, et al., Nature Physics 8, 871 (2012).
  * Emery and Kivelson (1995) V. Emery and S. Kivelson, Nature 374, 434 (1995).
  * Corson et al. (1999) J. Corson, R. Mallozzi, J. Orenstein, J. N. Eckstein, and I. Bozovic, Nature 398, 221 (1999).
  * Li et al. (2010) L. Li, Y. Wang, S. Komiya, S. Ono, Y. Ando, G. Gu, and N. Ong, Physical Review B 81, 054510 (2010).
  * Dubroka et al. (2011) A. Dubroka, M. Rössle, K. W. Kim, V. K. Malik, D. Munzar, D. N. Basov, A. A. Schafgans, S. Moon, C. Lin, D. Haug, et al., Physical review letters 106, 047006 (2011).
  * White (1992) S. R. White, Physical review letters 69, 2863 (1992).
  * Schollwöck (2005) U. Schollwöck, Reviews of modern physics 77, 259 (2005).
  * Ceperley and Alder (1986) D. Ceperley and B. Alder, Science 231, 555 (1986).
  * Foulkes et al. (2001) W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Reviews of Modern Physics 73, 33 (2001).
  * Gao et al. (1988) Y. Gao, P. Lee, P. Coppens, M. Subramania, and A. Sleight, Science 241, 954 (1988).
  * Mans et al. (2006) A. Mans, I. Santoso, Y. Huang, W. Siu, S. Tavaddod, V. Arpiainen, M. Lindroos, H. Berger, V. Strocov, M. Shi, et al., Physical Review Letters 96, 1103 (2006).
  * Andersen et al. (2007) B. M. Andersen, P. Hirschfeld, and J. A. Slezak, Physical Review B 76, 020507 (2007).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Physical review 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Physical review 140, A1133 (1965).
  * Hybertsen et al. (1989) M. S. Hybertsen, M. Schlüter, and N. E. Christensen, Physical Review B 39, 9028 (1989).
  * Martin (1996) R. L. Martin, Physical Review B 53, 15501 (1996).
  * Markiewicz et al. (2005) R. Markiewicz, S. Sahrakorpi, M. Lindroos, H. Lin, and A. Bansil, Physical Review B 72, 054519 (2005).
  * Korshunov et al. (2005) M. Korshunov, V. Gavrichkov, S. Ovchinnikov, I. Nekrasov, Z. Pchelkina, and V. Anisimov, Physical Review B 72, 165104 (2005).
  * Delannoy et al. (2009) J.-Y. Delannoy, M. Gingras, P. C. Holdsworth, and A.-M. Tremblay, Physical Review B 79, 235130 (2009).
  * Hanke et al. (2010) W. Hanke, M. L. Kiesel, M. Aichhorn, S. Brehm, and E. Arrigoni, The European Physical Journal Special Topics 188, 15 (2010).
  * Varignon et al. (2019) J. Varignon, M. Bibes, and A. Zunger, Nature communications 10, 1658 (2019).
  * Wang et al. (2021) Z. Wang, O. I. Malyi, X. Zhao, and A. Zunger, Physical Review B 103, 165110 (2021).
  * Zunger (2022) A. Zunger, Nature Computational Science 2, 529 (2022).
  * Maeda et al. (1988) H. Maeda, Y. Tanaka, M. Fukutomi, and T. Asano, Japanese Journal of Applied Physics 27, L209 (1988).
  * Petricek et al. (1990) V. Petricek, Y. Gao, P. Lee, and P. Coppens, Physical Review B 42, 387 (1990).
  * Levin et al. (1994) A. Levin, Y. I. Smolin, and Y. F. Shepelev, Journal of Physics: Condensed Matter 6, 3539 (1994).
  * Slezak et al. (2008) J. Slezak, J. Lee, M. Wang, K. McElroy, K. Fujita, B. Andersen, P. Hirschfeld, H. Eisaki, S. Uchida, and J. Davis, Proceedings of the National Academy of Sciences 105, 3203 (2008).
  * He et al. (2008) Y. He, S. Graser, P. Hirschfeld, and H.-P. Cheng, Physical Review B 77, 220507 (2008).
  * Wang et al. (2005) L.-L. Wang, P. Hirschfeld, and H.-P. Cheng, Physical Review B 72, 224516 (2005).
  * Lin et al. (2006) H. Lin, S. Sahrakorpi, R. Markiewicz, and A. Bansil, Physical review letters 96, 097001 (2006).
  * Foyevtsova et al. (2010) K. Foyevtsova, H. C. Kandpal, H. O. Jeschke, S. Graser, H.-P. Cheng, R. Valentí, and P. Hirschfeld, Physical Review B 82, 054514 (2010).
  * Fan and Zeng (2011) W. Fan and Z. Zeng, Superconductor Science and Technology 24, 105007 (2011).
  * He et al. (2006) Y. He, T. Nunner, P. Hirschfeld, and H.-P. Cheng, Physical review letters 96, 197002 (2006).
  * Song et al. (2019) D. Song, X. Zhang, C. Lian, H. Liu, I. Alexandrou, I. Lazić, E. G. Bosch, D. Zhang, L. Wang, R. Yu, et al., Advanced Functional Materials 29, 1903843 (2019).
  * Zhang et al. (2020) Y. Zhang, C. Lane, J. W. Furness, B. Barbiellini, J. P. Perdew, R. S. Markiewicz, A. Bansil, and J. Sun, Proceedings of the National Academy of Sciences 117, 68 (2020).
  * Liang et al. (1988) J. Liang, S. Xie, G. Che, J. Huang, Y. Zhang, and Z. Zhao, Modern Physics Letters B 2, 483 (1988).
  * Jain et al. (2013) A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, et al., APL Materials 1 (2013), ISSN 2166-532X, 011002.
  * (59) https://next-gen.materialsproject.org/materials/mp-1218930?chemsys=Bi-Sr-Cu-Ca-O, data retrieved from the Materials Project for Sr2CaCu2(BiO4)2 (mp-1218930) from database version v2022.10.28.
  * Yelpo et al. (2021) C. Yelpo, R. Faccio, D. Ariosa, and S. Favre, Journal of Physics: Condensed Matter 33, 185705 (2021).
  * Subramanian et al. (1988) M. Subramanian, C. Torardi, J. Calabrese, J. Gopalakrishnan, K. Morrissey, T. Askew, R. Flippen, U. Chowdhry, and A. Sleight, Science 239, 1015 (1988).
  * Schrieffer and Brooks (2007) J. R. Schrieffer and J. S. Brooks, _Handbook of high-temperature superconductivity: theory and experiment_ (Springer, 2007).
  * Zeljkovic et al. (2012) I. Zeljkovic, E. J. Main, T. L. Williams, M. Boyer, K. Chatterjee, W. Wise, Y. Yin, M. Zech, A. Pivonka, T. Kondo, et al., Nature materials 11, 585 (2012).
  * Wang et al. (2023) Z. Wang, C. Zou, C. Lin, X. Luo, H. Yan, C. Yin, Y. Xu, X. Zhou, Y. Wang, and J. Zhu, Science 381, 227 (2023).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Physical review B 56, 12847 (1997).
  * Pizzi et al. (2020) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, et al., Journal of Physics: Condensed Matter 32, 165902 (2020).
  * (67) Z. Jin and S. Ismail-Beigi, Supplementary Information (????).
  * Massee et al. (2020) F. Massee, Y. Huang, and M. Aprili, Science 367, 68 (2020).
  * Anisimov et al. (2004) V. Anisimov, M. Korotin, A. Mylnikova, A. Kozhevnikov, D. M. Korotin, and J. Lorenzana, Physical Review B 70, 172501 (2004).
  * Pesant and Côté (2011) S. Pesant and M. Côté, Physical Review B 84, 085104 (2011).
  * Moriya (2006) T. Moriya, Proceedings of the Japan Academy, Series B 82, 1 (2006).
  * Dahm et al. (2009) T. Dahm, V. Hinkov, S. Borisenko, A. Kordyuk, V. Zabolotnyy, J. Fink, B. Büchner, D. Scalapino, W. Hanke, and B. Keimer, Nature Physics 5, 217 (2009).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Computational materials science 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Physical review B 54, 11169 (1996b).
  * Damascelli et al. (2003) A. Damascelli, Z. Hussain, and Z.-X. Shen, Reviews of modern physics 75, 473 (2003).
  * Brouet et al. (2008) V. Brouet, W. Yang, X. Zhou, Z. Hussain, R. Moore, R. He, D. Lu, Z. Shen, J. Laverock, S. Dugdale, et al., Physical Review B 77, 235104 (2008).
  * Rotenberg et al. (2005) E. Rotenberg, B. Freelon, H. Koh, A. Bostwick, K. Rossnagel, A. Schmid, and S. Kevan, New Journal of Physics 7, 114 (2005).
  * Kampf and Schrieffer (1990) A. P. Kampf and J. Schrieffer, Physical Review B 42, 7967 (1990).
  * Langer et al. (1995) M. Langer, J. Schmalian, S. Grabowski, and K. Bennemann, Physical review letters 75, 4508 (1995).
  * Chubukov (1995) A. V. Chubukov, Physical Review B 52, R3840 (1995).
  * Koitzsch et al. (2004) A. Koitzsch, S. V. Borisenko, A. Kordyuk, T. K. Kim, M. Knupfer, J. Fink, M. Golden, W. Koops, H. Berger, B. Keimer, et al., Physical Review B 69, 220505 (2004).
  * Nakayama et al. (2006) K. Nakayama, T. Sato, T. Dobashi, K. Terashima, S. Souma, H. Matsui, T. Takahashi, J. Campuzano, K. Kudo, T. Sasaki, et al., Physical review B 74, 054505 (2006).
  * Yu et al. (2020) T. Yu, C. E. Matt, F. Bisti, X. Wang, T. Schmitt, J. Chang, H. Eisaki, D. Feng, and V. N. Strocov, npj Quantum Materials 5, 46 (2020).
  * Stepanov et al. (2018) E. A. Stepanov, L. Peters, I. S. Krivenko, A. I. Lichtenstein, M. I. Katsnelson, and A. N. Rubtsov, npj Quantum Materials 3, 54 (2018).
  * Zhu et al. (2023) M. Zhu, D. Voneshen, S. Raymond, O. Lipscombe, C. Tam, and S. Hayden, Nature Physics 19, 99 (2023).
  * Sobota et al. (2021) J. A. Sobota, Y. He, and Z.-X. Shen, Reviews of Modern Physics 93, 025006 (2021).
  * Ku et al. (2010) W. Ku, T. Berlijn, C.-C. Lee, et al., Physical review letters 104, 216401 (2010).
  * Brouet et al. (2012) V. Brouet, M. F. Jensen, P.-H. Lin, A. Taleb-Ibrahimi, P. Le Fèvre, F. Bertran, C.-H. Lin, W. Ku, A. Forget, and D. Colson, Physical Review B 86, 075123 (2012).
  * Lin et al. (2011) C.-H. Lin, T. Berlijn, L. Wang, C.-C. Lee, W.-G. Yin, and W. Ku, Physical review letters 107, 257001 (2011).
  * Medeiros et al. (2014) P. V. Medeiros, S. Stafström, and J. Björk, Physical Review B 89, 041407 (2014).
  * Tomić et al. (2014) M. Tomić, H. O. Jeschke, and R. Valentí, Physical Review B 90, 195121 (2014).
  * Zhu et al. (2018) Z. Zhu, T.-R. Chang, C.-Y. Huang, H. Pan, X.-A. Nie, X.-Z. Wang, Z.-T. Jin, S.-Y. Xu, S.-M. Huang, D.-D. Guan, et al., Nature communications 9, 4153 (2018).
  * He et al. (2018) Y. He, M. Hashimoto, D. Song, S.-D. Chen, J. He, I. Vishik, B. Moritz, D.-H. Lee, N. Nagaosa, J. Zaanen, et al., Science 362, 62 (2018).
  * Liu et al. (2019) J. Liu, L. Zhao, Q. Gao, P. Ai, L. Zhang, T. Xie, J.-W. Huang, Y. Ding, C. Hu, H.-T. Yan, et al., Chinese Physics B 28, 077403 (2019).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical review b 59, 1758 (1999).
  * Anisimov et al. (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Physical Review B 44, 943 (1991).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and Y. Wang, Physical review B 54, 16533 (1996).
  * Gotlieb et al. (2018) K. Gotlieb, C.-Y. Lin, M. Serbyn, W. Zhang, C. L. Smallwood, C. Jozwiak, H. Eisaki, Z. Hussain, A. Vishwanath, and A. Lanzara, Science 362, 1271 (2018).
  * Mattheiss and Hamann (1988) L. Mattheiss and D. Hamann, Physical Review B 38, 5012 (1988).
  * Hybertsen and Mattheiss (1988) M. S. Hybertsen and L. Mattheiss, Physical review letters 60, 1661 (1988).
  * Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Physical Review B 23, 5048 (1981).
  * Wang et al. (2006) L. Wang, T. Maxisch, and G. Ceder, Physical Review B 73, 195107 (2006).
  * Deng et al. (2019) L. Deng, Y. Zheng, Z. Wu, S. Huyan, H.-C. Wu, Y. Nie, K. Cho, and C.-W. Chu, Proceedings of the National Academy of Sciences 116, 2004 (2019).
  * Nokelainen et al. (2020) J. Nokelainen, C. Lane, R. S. Markiewicz, B. Barbiellini, A. Pulkkinen, B. Singh, J. Sun, K. Pussi, and A. Bansil, Physical Review B 101, 214523 (2020).
  * Sun et al. (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Physical review letters 115, 036402 (2015).
  * Gautam and Carter (2018) G. S. Gautam and E. A. Carter, Physical Review Materials 2, 095401 (2018).
  * Long et al. (2020) O. Y. Long, G. S. Gautam, and E. A. Carter, Physical Review Materials 4, 045401 (2020).
  * Soler et al. (2002) J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, Journal of Physics: Condensed Matter 14, 2745 (2002).