##  1 Introduction

In General Relativity, the simplest black hole solution is that provided by
the Schwarzschild metric [1], where the only parameter of the solution is the
mass $m$. Less then half a century later a more general solution was obtained
[2], which describes a black hole with three parameters: mass, angular
momentum, and charge. The introduction of a cosmological constant admits
further solutions, as do other matter sources.

In 4 dimensions the assumption that the metric has radial symmetry restricts
black hole solutions to have horizons of constant curvature in Einstein
gravity. These can be spherical, planar, or hyperbolic if a negative
cosmological constant $\Lambda<0$ is introduce [3, 4, 5]. If the geometry is
hyperbolic then these black holes can have negative mass [6] and can form from
gravitational collapse [7].

However if $\Lambda\geq 0$ then the situation is notably different. While it
has recently been shown that de Sitter Schwarzschild spacetimes with negative
mass in [8, 9] can exist in Einstein gravity, these solutions consist of
spacetime bubbles. No negative mass black hole vacuum solutions have been
obtained.

We show in this paper that asymptotically de Sitter black hole vacuum
solutions of negative mass exist in Lovelock gravity. This somewhat unexpected
finding arises from a feature recently pointed out for black holes in higher-
curvature gravity: their event horizons need not be of constant curvature [10,
11, 12, 13, 14, 15]. The transverse space of the (potential) black hole
solution need only be an Einstein space, admitting horizon geometries of
increasing complexity as the spacetime dimension increases. For this reason
these have been referred to as Exotic Black Holes (EBHs).

We demonstrate that static asymptotically de Sitter Lovelock EBHs can have
negative mass. We explicitly illustrate this for Gauss-Bonnet gravity, but it
is clear that this will occur in higher-order Lovelock gravity and we expect
in more general higher-curvature theories as well. Lovelock gravity is a
particular subclass of higher curvature gravity theories, which have garnered
much attention since quantum gravity induces such higher-order corrections to
the standard Einstein-Hilbert gravitational action [16]. Furthermore,
renormalization properties are improved for such theories [17]. The Lovelock
class [18] is of particular interest: it is regarded as the most natural
higher-curvature generalization of Einstein gravity whose field equations are
of second order in the metric functions.

Research into EBHs has been of recent interest particularly in the
thermodynamics of asymptotically Anti de Sitter black holes. The geometry of
event horizon was shown to greatly affect the thermodynamic phenomena [19, 20,
21, 22]. While investigation of black holes in AdS space has been of great
interest since the discovery the AdS/Conformal Field Theory correspondence
[23, 24], interest in de Sitter black holes is motivated primarily by
cosmological considerations [25], though some motivation comes from
considerations of holographic duality [26].

We begin with a brief review of Lovelock gravity and EBHs before moving onto
Gauss-Bonnet solutions. We provide general polynomial equations for
determining horizon locations as well as the Kretschmann scalar, which will
function as a diagnostic for finding spacetime singularties. We then discuss
the positive mass case and how horizon geometry can dictate whether or not a
black hole may exist. We then provide an analysis of the negative mass case,
where we obtain a bound on the minimum allowed mass for obtaining a black
hole. We also present a closed form expression that dictates the bounds on the
horizon geometry for which a negative mass black hole is admitted as a vacuum
solution. We have shown that it is possible to obtain massless de Sitter black
holes with the proper set of parameters, and discuss this in the penultimate
section before summarizing our work.

##  2 Lovelock Gravity & Exotic Black Holes

For a Lovelock theory of gravity the Lagrangian density is given by

$\mathcal{L}=\frac{1}{16\pi
G_{N}}\sum_{k=0}^{K}\hat{\alpha}_{k}\mathcal{L}^{(k)}$ (2.1)

where $\hat{\alpha}_{k}$ are the Lovelock coupling constants and
$\mathcal{L}^{(k)}$ are the Euler densities with dimension $2k$, which are

$\mathcal{L}^{(k)}=\frac{1}{2^{k}}\delta_{c_{1}d_{1}\ldots
c_{k}d_{k}}^{a_{1}b_{1}\ldots a_{k}b_{k}}R_{a_{1}b_{1}}^{c_{1}d_{1}}\ldots
R_{a_{k}b_{k}}^{c_{k}d_{k}}$ (2.2)

with $\delta$ representing the generalized fully anti-symmetric Kronecker
delta. The first few terms are

$\mathcal{L}^{(0)}=1\qquad\mathcal{L}^{(1)}=R\qquad\mathcal{L}^{(2)}=R^{2}-4R_{ab}R^{ab}+R_{abcd}R^{abcd}$

with $R$ the Ricci Scalar and $\mathcal{L}^{(2)}$ the Gauss-Bonnet term. From
the Lagrangian density (2.1) we may construct an action

$S=\int d^{d}x\sqrt{-g}\left(\frac{1}{16\pi
G_{N}}\sum_{k=0}^{K}\alpha_{k}\mathcal{L}^{(k)}+\mathcal{L}_{m}\right)$ (2.3)

including a matter term $\mathcal{L}_{m}$. Variation with respect to the
metric $g^{ab}$ yields the field equations

$\sum_{k=0}^{K}\hat{\alpha}_{(k)}\mathcal{G}_{ab}^{(k)}=8\pi G_{N}T_{ab}.$
(2.4)

where $T_{ab}$ is the stress energy tensor of the matter field and
$\mathcal{G}_{ab}^{(k)}$ are the Lovelock tensors

$\mathcal{G}_{ab}^{(k)}=-\frac{1}{2^{(k+1)}}g_{ab}\delta_{e_{1}f_{1}\ldots
e_{k}f_{k}}^{c_{1}d_{1}\ldots c_{k}d_{k}}R_{c_{1}d_{1}}^{e_{1}f_{1}}\ldots
R_{c_{k}d_{k}}^{e_{k}f_{k}}$ (2.5)

We will consider only vacuum solutions with $\mathcal{L}_{m}=0$ henceforth, so
the field equations will be (2.4) with the stress energy being set to zero.
The parameter $K\leq\frac{d-1}{2}$, where $d$ is the dimension of spacetime,
and sets the maximal degree of non-linearity in the curvature. Values of $K$
larger than this make no contribution to the field equations and are
topological invariants in the action.

The metric ansatz we will make use of is

$\displaystyle ds^{2}$
$\displaystyle=\textsf{g}_{ij}dy^{i}dy^{j}+\gamma_{\alpha\beta}dx^{\alpha}dx^{\beta}$
(2.6) $\displaystyle=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Sigma_{d-2}^{2}$

where $y^{i}=(t,r)$ and $\textsf{g}_{ij}=\mathrm{diag}(-f(r),\frac{1}{f(r)})$.
The line element $d\Sigma^{2}_{d-2}$ of the base manifold (or transverse
space) has metric $\gamma_{\alpha\beta}$ and coordinates $x^{\alpha}$. We
shall assume it to be compact, with volume $\Sigma_{d-2}$, the analogue of the
volume of a unit sphere. If this space is not compact then quantities such as
$M/\Sigma_{d-2}$ will be regarded as finite.

The only other condition on the transverse space is that it is an Einstein
space, $R_{\alpha\beta}\propto\gamma_{\alpha\beta}$. The most commonly studied
special case is when the base space has constant curvature, with curvature
parameter $\kappa=(-1,0,1)$, which corresponds to negative, flat and positive
curvature respectively. Solutions in which this latter condition is relaxed,
so that the transverse space is only an Einstein space [11], shall be referred
to as Exotic Black Holes [10, 27].

Introducing this metric into the field equations yields a resulting polynomial
equation for $f(r)$ [10]

$\sum_{n=0}^{K}\frac{b_{n}}{r^{2n}}\left(\sum_{k=n}^{K}\alpha_{k}\binom{k}{n}\left(\frac{-f(r)}{r^{2}}\right)^{k-n}\right)=\frac{16\pi
G_{N}M}{(d-2)\Sigma_{d-2}r^{d-1}}$ (2.7)

where $M$ represents the mass of the black hole and $\Sigma_{d-2}$ the volume
of the base space. The quantities $\alpha_{k}$ are rescaled Lovelock
couplings, defined as

$\alpha_{0}=\frac{\hat{\alpha}_{0}}{(d-1)(d-2)},\quad\alpha_{1}=\hat{\alpha}_{1},\quad\alpha_{k}=\hat{\alpha}_{k}\prod_{n=1}^{2k}(d-n)\quad\text{
for }k\geq 2.$ (2.8)

The cosmological constant is defined through the zeroth Lovelock coupling
constant $\Lambda=-\frac{\hat{\alpha}_{0}}{2}$.

Imposing only the condition that the transverse space is an Einstein space
admits the following possibilities

$\displaystyle\hat{\mathcal{G}}_{\beta}^{(n)\alpha}$
$\displaystyle=-\frac{(d-3)!b_{n}}{2(d-2n-3)!}\delta_{\beta}^{\alpha}$ (2.9)
$\displaystyle\hat{\mathcal{L}}^{(n)}$
$\displaystyle=\frac{(d-2)!b_{n}}{(d-2n-2)!}$

for its Lovelock tensor and associated intrinsic Euler density respectively.
The constants $b_{n}$ we refer to as topological terms, and can take any value
in $\mathbb{R}$. They define the topology of the base manifold and the
geometry of the event horizon from the condition $f(r_{+})=0$ that defines the
horizon location $r_{+}$. Without loss of generality we set $b_{0}=1$; to
recover Einstein gravity in the proper limit of vanishing $\alpha_{k\geq 2}$,
we will set $\alpha_{1}=1$. The mass $M$ is defined through the Hamiltonian
formalism, and is the conserved charge of the timelike Killing vector of the
background spacetime. When the dimension is $d=2K+1$ the transverse space is
restricted to be of constant curvature [28, 29].

The Kretschmann scalar for a metric of the form (2.6) is

$R^{abcd}R_{abcd}=\left(\frac{d^{2}f(r)}{dr^{2}}\right)^{2}+2\frac{(d-2)}{r^{2}}\left(\frac{df(r)}{dr}\right)^{2}+2\frac{(d-2)(d-3)f(r)^{2}}{r^{4}}-4\frac{R[\gamma]f(r)}{r^{4}}+\frac{\mathcal{K}[\gamma]}{r^{4}}$
(2.10)

where $R[\gamma]$ and $\mathcal{K}[\gamma]$ are the Ricci and Kretschmann
scalars of the tranverse space respectively.

##  3 Gauss-Bonnet Solutions

Setting $K=2$ in (2.7) yields

$\frac{\alpha_{2}f^{2}}{r^{4}}+\left(-\frac{1}{r^{2}}-\frac{2b_{1}\alpha_{2}}{r^{4}}\right)f+\alpha_{0}+\frac{b_{1}}{r^{2}}+\frac{b_{2}\alpha_{2}}{r^{4}}=\frac{16\pi
M}{(d-2)\Sigma_{d-2}r^{d-1}}$ (3.1)

with the solutions

$f=f_{\pm}(\mathrm{m})\equiv\frac{r^{2}+2b_{1}\alpha_{2}\pm\sqrt{\left(b_{1}^{2}-b_{2}\right)4\alpha_{2}^{2}+r^{4}\left(1-4\alpha_{2}\alpha_{0}\right)+\frac{8\mathrm{m}\alpha_{2}}{r^{d-5}}}}{2\alpha_{2}}$
(3.2)

where we have written

$\mathrm{m}\equiv\frac{8\pi M}{(d-2)\Sigma_{d-2}}.$ (3.3)

The two solutions are distinguished by their behaviour in the
$\alpha_{2}\rightarrow 0$ limit. The $f_{-}$ branch is referred to as the
Einstein branch as it has a smooth $\alpha_{2}\rightarrow 0$ limit, giving

$\lim_{\alpha_{2}\rightarrow
0}f_{-}(\mathrm{m})=\alpha_{0}r^{2}+b_{1}-\frac{2\mathrm{m}}{r^{d-3}}$ (3.4)

whereas the $f_{+}$ solution does not have a smooth limit as
$\alpha_{2}\rightarrow 0$; it is referred to as the Gauss-Bonnet branch.

For the remainder of this paper we will only consider the Einstein branch for
analysis; we will also choose $b_{1}=1$, so that as $\alpha_{2}\rightarrow 0$
we recover the standard Schwarzschild (Anti) de Sitter solution. We will only
examine solutions with $d\geq 6$, since the $d=5$ case only admits constant
curvature solutions in the transverse space. The cosmological constant is

$\displaystyle\Lambda$
$\displaystyle=-\frac{\hat{\alpha}_{0}}{2}=-\frac{\alpha_{0}(d-1)(d-2)}{2}$
(3.5)

and we shall set $\alpha_{0}<0$ since we are interested in asymptotically de
Sitter solutions. This allows us to write

$\displaystyle f_{-}$
$\displaystyle=\frac{r^{2}+2\alpha_{2}-\sqrt{\left(1-b_{2}\right)4\alpha_{2}^{2}+r^{4}\left(1+4\alpha_{2}\frac{2\Lambda}{(d-1)(d-2)}\right)+\frac{8\mathrm{m}\alpha_{2}}{r^{d-5}}}}{2\alpha_{2}}$
(3.6)

where we note $\Lambda>0$.

It is convenient at this point to introduce a set of dimensionless variables

$r=x\sqrt{\alpha_{2}},\quad\Lambda=\frac{(d-1)(d-2)z}{2\alpha_{2}},\quad\mathrm{m}=m\alpha_{2}^{\frac{d-3}{2}},\quad
z>0$ (3.8)

so that

$f=1+\frac{x^{2}}{2}-\frac{\sqrt{(1+4z)x^{4}+4(1-b_{2})+\frac{8m}{x^{d-5}}}}{2}.$
(3.9)

The horizon(s) are located at the roots of $f$, which can be found from
solving the polynomial:

$x^{d-1}z-x^{d-3}-b_{2}x^{d-5}+2m=0\quad\mathrm{with}\quad d\geq 6,\quad z>0$
(3.10)

For a black hole solution we must have $f(x)>0\;\mathrm{for}\;x>x_{+}$ for
some range of $x$, where $x_{+}$ locates the outermost event horizon of the
black hole.

Superfically it appears that there will be 2 or 0 positive roots to this
polynomial from Descartes’ rule of signs. However if we allow both $b_{2}$ and
$m$ to be either positive or negative the number of possible positive roots
are then the following:

$\displaystyle\mathrm{If}\ m>0$ $\displaystyle\mathrm{then}\
2\;\textrm{roots}\rightarrow(x_{c},x_{+})\ \mathrm{or}\ 0\ \textrm{roots for
any}\;b_{2}$ $\displaystyle\mathrm{If}\ m<0\ \textrm{and}\ b_{2}>0,$
$\displaystyle\mathrm{then}\ 1\ \mathrm{root}\rightarrow x_{c}\hfill$
$\displaystyle\mathrm{If}\ m<0\ \textrm{and}\ b_{2}<0,$
$\displaystyle\mathrm{then}\ 3\ \mathrm{roots}\rightarrow(x_{-},x_{+},x_{c})\
\mathrm{or}\ 1\ \mathrm{root}\rightarrow x_{c}.$ (3.11)

where $x_{c}$ is the cosmological horizon, $x_{-}$ is the inner horizon, and
$x_{+}$ is the outer horizon of the black hole.

###  3.1 Positive Mass

Consider first positive mass. We plot the behaviour of $f$ in Figure 1 for
sample values $b_{2}=-2.5$, $z=0.05$, and $d=6$, varying $m$. For small
nonzero $m$ we have an asymptotically de Sitter black hole (left diagram in
Figure 1) where as for sufficiently large $m$ there is no static black hole
solution and we have an expanding spacetime with $r$ and $t$ interchanging
roles (middle diagram). As $m$ increases from small to large values, there is
one value of $m$ where $f$ has a double root and the horizons merge,
illustrated in the right diagram in Figure 1. This is a Gauss-Bonnet
generalization of the Nariai solution. This structure is the same in all
dimensions $d\geq 6$. Only quantitative differences are seen as we increase
the dimension: the spread between $x_{+}$ and $x_{c}$ increases and the
maximum value of $f$ slightly increases. Consequently we will illustrate our
results in 6 dimensions unless otherwise stated.

Figure 1: Positive Mass Black Hole $b_{2}=-2.5,z=0.5,d=6$. Left: Black hole
with mass $m=1$ showing two distinct real roots. Center: Mass $m=5$, showing
no roots, and a naked singularity. Right: Mass $m=4.0638$, showing a double
root.

We may also examine what happens when we hold all the parameters fixed except
for $b_{2}$, which is seen in Figure 2. For a fixed mass $m$ and cosmological
parameter $z$ increasing the value of $b_{2}$ will shift $f$ upwards. In this
manner the geometry of the horizon plays an important role regarding the
existence of a black hole solution, namely that a minimum allowed value of
$b_{2}$ is required. This minimum value occurs when $f$ has a double root,
corresponding to a Nariai type of black hole solution. In de Sitter space this
is referred to as the maximum allowed mass a black hole can have, other
parameters being fixed. To determine the value of $x_{+}=x_{c}$ for the
Nariai-type solution, the equations

$f=0,\quad\frac{\partial f}{\partial x}=0$ (3.12)

must be simultaneously solved. In 6 dimensions there is no closed form
solution since the latter equation is quintic polynomal is $x$. For the
parameters in Figure 2 we obtain numerically $b_{{2}_{min}}\approx-1.62$. For
$b_{2}$ less than this the spacetime will have a naked singularity at the
origin.

We can also make another observation, shown in Figure 3. That is, for a
certain set of parameters $f$ can become discontinuous for some range of $x$.
This discontinuity is due to the argument inside the square root becoming
negative. Whenever this occures the Kretschmann scalar (2.10) diverges,
yielding a spacetime singularity. Solutions to the following polynomial

$\left(4-4b_{2}\right)x^{d-5}+(4z+1)x^{d-1}+8m=0$ (3.13)

yield the singularity location(s). From Decartses’ rule of signs this will
only have a root when $b_{2}>1$. This is a necessary but not sufficient
condition; as illustrated in the right image of Figure 3, increasing the value
of $z$ allows the two discontinuous branches to join, yielding a smooth
continuous metric function $f$. If we hold all parameters fixed, increasing
$z$ shifts $f$ down vertically and reduces the spread between the black hole
and cosmological horizons, allowing the possibility of a smooth metric
function between them.

Figure 2: Positive Mass $m=1,z=0.1,d=6$. Left: Naked singularity with
$b_{2}=-2$ showing no real roots. Right: de Sitter black hole with $b_{2}=2$
show a cosmological horizon and event horizon $x_{c}>x_{+}$

Figure 3: Positve mass  $m=1$,$d=6$,$b_{2}=3.5$. Left:$z=0.1$ Two
discontinuous branches with two spacetime singularities. Right: $z=0.3$ de
Sitter black hole.

###  3.2 Negative Mass

Negative mass solutions have more interesting behaviour, since if $b_{2}<0$
there can now be three roots to (3.10). We illustrate this in Figure 4. The
left image depicts three horizons, while the center and right image show the
two possible cases where horizons can merge. In the center image $x_{c}=x_{+}$
(a Nariai-type solution) and in the right $x_{-}=x_{+}$ (an extremal black
hole with a cosmological horizon). These latter two extremes constrain the
range of allowed values of $b_{2}<0$ that yield black hole solutions of
negative mass. For the case illustrated in Figure 4, we approximately have
$-1.3\geq b_{2}>-3$.

We also find that as $x\to 0$, $f$ does not have a smooth limit, and instead
becomes complex at some finite positive value of $x$. For $m<0$ and $b_{2}<0$
there is only one sign change in the coefficients of the terms in (3.13) and
hence this equation has only one positive root, at which point the Kretschmann
scalar (2.10) diverges. For the admissible range of $b_{2}$, this singularity
is always behind an horizon as can be seen in Figure 4. For the more generic
3-horizon case shown in the left diagram in Figure 4, the solution to (3.13)
is $x=0.1999637662$.

Figure 4: Negative Mass: $m=-0.3$, $z=0.09$, $d=6$ Left: $b_{2}=-2$, showing
three distinct horizons. Center: $b_{2}=-3$, showing one double root depicting
a Nariai like solution Right: $b_{2}=-1.3$, showing an extremal black hole
with a distinct cosmological horizon.

####  3.2.1 Allowed values of $b_{2}$ for valid negative mass black hole
solutions

To determine the the range of $b_{2}$ that allows for a valid black hole
solution of negative mass in any dimension, we must solve the following set of
equations

$\displaystyle f=0,\quad\frac{\partial f}{\partial x}=0\ \mathrm{for}\
b_{2},x\quad\mathrm{with}\quad\\{x,b_{2}|\in\mathbb{R}>0\\}.$ (3.14)

The first equation is given by (3.10) and yields

$b_{2}=x_{+}^{4}z-x_{+}^{2}+2x_{+}^{5-d}m$ (3.15)

where $x_{+}$ solves

$2x_{+}^{d-1}z-x_{+}^{d-3}-m(d-5)=0.$ (3.16)

from the latter equation. There is no general solution to this set of
equations except in certain dimensions.

For $d=6$ equation (3.16) is a quintic, and there is no analytic solution for
arbitrary $m,z$. For the parameters used in Figure 4, we obtain numerically

$\left\\{x_{{1+}}=0.6897087623,\;b_{{2}_{max}}=-1.325264590\right\\},\left\\{x_{{2+}}=2.328863088,\;b_{{2}_{min}}=-3.033847194\right\\}.$
(3.17)

refining the values given above.

In 7 dimensions however, we may find a closed form expression. The equations
become

$\displaystyle b_{2}=\frac{x_{+}^{6}z-x_{+}^{4}+2m}{x_{+}^{2}}$ (3.18)
$\displaystyle 2x_{+}^{6}z-x_{+}^{4}-2m=0.$ (3.19)

which admits analytic solutions for arbitrary $m,z$. The latter equation is a
cubic in $x_{+}^{2}$; writing $m=-|\textbf{m}|$ we get

$2y^{3}z-y^{2}+2|\textbf{m}|=0\quad\text{with}\quad y=x_{+}^{2},\ z>0$ (3.20)

Now through Decartes’ rule of signs, the possible number of distinct positive
roots is either 2 or zero, while the number of negative roots will always be
one; this latter case is inadmissible since we must have $y>0$.

To have two distinct positive real roots for $y$ (and thus the same number of
distinct positive values for $x_{+}$), we require the discriminant

$\Delta_{3}=8|\textbf{m}|\left(1-54z^{2}|\textbf{m}|\right)>0$ (3.21)

for (3.20), implying

$-\frac{1}{54z^{2}}<m<0.$ (3.22)

This condition must be satisfied to have three distinct roots. Writing (3.20)
as a depressed cubic, the solutions (for $x$) are easily written as

$x_{k}=\sqrt{\frac{1}{3z}\cos\left(\frac{1}{3}\arccos\left(1-108\mathbf{|m|}z^{2}\right)-\frac{2\pi
k}{3}\right)-\frac{1}{6z}}\quad k=0,1,2$ (3.23)

For all values of $m$ satisfying (3.21), the only $k=0,1$ yield real
solutions, where $x_{0}>x_{1}$. The value $x_{0}$ corresponds to the location
of the location of the Nariai-like black hole solution whereas $x_{1}$ locates
the extremal black hole of negative mass. From (3.18), the respective minimal
and maximal values of $b_{2}$ are

$\begin{gathered}b_{{2}_{min}}=\frac{8\cos(A)^{3}-432\mathbf{|m|}z^{2}-36\cos(A)^{2}+30\cos(A)-7}{72z^{2}\cos(A)-36z}\\\
A=\frac{1}{3}\arccos\left(1-108\mathbf{|m|}z^{2}\right)\\\ \\\
b_{{2}_{max}}=\frac{8\cos(B)^{3}-432\mathbf{|m|}z^{2}-36\cos(B)^{2}+30\cos(B)-7}{72z^{2}\cos(B)-36z}\\\
B=\frac{1}{3}\arccos\left(1-108\mathbf{|m|}z^{2}\right)-\frac{2\pi}{3}\end{gathered}$
(3.24)

and provided $b_{{2}_{max}}>b_{2}>b_{{2}_{min}}$, there will be a valid
negative mass black hole with two horizons and a de Sitter cosmological
horizon.

We can make a further interesting observation about (3.21). If we saturate the
inequality (3.22) from above, setting set $m=0$, a range of interesting
phenomena can occur. It is possible to have a black hole, a Nariai like black
hole, pure de Sitter space or a naked singularity. These scenarios are
analyzed in the next section.

##  4 Massless Black Holes

The saturation of (3.22) from above allows for the possibility of zero mass
exotic black holes. A intriguing feature of these possible black holes is that
they are independent of the dimension of the spacetime. Taking (3.9) and
setting the mass to zero gives us

$f=1+\frac{x^{2}}{2}-\frac{\sqrt{(1+4z)x^{4}+4-4b_{2}}}{2}.$ (4.1)

The root equation (3.16) is now

$x^{4}z-x^{2}-b_{2}=0$ (4.2)

which is a biquadratic equation in $x^{2}$ so it can be written as

$y^{2}z-y-b_{2}=0.$ (4.3)

We may again analyze the roots of this polynomial by Decartses’ rule of signs.
Assuming $z>0$, then

$\displaystyle\text{for}\ b_{2}<0\qquad 2\ \text{or}\ 0\ \text{roots}$ (4.4)
$\displaystyle\text{for}\ b_{2}\geq 0\qquad 1\ \text{root}$

As a quadratic the discriminant is immediately written down as:
$\Delta_{2}\equiv 4zb_{2}+1$. For a valid black hole solution, we must have
$\Delta_{2}\geq 0$. The saturation of this inequality suggests the Nariai
solution. The unsaturated inequality gives us, similar to the negative mass
case, a condition on $b_{2}$ as a function of $z$ that must be satisfied for a
black hole to exist: $-\frac{1}{4z}<b_{2}<0$. If $b_{2}$ is not within these
bounds there will be a naked singularity located at the origin.

For $b_{2}>0$ the space is pure de Sitter, provided (3.13) has no solutions so
as to avoid singularities. These equations are easily solved for any $d$,
yielding

$x=\frac{\sqrt{2}\left(b_{2}-1\right)^{1/4}}{(1+4z)^{1/4}}.$ (4.5)

and so we conclude that we must have $0\leq b_{2}<1$ to obtain a pure de
Sitter space without any naked singularities. The naked singularity that would
occur for $b_{2}\geq 1$ is at a finite value of $x$ given by (4.5). The
different scenarios are outlined in Figure 5. The 3 leftmost images display
negative values of $b_{2}$ respectively showing a naked singularity, Nariai
black hole and de Sitter black hole. The two right most images are for zero
and positive values of $b_{2}$ respectively, with the former being pure de
Sitter space and the latter a naked singularity outside of the origin. For
increasing values of $b_{2}>0$ the termination point of $f$ at $x=0$ moves up
the vertical axes until it reaches a maximum point at $f(x,b_{2},z)\rightarrow
f(0,1,0.1)=1$.

Figure 5: Massless Soltuions $z=0.1$ From left to right the $b_{2}$ values
are:-3, -2.5, -1, 0, 2. It depicts a naked singularity at the origin, a Nariai
like black hole, a de Sitter black hole, pure de Sitter space, and then a non-
trivial naked singularity outside of the origin, respectively.

##  5 Conclusion

We have shown that metrics of the form (2.6), with $\gamma_{\alpha\beta}$
being the metric of an arbitrary Einstein space, can possess both negative
mass and massless asymptotically de Sitter black hole solutions in Gauss-
Bonnet gravity. These $M\leq 0$ black holes are all exotic black holes, and
their existence depends on the relative values of the parameters $b_{2}$, $z$
and $m$. For the massless case, for there is a direct relationship between
$b_{2}$ and $z$ with $b_{2}$ being bounded from below by $z$. For $M<0$ black
holes there are more restrictions on the parameters. The range of $m$ is
bounded by $z$, as given by (3.22). There is also a minimum and maximum
allowed value of $b_{2}$ admitting such solutions.

We close by noting that the existence of such horizon geometries is still an
open question. This is due to the fact that field equations for the transverse
space are an over constrained PDE system, where certain topological parameters
$b_{n}$ may not provide a solution. However we are not aware of any theorem
forbidding the existence of non-trivial solutions to (2.9).

With the discovery of negative mass black holes in de Sitter space interesting
questions arise. These include the behaviour of geodesics, formation from
gravitational collapse, and possible pair production in the early universe of
such objects. Their thermodynamic behaviour is also likely to yield
interesting surprises in comparison to their Anti de Sitter counterparts [19,
20].

## Acknowledgements

This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada.

## References

  * [1] K. Schwarzschild.  Uber das Gravitationsfeld eines Massenpunktes nach der Einsteinschen Theorie.  Deutsch. Akad. Wiss. Berlin, Kl. Math. Phys. Tech., 1916:189–196, 1916.
  * [2] E. T. Newman, E. Couch, K. Chinnapared, A. Exton, A. Prakash, and R. Torrence.  Metric of a rotating, charged mass.  Journal of Mathematical Physics, 6(6):918–919, 1965.
  * [3] Stefan Aminneborg, Ingemar Bengtsson, Soren Holst, and Peter Peldan.  Making anti-de Sitter black holes.  Class. Quant. Grav., 13:2707–2714, 1996.
  * [4] Robert B. Mann.  Pair production of topological anti-de Sitter black holes.  Class. Quant. Grav., 14:L109–L114, 1997.
  * [5] Robert B. Mann.  Topological black holes: Outside looking in.  Annals Israel Phys. Soc., 13:311, 1997.
  * [6] Robert Mann.  Black holes of negative mass.  Classical and Quantum Gravity, 14(10):2927–2930, oct 1997.
  * [7] W. L. Smith and Robert B. Mann.  Formation of topological black holes from gravitational collapse.  Phys. Rev. D, 56:4942–4947, 1997.
  * [8] Matthew C Johnson, MB Paranjape, Antoine Savard, and Natalia Tapia-Arellano.  Stable, thin wall, negative mass bubbles in de sitter space-time.  General Relativity and Gravitation, 52(8):1–18, 2020.
  * [9] Saoussen Mbarek and MB Paranjape.  Negative mass bubbles in de sitter spacetime.  Physical Review D, 90(10):101502, 2014.
  * [10] Sourya Ray.  Birkhoff’s theorem in lovelock gravity for general base manifolds.  Classical and Quantum Gravity, 32(19):195022, Sep 2015.
  * [11] Gustavo Dotti and Reinaldo J. Gleiser.  Obstructions on the horizon geometry from string theory corrections to einstein gravity.  Physics Letters B, 627(1-4):174–179, Oct 2005.
  * [12] Gustavo Dotti, Julio Oliva, and Ricardo Troncoso.  Exact solutions for the einstein-gauss-bonnet theory in five dimensions: Black holes, wormholes, and spacetime horns.  Physical Review D, 76(6), Sep 2007.
  * [13] GUSTAVO DOTTI, JULIO OLIVA, and RICARDO TRONCOSO.  Vacuum solutions with nontrivial boundaries for the einstein-gauss-bonnet theory.  International Journal of Modern Physics A, 24(08n09):1690–1694, Apr 2009.
  * [14] Gustavo Dotti, Julio Oliva, and Ricardo Troncoso.  Static solutions with nontrivial boundaries for the einstein-gauss-bonnet theory in vacuum.  Physical Review D, 82(2), Jul 2010.
  * [15] Julio Oliva.  All the solutions of the form m2(warped) x $\Sigma$(d-2) for lovelock gravity in vacuum in the chern-simons case.  Journal of Mathematical Physics, 54(4):042501, Apr 2013.
  * [16] N. D. Birrell and P. C. W. Davies.  Quantum Fields in Curved Space.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 1982.
  * [17] K. S. Stelle.  Renormalization of higher-derivative quantum gravity.  Phys. Rev. D, 16:953–969, Aug 1977.
  * [18] D. Lovelock.  The Einstein tensor and its generalizations.  J. Math. Phys., 12:498–501, 1971.
  * [19] Brayden R. Hull and Robert B. Mann.  Thermodynamics of exotic black holes in Lovelock gravity.  Phys. Rev. D, 104(8):084032, 2021.
  * [20] Fil. Hull, Brayden R. Simovi.  Exotic Black Hole Thermodynamics in Third-Order Lovelock Gravity.  In Writing, 2022.
  * [21] N. Farhangkhah and M. H. Dehghani.  Lovelock black holes with nonmaximally symmetric horizons.  Phys. Rev. D, 90(4):044014, 2014.
  * [22] N. Farhangkhah and Z. Dayyani.  Extended phase space thermodynamics for third-order Lovelock black holes with nonmaximally symmetric horizons.  Phys. Rev. D, 104(2):024068, 2021.
  * [23] Edward Witten.  Anti-de Sitter space, thermal phase transition, and confinement in gauge theories.  Adv. Theor. Math. Phys., 2:505–532, 1998.
  * [24] Juan M. Maldacena.  The Large N Limit of Superconformal Field Theories and Supergravity.  Advances in Theoretical and Mathematical Physics, 2:231, January 1998.
  * [25] Planck Collaboration 2018F.  Planck 2018 results. VI. Cosmological parameters.  Astronomy & Astrophysics, 641, 2020.
  * [26] Andrew Strominger.  The ds/cft correspondence.  Journal of High Energy Physics, 2001(10):034, 2001.
  * [27] S. Ray.  Exotic black holes in lovelock gravity.  2015\.
  * [28] Gustavo Dotti, Julio Oliva, and Ricardo Troncoso.  Static solutions with nontrivial boundaries for the Einstein-Gauss-Bonnet theory in vacuum.  Phys. Rev. D, 82:024002, 2010.
  * [29] Gustavo Dotti, Julio Oliva, and Ricardo Troncoso.  Exact solutions for the Einstein-Gauss-Bonnet theory in five dimensions: Black holes, wormholes and spacetime horns.  Phys. Rev. D, 76:064038, 2007.


# A Coupled Random Projection Approach to Large-Scale Canonical Polyadic
Decomposition

Lu-Ming Wang, Ya-Nan Wang, Xiao-Feng Gong, , Qiu-Hua Lin, , and Fei Xiang This
work was supported in part by National Natural Science Foundation of China
under grants 62071082, 61671106, 61871067, and in part by Fundamental Research
Funds for Central Universities, China, under grant DUT20LAB120, DUT20ZD220.
(Corresponding author: Xiao-Feng Gong.)Lu-Ming Wang, Ya-Nan Wang, Xiao-Feng
Gong, and Qiu-Hua Lin were with the School of Information and Communication
Engineering, Dalian University of Technology, Dalian 116024, China
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>qhlin@dlut.edu.cn).Fei Xiang, was with Xiaomi Inc. Artificial Intelligence
Department, AI Lab, Beijing 100085, Chian (e-mail: xiangfei@xiaomi.com).

###### Abstract

We propose a novel algorithm for the computation of canonical polyadic
decomposition (CPD) of large-scale tensors. The proposed algorithm generalizes
the random projection (RAP) technique, which is often used to compute large-
scale decompositions, from one single projection to multiple but coupled
random projections (CoRAP). The proposed CoRAP technique yields a set of
tensors that together admits a coupled CPD (C-CPD) and a C-CPD algorithm is
then used to jointly decompose these tensors. The results of C-CPD are finally
fused to obtain factor matrices of the original large-scale data tensor. As
more data samples are jointly exploited via C-CPD, the proposed CoRAP based
CPD is more accurate than RAP based CPD. Experiments are provided to
illustrate the performance of the proposed approach.

###### Index Terms:

Large-scale tensor, coupled canonical polyadic decomposition, coupled random
projection.

##  I Introduction

With the size and dimension of datasets growing much faster than ever,
decomposition of large-scale datasets has become an important issue in signal
processing and machine learning. As such, large-scale tensor decomposition has
attracted much attention recently, and a number of algorithms[1, 2, 3, 4, 5,
6, 7] have been proposed so far. In the above-mentioned algorithms, various
randomization methods play an important role, and on top of that, random
projection (RAP)[4, 3] has been widely adopted in the decomposition of large-
scale tensors into various models, including canonical polyadic decomposition
(CPD)[3, 4], Tucker[7], tensor train (TT)[6] and tensor ring (TR)[5].

Besides the above works on large-scale tensor decompositions, coupled
decompositions of multiple datasets, including coupled matrix-tensor
decomposition[8], coupled and double coupled CPD[9, 10, 11, 12, 13, 14] and
coupled block term decomposition[15, 16], have also received much attention in
the recent decade. These coupled decomposition techniques are important tools
in multi-set data fusion and are shown to have better performance than their
uncoupled counterpart with regards to both identifiability and accuracy.

In this paper, we combine the idea of RAP and that of coupled tensor
decomposition to develop a novel algorithm for the CPD of a large-scale
tensor. We will propose a new coupled RAP (CoRAP) approach that can generate
multiple coupled sets of projection matrices, and applying these sets of
projection matrices to the data tensor yields multiple compressed tensors. Due
to the CPD structure of the original tensor and the coupling among different
sets of projection matrices, the set of compressed tensors together admits a
coupled CPD (C-CPD). Therefore, the compressed tensors can be jointly
decomposed via a C-CPD algorithm and the results can then be fused to obtain
the factor matrices of the original large-scale data tensor. Noting that the
proposed CoRAP based CPD algorithm exploits more structure, i.e., multiple
compressed tensors as well as their coupling, than its RAP based counterpart,
it is expected to have better performance with regards to accuracy. For
convenience, we limit ourselves to third-order real-valued tensors in this
paper, although the presented results can be analogously extended to complex-
valued tensors with order higher than three.

Notations: vectors, matrices and tensors are denoted by lowercase boldface,
uppercase boldface and uppercase calligraphic letters, respectively. The $r$th
column vector and the $(i,j)$th entry of $\mathbf{A}$ are denoted by
$\mathbf{a}_{r}$ and $a_{i,j}$ , respectively. Symbols ‘$\otimes$’, ‘$\odot$’
and ‘$\times_{n}$’ denote Kronecker product, Khatri-Rao product, mode-$n$
product, and outer product, respectively, defined as:

$\displaystyle\mathbf{A}\otimes\mathbf{B}\triangleq\begin{bmatrix}a_{11}\mathbf{B}&a_{12}\mathbf{B}&\cdots\\\
a_{21}\mathbf{B}&a_{22}\mathbf{B}&\cdots\\\
\vdots&\vdots&\ddots\end{bmatrix},$
$\displaystyle\mathbf{A}\odot\mathbf{B}\triangleq\left[\mathbf{a}_{1}\otimes\mathbf{b}_{1},\mathbf{a}_{2}\otimes\mathbf{b}_{2},\cdots\right],$
$\displaystyle\left(\boldsymbol{\mathcal{T}}\times_{n}\mathbf{G}\right)_{i_{1}...i_{n-1},j,i_{n+1}...i_{N}}\triangleq\sum_{n}t_{i_{1}...i_{N}}g_{j,i_{n}},$
$\displaystyle\left(\mathbf{a}\circ\mathbf{b}\circ\mathbf{c}\right)_{i,j,k}\triangleq
a_{i}b_{j}c_{k}.$

For the mode-$n$ product we assume that the $n$th dimension of
$\boldsymbol{\mathcal{T}}$ is equal to the number of columns of $\mathbf{G}$.
Transpose, Moore-Penrose pseudo inverse, and Frobenius norm are denoted as
$\left(\cdot\right)^{T},\left(\cdot\right)^{\dagger},\left\|\cdot\right\|_{F}$
respectively. Matlab notations will be used to denote submatrices of a tensor.
For instance, we use $\boldsymbol{\mathcal{T}}\left(:,:,k\right)$ to denote
the frontal slice of a tensor by fixing the third index to $k$.

The mode-$n$ vectors of $\boldsymbol{\mathcal{T}}$ are obtained by fixing all
but the $n$th index of $\boldsymbol{\mathcal{T}}$. For a given matrix
$\mathbf{T}\in\mathbb{R}^{I\times J}$, ${\rm
vec}\left(\mathbf{T}\right)\triangleq\left[t_{1}^{T},\cdots,t_{I}^{T}\right]\in\mathbb{R}^{IJ}$
denotes column-wise vectorization of $\mathbf{T}$ and ${\rm
unvec}\left(\cdot\right)$ performs the inverse. The mode-$n$ matricization of
a third-order tensor $\boldsymbol{\mathcal{T}}\in\mathbb{R}^{I\times J\times
K}$ is denoted as $\mathbf{T}_{n},n=1,2,3$, and defined by:

$(\mathbf{T}_{1})_{i,(j-1)K+k}=(\mathbf{T}_{2})_{j,(i-1)K+k}=(\mathbf{T}_{3})_{k,(i-1)J+j}=t_{i,j,k}.$

A polyadic decomposition (PD) of $\boldsymbol{\mathcal{T}}$ expresses
$\boldsymbol{\mathcal{T}}$ as the sum of rank-1 terms:

$\boldsymbol{\mathcal{T}}=\left[\kern-1.49994pt\left[\mathbf{A},\mathbf{B},\mathbf{C}\right]\kern-1.49994pt\right]_{R}=\sum_{r=1}^{R}\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r}\in\mathbb{R}^{I\times
J\times K}.$

where
$\mathbf{A}\triangleq\left[\mathbf{a}_{1},\cdots,\mathbf{a}_{R}\right]\in\mathbb{R}^{I\times
R}$,
$\mathbf{B}\triangleq\left[\mathbf{b}_{1},\cdots,\mathbf{b}_{R}\right]\in\mathbb{R}^{J\times
R}$, and
$\mathbf{C}\triangleq\left[\mathbf{c}_{1},\cdots,\mathbf{c}_{R}\right]\in\mathbb{R}^{K\times
R}$. We call it a canonical PD (CPD) if $R$ is minimal, and the $R$ is the
rank of $\boldsymbol{\mathcal{T}}$.

##  II Problem Formulation

We consider the computation of an approximate CPD of a large-scale third-order
tensor $\boldsymbol{\mathcal{T}}$ of size $I\times J\times K$. The term
“approximate” indicates that the CPD model does not hold precisely for
$\boldsymbol{\mathcal{T}}$, due to practical model defect such as model
mismatch or noise, and we formulate this model defect as an additive noise
term to an exact CPD model. That is to say, the data tensor
$\boldsymbol{\mathcal{T}}$ is formulated as a structured tensor
$\boldsymbol{\mathcal{S}}$ plus an unstructured noise tensor
$\boldsymbol{\mathcal{N}}$:

$\boldsymbol{\mathcal{T}}=\boldsymbol{\mathcal{S}}+\boldsymbol{\mathcal{N}},$
(1)

where $\boldsymbol{\mathcal{S}}$ admits a third-order CPD of rank $R$:

$\boldsymbol{\mathcal{S}}=\left[\kern-1.49994pt\left[\mathbf{A},\mathbf{B},\mathbf{C}\right]\kern-1.49994pt\right]_{R}.$
(2)

with $\mathbf{A}\in\mathbb{R}^{I\times R},\mathbf{B}\in\mathbb{R}^{J\times
R},\mathbf{C}\in\mathbb{C}^{K\times R}$ being the first, second, and third
factor matrices, respectively.

We assume that the data tensor has large size and low rank, that is to say,
${\rm min}(I,J,K)\gg R$. We also assume that the factor matrices
$\mathbf{A},\mathbf{B},\mathbf{C}$ all have full column rank, and the
approximate CPD problem in this case is labelled as overdetermined. Note that
we do not consider the underdetermined case where one or more factor matrices
do not have full column rank. We can see from (2) that
$\boldsymbol{\mathcal{S}}$ has exact multilinear ${\rm rank}\ (R,R,R)$, i.e.,
$R_{n}={\rm rank}\left(\mathbf{S}_{n}\right)=R,n=1,2,3$.

As such, the approximate CPD model (1) can be written in the following Tucker
form:

$\boldsymbol{\mathcal{T}}=\boldsymbol{\mathcal{I}}\times_{1}\mathbf{A}\times_{2}\mathbf{B}\times_{3}\mathbf{C}+\boldsymbol{\mathcal{N}}.$
(3)

The problem of approximate CPD can be formulated as the following least
squares (LS) based optimization problem:

$\left\\{\tilde{\mathbf{A}},\tilde{\mathbf{B}},\tilde{\mathbf{C}}\right\\}=\underset{\mathbf{A},\mathbf{B},\mathbf{C}}{\arg\min}\left(\boldsymbol{\mathcal{T}}-\boldsymbol{\mathcal{I}}\times_{1}\mathbf{A}\times_{2}\mathbf{B}\times_{3}\mathbf{C}\right).$
(4)

where $\tilde{\mathbf{A}}$, $\tilde{\mathbf{B}}$, $\tilde{\mathbf{C}}$ denote
the estimates of $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, respectively.

##  III  Proposed Algorithm

###  III-A CPD Based on Tensor Compression

Let us assume that the structured part $\boldsymbol{\mathcal{S}}$ in the data
tensor has exact rank $R$, and hence multilinear ${\rm rank}\ (R,R,R)$, as
explained in Section II. We construct unitary matrices
$\mathbf{U},\mathbf{V},\mathbf{W}$. The columns of each matrix form a set of
orthonormal bases of the vector space in corresponding mode.

In practice, noise is usually present and we do not directly know
$\boldsymbol{\mathcal{S}}$. In this case, $\mathbf{U},\mathbf{V},\mathbf{W}$
need to be estimated from the noisy tensor $\boldsymbol{\mathcal{T}}$. Methods
for estimating these matrices will be briefly explained later. Note that we
often let $\mathbf{U},\mathbf{V},\mathbf{W}$ have more than $R$ columns to
oversample the vector space of each tensor mode, such that the major structure
of the tensor in that mode is well preserved. We denote the number of columns
of $\mathbf{U},\mathbf{V},\mathbf{W}$ as $R^{\prime}$, $R^{\prime}>R$.

We can use $\mathbf{U},\mathbf{V},\mathbf{W}$ to project the tensor
$\boldsymbol{\mathcal{T}}$ in all modes into a compressed core tensor
$\boldsymbol{\mathcal{G}}\in\mathbb{R}^{R^{\prime}\times R^{\prime}\times
R^{\prime}}$ as follows:

$\boldsymbol{\mathcal{G}}\triangleq\boldsymbol{\mathcal{T}}\times_{1}\mathbf{U}^{T}\times_{2}\mathbf{V}^{T}\times_{3}\mathbf{W}^{T}.$
(5)

Substituting (3) into (5) yields:

$\boldsymbol{\mathcal{G}}=\boldsymbol{\mathcal{I}}\times_{1}\mathbf{A^{\prime}}\times_{2}\mathbf{B^{\prime}}\times_{3}\mathbf{C}^{\prime}+\boldsymbol{\mathcal{N}}^{\prime}.$
(6)

where
$\mathbf{A^{\prime}}\triangleq\mathbf{U}^{T}\mathbf{A}\in\mathbb{R}^{R^{\prime}\times
R},\mathbf{B^{\prime}}\triangleq\mathbf{V}^{T}\mathbf{B}\in\mathbb{R}^{R^{\prime}\times
R},\mathbf{C^{\prime}}\triangleq\mathbf{W}^{T}\mathbf{C}\in\mathbb{R}^{R^{\prime}\times
R},$ and
$\boldsymbol{\mathcal{N}}^{\prime}=\boldsymbol{\mathcal{N}}\times_{1}\mathbf{U}^{T}\times_{2}\mathbf{V}^{T}\times_{3}\mathbf{W}^{T}$.

We note that the core tensor $\boldsymbol{\mathcal{G}}$ also admits an
approximate CPD but with smaller size than $\boldsymbol{\mathcal{T}}$.
Therefore, by computing the approximate CPD of $\boldsymbol{\mathcal{G}}$ we
obtain estimates of its factor matrices
$\mathbf{A}^{\prime},\mathbf{B}^{\prime}$ and $\mathbf{C}^{\prime}$. Then the
factor matrices $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$ can be computed as:

$\mathbf{A}=\mathbf{U}\mathbf{A^{\prime}},\mathbf{B}=\mathbf{V}\mathbf{B^{\prime}},\mathbf{C}=\mathbf{W}\mathbf{C^{\prime}}.$
(7)

We note that the above idea of compression based CPD has been adopted in
several existing works such as the classical complex parallel factor analysis
(COMFAC) algorithm[17] and the CPD function in Tensorlab 3.0 software
package[18]. In these works, the projection matrices $\mathbf{U},\mathbf{V}$
and $\mathbf{W}$ are computed via low multilinear rank approximation (LMLRA),
which is usually computed iteratively, i.e., via the higher-order orthogonal
iteration. However, a truncated multilinear singular value decomposition
(MLSVD) may still offer near-to-optimal results and thus is often adopted to
compute LMLRA.

###  III-B Random Projection for Tensor Compression

For a large-scale tensor with low rank, the LMLRA based tensor compression
becomes computationally prohibited. Therefore, RAP based methods were
proposed. We note that RAP was originally used to compute the SVD of a large-
scale matrix[19, 20] and was later adapted for the computation of various
large-scale tensor models[3, 4, 5, 6, 7]. Here we explain RAP with power
iterations proposed in [4].

First, a sketch matrix is constructed that holds a set of approximate base
vectors of the vector space of each mode of the tensor as follows (we take the
construction of the sketch matrix in the first mode as an example):

$\mathbf{X}^{(m)}=\mathbf{T}_{1}^{(m)}\mathbf{\Psi}_{1}.$ (8)

where $\mathbf{T}_{1}^{(m)}$ is defined as:
$\mathbf{T}^{\left(m\right)}_{1}\triangleq\left(\mathbf{T}_{1}\mathbf{T}_{1}^{T}\right)^{m}\mathbf{T}_{1}$,
$\mathbf{T}_{1}$ is the mode-1 matricization of tensor
$\boldsymbol{\mathcal{T}}$ and $\mathbf{\Psi}_{1}$ is a random test matrix
used to sample the mode-1 vector space of $\boldsymbol{\mathcal{T}}$, i.e.,
the column space of $\mathbf{T}_{1}$. Note that the randomly generated columns
of $\mathbf{\Psi}_{1}$ are linearly independent in the generic sense, and that
$\mathbf{T}_{1}^{(m)}$ and $\mathbf{T}_{1}$ have identical column space.
Hence, the set of column vectors
$\left\\{\mathbf{x}_{1},\cdots,\mathbf{x}_{R^{\prime}}\right\\}$ of
$\mathbf{X}^{(m)}$ provides an efficient sample of the mode-1 vector space of
$\boldsymbol{\mathcal{T}}$. In the noiseless case, an economic QR
decomposition of the sketch matrix $\mathbf{X}^{(m)}$ provides a set of
orthonormal bases of that vector space and thus the projection matrix
$\mathbf{U}^{(m)}$ is constructed using these basis vectors as columns. When
noise is present, we can perform a truncated SVD(T-SVD) of the sketch matrix
to estimate $\mathbf{U}^{(m)}$. Note that we use $\mathbf{T}_{1}^{(m)}$
instead of $\mathbf{T}_{1}$ to construct $\mathbf{X}^{(m)}$ because the
singular values of $\mathbf{T}_{1}^{(m)}$ are
$\left\\{\sigma_{r}^{2m+1}\right\\}$ with $\left\\{\sigma_{r}\right\\}$ being
the singular values of $\mathbf{T}_{1}$ such that the singular vectors of
$\mathbf{T}_{1}^{(m)}$ associated with major singular values are more
identifiable.

In practice, the sketch matrices are usually constructed via normalized
subspace iterations[4]. After the construction of the sketch matrices,
$\mathbf{X}^{(m)}$ and $\mathbf{Y}^{(m)}$, $\mathbf{Z}^{(m)}$ that correspond
to the first, second, and third mode of the tensor, respectively, the
projection matrices $\mathbf{U}^{(m)},\mathbf{V}^{(m)},\mathbf{W}^{(m)}$ are
computed via the T-SVD of $\mathbf{X}^{(m)}$ and $\mathbf{Y}^{(m)}$,
$\mathbf{Z}^{(m)}$ , respectively.

###  III-C Coupled Random Projection Based CPD

Now we propose a CPD algorithm for large-scale tensors based on coupled RAP
(CoRAP). First, we construct sketch matrices $\mathbf{X}^{(m)}$ and
$\mathbf{Y}^{(m)}$ in the first and second mode, respectively, according to
(8), where $m$ is a fixed integer denoting the order of power iterations.
Performing the T-SVD of $\mathbf{X}^{(m)}$ and $\mathbf{Y}^{(m)}$ yields the
projection matrices $\mathbf{U}^{(m)}\in\mathbb{R}^{I\times R^{\prime}}$ and
$\mathbf{V}^{(m)}\in\mathbb{R}^{I\times R^{\prime}}$, respectively.

Then we let $m$ vary and repeat the above procedure for each $m$. As such, we
obtain two sets of projection matrices
$\left\\{\mathbf{U}^{(m)}\right\\},\left\\{\mathbf{V}^{(m)}\right\\}$ in the
first and second mode, respectively, where $m=1,...,M$ and $M$ is the maximal
order of power iterations. In the third mode, we construct one sketch matrix
$\mathbf{Z}$ and calculate the projection matrix
$\mathbf{W}\in\mathbb{R}^{K\times R^{\prime}}$.

Finally, we obtain $M$ projection matrices in the first mode and the second
mode, respectively, as well as one projection matrix in the third mode. We
group all the projection matrices into $M$ triads:
$\boldsymbol{\Omega}^{(m)}\triangleq\left\\{\mathbf{U}^{(m)},\mathbf{V}^{(m)},\mathbf{W}\right\\},\
m=1,...,M$. Each triad $\boldsymbol{\Omega}^{(m)}$ contains projection
matrices in all three modes, while different triads share a common projection
matrix $\mathbf{W}$ in the third mode.

We perform the following projection using each $\boldsymbol{\Omega}^{(m)}$
(5):

$\boldsymbol{\mathcal{G}}^{(m)}\triangleq\boldsymbol{\mathcal{T}}\times_{1}\mathbf{U}^{(m)T}\times_{2}\mathbf{V}^{(m)T}\times_{3}\mathbf{W}^{T},$
(9)

According to (6), each $\boldsymbol{\mathcal{G}}^{(m)}$ admits an approximate
CPD:

$\boldsymbol{\mathcal{G}}^{(m)}=\left[\kern-2.5pt\left[\mathbf{A}^{(m)},\mathbf{B}^{(m)},\mathbf{C}^{\prime}\right]\kern-2.5pt\right]+\boldsymbol{\mathcal{N}}^{(m)}.$
(10)

where
$\,\mathbf{A}^{(m)}\triangleq\mathbf{U}^{(m)T}\mathbf{A},\mathbf{B}^{(m)}\triangleq\mathbf{R}^{(m)T}\mathbf{B},\mathbf{C}^{\prime}\triangleq\mathbf{W}^{T}\mathbf{A}$
and
$\,\boldsymbol{\mathcal{N}}^{(m)}\triangleq\boldsymbol{\mathcal{N}}\times_{1}\mathbf{U}^{(m)T}\times_{2}\mathbf{V}^{(m)T}\times_{3}\mathbf{W}^{T}$.

Note that the set of tensors $\\{\boldsymbol{\mathcal{G}}^{(m)},m=1,...,M\\}$
together admits an approximate C-CPD, with a common factor matrix
$\mathbf{C}^{\prime}$ in the third mode. Therefore, instead of performing the
uncoupled CPD for each tensor $\boldsymbol{\mathcal{G}}^{(m)}$, we can perform
a C-CPD for all the tensors $\\{\boldsymbol{\mathcal{G}}^{(m)}\\}$. The
projections (9) for all $m=1,...,M$ are together labelled as CoRAP.

As long as $\mathbf{A}^{(m)}$, $\mathbf{B}^{(m)}$, and $\mathbf{C}^{\prime}$
are computed by C-CPD of $\\{\boldsymbol{\mathcal{G}}^{(m)}\\}$, $m=1,...,M$,
we recover the factor matrices $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$ by
(7). Note that as $m$ varies from 1 to $M$, we indeed have $M$ estimates of
$\mathbf{A}$ and $\mathbf{B}$, respectively. Here we select the optimal value
of $m$, $m_{opt}$, such that
$\left[\kern-2.5pt\left[\mathbf{A}^{(m_{opt})},\mathbf{B}^{(m_{opt})},\mathbf{C}^{\prime}\right]\kern-2.5pt\right]$
offers the best fit of the data tensor in the least squares sense:

$m_{opt}=\underset{m}{\arg\min}\left\|\boldsymbol{\mathcal{G}}^{(m)}-\left[\kern-2.5pt\left[\mathbf{A}^{(m)},\mathbf{B}^{(m)},\mathbf{C}^{\prime}\right]\kern-2.5pt\right]\right\|^{2}_{F}.$
(11)

The factor matrices of the original tensor $\boldsymbol{\mathcal{T}}$ are then
computed as:

$\,\mathbf{A}\\!=\\!\mathbf{U}^{(m_{opt})}\mathbf{A}^{(m_{opt})},\mathbf{B}\\!=\\!\mathbf{V}^{(m_{opt})}\mathbf{B}^{(m_{opt})},\mathbf{C}\\!=\\!\mathbf{W}\mathbf{C}^{\prime}.$
(12)

###  III-D Coupled CPD

As explained in Subsection III.C, C-CPD plays an important role in CoRAP based
CPD. In fact, C-CPD algorithms have been studied tremendously in the past
decade, including the algebraic algorithms[10, 12, 11] and optimization based
algorithms[21, 11], which are shown to have their own pros and cons,
respectively. In particular, the algebraic algorithms are guaranteed to return
the exact solution in the noiseless case. However, they only return suboptimal
solution when noise is present. On the other hand, the optimization based
algorithms can provide optimal results in the LS sense if they converge to the
global minimum. However, they are likely to converge to a local minimum and
are thus sensitive to initialization. As such, we can use the algebraic
algorithms to efficiently initialize the optimization based algorithms.

We note that the C-CPD problem in our setting (10) is overdetermined in the
sense that all the factor matrices have full column rank. Therefore, instead
of using the algebraic algorithms in [10] and [12] that are designed for
underdetermined problems, we can use a simple algebraic C-CPD algorithm,
summarized as follows:

  * 1)

Calculate the CPD of one tensor $\boldsymbol{\mathcal{G}}^{(m)}$ with fixed
$m$. We assume without loss of generality that CPD of
$\boldsymbol{\mathcal{G}}^{(1)}$ is performed to compute estimates of
$\mathbf{A}^{(1)}$, $\mathbf{B}^{(1)}$, and $\mathbf{C}^{\prime}$;

  * 2)

Since $\mathbf{C}^{\prime}$ is already obtained, we have:

$\mathbf{G}_{3}^{(m)}\left(\mathbf{C}^{\prime
T}\right)^{\dagger}=\mathbf{A}^{(m)}\odot\mathbf{B}^{(m)}+\mathbf{N}^{\prime(m)},$
(13)

where $m\in\left[2,M\right]$, and $\mathbf{N}^{\prime(m)}$ is the noise term.
It is shown in (16) that each column of
$\mathbf{G}_{3}^{(m)}\left(\mathbf{C}^{\prime T}\right)^{\dagger}$, denoted as
$\mathbf{g}_{r}^{(m)}$, is approximately a vectorized rank-1 matrix:

${\rm
unvec}(\mathbf{g}_{r}^{(m)})\approx\mathbf{a}_{r}^{(m)}\cdot\mathbf{b}_{r}^{(m)T},\
r=1,...,R.$ (14)

Therefore, by performing rank-1 approximation of ${\rm
unvec}(\mathbf{g}_{r}^{(m)})$ with varying $r$ and fixed $m$ we obtain factor
matrices $\mathbf{A}^{(m)}$ and $\mathbf{B}^{(m)}$. We let $m$ vary and repeat
the above procedure (13) and (14) for each $m$, then all the factor matrices
can be computed.

Now that we obtain the initial results via an algebraic algorithm, we use them
as initialization of an optimization based C-CPD algorithm to obtain the final
results. Possible candidates include the alternating least squares (ALS)[11]
and nonlinear least squares (NLS)[22] algorithms.

##  IV Numerical Experiments

In this section we use numerical experiments to illustrate the performance of
the proposed CoRAP based CPD algorithm in comparison with RAP based CPD [4]
and direct CPD. In all the experiments the data tensor is generated as
follows:

$\boldsymbol{\mathcal{T}}=P_{s}\boldsymbol{\mathcal{S}}+P_{n}\boldsymbol{\mathcal{N}}_{n},$
(15)

where tensor $\boldsymbol{\mathcal{S}}$ is constructed according to (2) with
the factor matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$ randomly drawn
from a Gaussion distribution with zero mean and unit variance.
$\boldsymbol{\mathcal{N}}$ denotes the additive noise term that is also drawn
from a Gaussion distribution with zero mean and unit variance. Parameters
$P_{s}$ and $P_{n}$ denote signal level and noise level, respectively. The
Signal-to-noise ratio(SNR) is defined as:

$\rm SNR\triangleq 10\rm
log_{10}(\textit{P}_{\textit{s}}/\textit{P}_{\textit{n}}),$ (16)

We evaluate the performance of compared algorithms in terms of mean relative
error, defined as follows:

$\epsilon\triangleq\frac{1}{3}\sum_{\mathbf{H}\in\\{\mathbf{A},\mathbf{B},\mathbf{C}\\}}\left\|\mathbf{H}\right\|_{F}^{-2}\left\|\mathbf{H}-\tilde{\mathbf{H}}\mathbf{P}\mathbf{S}\right\|_{F}^{2}.$
(17)

where $\mathbf{H}$ and $\tilde{\mathbf{H}}$ denote the true factor matrix and
its estimate, respectively. $\mathbf{P}$ and $\mathbf{S}$ are permutation and
scaling matrices, respectively, such that the estimated factor matrix
$\tilde{\mathbf{H}}$ is permuted and scaled to optimally fit $\mathbf{H}$ in
the LS sense. In addition, we use CPU time to evaluate how fast the algorithms
completes the calculation. The mean relative error and CPU time are calculated
as the average of the results of 200 Monte Carlo runs. Note that the direct
CPD is implemented using ‘cpd.m’ function in Tensorlab 3.0 with default
setting[18].

Figure 1: The performance of CoRAP-CPD, RAP-CPD and direct CPD versus SNR of a
noisy tensor $\boldsymbol{\mathcal{T}}$ of size $500\times 500\times 500$ and
rank $R=50$. The left sub-figure shows the overall error, and the right shows
the CPU time

In the first experiment, we mainly study the impact of noise on the compared
algorithms. The data tensor $\boldsymbol{\mathcal{S}}$ is of size $500\times
500\times 500$ with rank $R=50$. For the proposed algorithm, we set $M=2$. In
addition, for the implementation of C-CPD in CoRAP-CPD, we use algebraic C-CPD
to efficiently initialize an optimization based C-CPD algorithm. Experiment
results, as demonstrated in Fig. 1, show that the performance of CoRAP-CPD, in
terms of mean relative error, is better than that of RAP-CPD, and it
approaches that of direct CPD as SNR increases. On the other hand, the CPU
time of CoRAP-CPD is slightly higher than that of RAP-CPD, and is still lower
than that of direct CPD. This is because more data tensors are generated and
processed for CoRAP than RAP. The observation in this experiment generally
suggests that CoRAP based CPD improves the accuracy of its RAP based
counterpart, while can still maintain the computational efficiency of RAP, in
comparison with direct CPD.

Figure 2: Performance of CoRAP-CPD for different $M$ of a noisy tensor of size
$200\times 200\times 200$ with rank $R=20$. The left and right sub-figures
correspond to SNR = $-2$dB and SNR = 10dB, respectively.

In the second experiment, we study the impact of $M$ on the accuracy of CoRAP-
CPD. The data tensor is generated analogously to the last experiment except
that we let $M$ vary while SNR is fixed. The results for SNR = $-2$dB and SNR
= 10dB are given in Fig. 2. The observation generally suggests that $M$ has
little impact on the accuracy of CoRAP-CPD, while a larger $M$ would lead to
higher complexity. Therefore, in practice we can generally use a small $M$,
i.e.,$M=2$.

Figure 3: The performance of CoRAP-CPD, RAP-CPD and direct CPD versus rank $R$
of a noisy tensor $\boldsymbol{\mathcal{T}}$ of size $500\times 500\times
500$. The left and right sub-figures correspond to SNR = $-2$dB and SNR = 6dB,
respectively.

In the third experiment, we study the impact of tensor rank, $R$, on the
performance of CoRAP-CPD and RAP-CPD. The data tensor is generated in an
analogous manner as the last two experiments. We fix SNR to $-2$dB and 6dB,
respectively, and let $R$ vary from 10 to 50. The results are given in Fig. 3.
It has been observed that in both low SNR (SNR = $-2$dB) and medium SNR (SNR =
6dB), the proposed CoRAP based CPD algorithm has remarkably better accuracy
than RAP-CPD when $R$ is large. This observation clearly shows the interests
of the proposed algorithm in difficult cases where the data tensor not only
has large size but also relatively high rank.

##  V Conclusion

A novel CPD algorithm for large-scale tensor is proposed by combining the idea
of random projection (RAP) and coupled tensor decomposition. The coupled
random projection (CoRAP) technique is developed and applied to project the
original large-scale tensor into a set of small tensors that together admits a
coupled CPD(C-CPD), and a C-CPD algorithm is used to jointly decompose these
tensors. Then, the results of C-CPD are back projected to obtain the factor
matrices of the original tensor. Through numerical experiments, we have shown
that the proposed CoRAP based CPD algorithm has better performance in terms of
accuracy than RAP based CPD, with slightly increased CPU time. In particular,
the proposed algorithm has interests in difficult cases where the data tensor
not only has large size but also relatively high rank.

## References

  * [1] N. Vervliet and L. De Lathauwer,  “A randomized block sampling approach to canonical polyadic decomposition of large-scale tensors,”  IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 2, pp. 284–295, 2015.
  * [2] A. H. Phan and A. Cichocki,  “PARAFAC algorithms for large-scale problems,”  Neurocomputing, vol. 74, no. 11, pp. 1970–1984, 2011.
  * [3] N. D. Sidiropoulos, E. E. Papalexakis, and C. Faloutsos,  “Parallel randomly compressed cubes: A scalable distributed architecture for big tensor decomposition,”  IEEE Signal Processing Magazine, vol. 31, no. 5, pp. 57–70, 2014\.
  * [4] N. B. Erichson, K. Manohar, S. L. Brunton, and J. N. Kutz,  “Randomized CP tensor decomposition,”  Machine Learning: Science and Technology, vol. 1, no. 2, pp. 025012, 2020.
  * [5] L.-H. Yuan, C. Li, J.-T Cao, and Q.-B. Zhao,  “Randomized tensor ring decomposition and its application to large-scale data reconstruction,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Brighton, U. K., May.12–17, 2019, pp. 2127–2131.
  * [6] M.-L Che and Y.-M Wei,  “Randomized algorithms for the approximations of Tucker and the tensor train decompositions,”  Advances in Computational Mathematics, vol. 45, no. 1, pp. 395–428, 2019.
  * [7] R. Minster, A. K. Saibaba, and M. E. Kilmer,  “Randomized algorithms for low-rank tensor decompositions in the Tucker format,”  SIAM Journal on Mathematics of Data Science, vol. 2, no. 1, pp. 189–215, 2020.
  * [8] E. Acar, R. Bro, and A. K. Smilde,  “Data fusion in metabolomics using coupled matrix and tensor factorizations,”  Proceedings of the IEEE, vol. 103, no. 9, pp. 1602–1620, 2015.
  * [9] M. Sørensen and L. De Lathauwer,  “Coupled canonical polyadic decompositions and (coupled) decompositions in multilinear rank-(${L}_{r,n},{L}_{r,n},1$) terms—part I: Uniqueness,”  SIAM Journal on Matrix Analysis and Applications, vol. 36, no. 2, pp. 496–522, 2015.
  * [10] M. Sørensen, I. Domanov, and L. De Lathauwer,  “Coupled canonical polyadic decompositions and (coupled) decompositions in multilinear rank-(${L}_{r,n},{L}_{r,n},1$) terms—part II: Algorithms,”  SIAM Journal on Matrix Analysis and Applications, vol. 36, no. 3, pp. 1015–1045, 2015.
  * [11] X.-F. Gong, Q.-H. Lin, F.-Y. Cong, and L. De Lathauwer,  “Double coupled canonical polyadic decomposition for joint blind source separation,”  IEEE Transactions on Signal Processing, vol. 66, no. 13, pp. 3475–3490, 2018.
  * [12] X.-F. Gong, Q.-H. Lin, F.-Y. Cong, and L. De Lathauwer,  “Double coupled canonical polyadic decomposition of third-order tensors: Algebraic algorithm and relaxed uniqueness conditions,”  Signal Processing: Image Communication, vol. 73, pp. 22–36, 2019\.
  * [13] X.-F. Gong, Y.-N. Hao, and Q.-H. Lin,  “Joint canonical polyadic decomposition of two tensors with one shared loading matrix,”  in Proceedings of IEEE International Workshop on Machine Learning for Signal Processing (MLSP), Southampton, U. K., Sep.22–25, 2013, pp. 1–6.
  * [14] R. C. Farias, J. E. Cohen, and P. Comon,  “Exploring multimodal data fusion through joint decompositions with flexible couplings,”  IEEE Transactions on Signal Processing, vol. 64, no. 18, pp. 4830–4844, 2016.
  * [15] X.-F Gong, Q. H Lin, O. Debals, N. Vervliet, and L. De Lathauwer,  “Coupled rank-(${L}_{m},{L}_{n},\cdot$) block term decomposition by coupled block simultaneous generalized Schur decomposition,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Shanghai, China, Mar.20–25, 2016, pp. 2554–2558.
  * [16] J.-X. Yang, X.-F. Gong, H. Li, Y.-G. Xu, and Z.-W. Liu,  “Using coupled multilinear rank-(${L},{L},1$) block term decomposition in multi-static-multi-pulse MIMO radar to localize targets,”  in Proceedings of International Symposium on Neural Networks, Moscow, Russia, Jul.10–12, 2019, pp. 565–574.
  * [17] R. Bro, N. D. Sidiropoulos, and G. B. Giannakis,  “A fast least squares algorithm for separating trilinear mixtures,”  in Proceedings of International workshop on independent component analysis and signal separation, Aussois, France, Jan.11–15, 1999, pp. 11–15.
  * [18] N. Vervliet, O. Debals, L. Sorber, M. Van Barel, and L. De Lathauwer,  “Tensorlab 3.0,” [Online], Mar. 2016,  Available: https://www.tensorlab.net.
  * [19] F. Woolfe, E. Liberty, V. Rokhlin, and M. Tygert,  “A fast randomized algorithm for the approximation of matrices,”  Applied and Computational Harmonic Analysis, vol. 25, no. 3, pp. 335–366, 2008.
  * [20] N. B. Erichson, S. Voronin, S. L. Brunton, and J. N. Kutz,  “Randomized matrix decompositions using R,”  Journal of Statistical Software, vol. 89, no. 11, pp. 1–48, 2019\.
  * [21] L. Sorber, M. Van Barel, and L. De Lathauwer,  “Structured data fusion,”  IEEE Journal of Selected Topics in Signal Processing, vol. 9, no. 4, pp. 586–600, 2015.
  * [22] L. Sorber, M. Van Barel, and L. De Lathauwer,  “Optimization-based algorithms for tensor decompositions: Canonical polyadic decomposition, decomposition in rank-(${L}_{r},{L}_{r},1$) terms, and a new generalization,”  SIAM Journal on Optimization, vol. 23, no. 2, pp. 695–720, 2013\.


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.††thanks: These authors contributed
equally to this work.††thanks: Corresponding authors:
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Corresponding authors:
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Corresponding authors:
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
# Evolution of Band Structure in a Kagome Superconductor
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$: Toward Universal
Understanding of CDW and Superconducting Phase Diagrams

Shuto Suzuki  Department of Physics, Graduate School of Science, Tohoku
University, Sendai 980-8578, Japan  Takemi Kato  Advanced Institute for
Materials Research (WPI-AIMR), Tohoku University, Sendai 980-8577, Japan
Yongkai Li  Centre for Quantum Physics, Key Laboratory of Advanced
Optoelectronic Quantum Architecture and Measurement (MOE), School of Physics,
Beijing Institute of Technology, Beijing 100081, P. R. China  Beijing Key Lab
of Nanophotonics and Ultrafine Optoelectronic Systems, Beijing Institute of
Technology, Beijing 100081, P. R. China  Material Science Center, Yangtze
Delta Region Academy of Beijing Institute of Technology, Jiaxing, 314011, P.
R. China  Kosuke Nakayama  Department of Physics, Graduate School of Science,
Tohoku University, Sendai 980-8578, Japan  Zhiwei Wang  Centre for Quantum
Physics, Key Laboratory of Advanced Optoelectronic Quantum Architecture and
Measurement (MOE), School of Physics, Beijing Institute of Technology, Beijing
100081, P. R. China  Beijing Key Lab of Nanophotonics and Ultrafine
Optoelectronic Systems, Beijing Institute of Technology, Beijing 100081, P. R.
China  Material Science Center, Yangtze Delta Region Academy of Beijing
Institute of Technology, Jiaxing, 314011, P. R. China  Beijing Institute of
Technology, Zhuhai, 519000, P. R. China  Seigo Souma  Advanced Institute for
Materials Research (WPI-AIMR), Tohoku University, Sendai 980-8577, Japan
Center for Science and Innovation in Spintronics (CSIS), Tohoku University,
Sendai 980-8577, Japan  Kenichi Ozawa  Institute of Materials Structure
Science, High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki
305-0801, Japan  Miho Kitamura  Institute of Materials Structure Science, High
Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801,
Japan  National Institutes for Quantum Science and Technology (QST), Sendai
980-8579, Japan  Koji Horiba  National Institutes for Quantum Science and
Technology (QST), Sendai 980-8579, Japan  Hiroshi Kumigashira  Institute of
Multidisciplinary Research for Advanced Materials (IMRAM), Tohoku University,
Sendai 980-8577, Japan  Takashi Takahashi  Department of Physics, Graduate
School of Science, Tohoku University, Sendai 980-8578, Japan  Yugui Yao
Centre for Quantum Physics, Key Laboratory of Advanced Optoelectronic Quantum
Architecture and Measurement (MOE), School of Physics, Beijing Institute of
Technology, Beijing 100081, P. R. China  Beijing Key Lab of Nanophotonics and
Ultrafine Optoelectronic Systems, Beijing Institute of Technology, Beijing
100081, P. R. China  Takafumi Sato  Department of Physics, Graduate School of
Science, Tohoku University, Sendai 980-8578, Japan  Advanced Institute for
Materials Research (WPI-AIMR), Tohoku University, Sendai 980-8577, Japan
Center for Science and Innovation in Spintronics (CSIS), Tohoku University,
Sendai 980-8577, Japan  International Center for Synchrotron Radiation
Innovation Smart (SRIS), Tohoku University, Sendai 980-8577, Japan
Mathematical Science Center for Co-creative Society (MathCCS), Tohoku
University, Sendai 980-8578, Japan

###### Abstract

Kagome superconductors $A\textrm{V}_{3}\textrm{Sb}_{5}$ ($A=$ K, Rb, Cs)
exhibit a characteristic superconducting and charge-density wave (CDW) phase
diagram upon carrier doping and chemical substitution. However, the key
electronic states responsible for such a phase diagram have yet to be
clarified. Here we report a systematic micro-focused angle-resolved
photoemission spectroscopy (ARPES) study of
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$ as a function of Cr
content $x$, where Cr substitution causes monotonic reduction of
superconducting and CDW transition temperatures. We found that the V-derived
bands forming saddle points at the M point and Dirac nodes along high-symmetry
cuts show an energy shift due to electron doping by Cr substitution, whereas
the Sb-derived electron band at the $\Gamma$ point remains almost unchanged,
signifying an orbital-selective band shift. We also found that band doubling
associated with the emergence of three-dimensional CDW identified at $x=0$
vanishes at $x\geq 0.25$, in line with the disappearance of CDW. A comparison
of band diagrams among Ti-, Nb-, and Cr-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ suggests the importance to
simultaneously take into account the two saddle points at the M point and
their proximity to the Fermi energy, to understand the complex phase diagram
against carrier doping and chemical pressure.

Two-dimensional (2D) kagome lattice provides a fertile ground for studying
exotic quantum phenomena originating from correlated and topological band
structures, featuring a Dirac-cone band, a flat band, and a saddle-point (SP)
band [1]. Among these bands, the SP band forms a van Hove singularity, leading
to a divergence in the density of states (DOS). This enhances many-body
interactions and Fermi-surface (FS) instability, inducing various ordered
phases, such as unconventional density waves, superconductivity, and magnetism
[2, 3, 4, 5]. However, achieving such phases requires a tuning of the SP band
close to the Fermi energy ($E_{\textrm{F}}$), which has been challenging in
kagome lattices.

Recently, $A\textrm{V}_{3}\textrm{Sb}_{5}$ ($A=$ K, Rb, Cs) [Fig. 1(a)] has
been identified as the first kagome-lattice family with electron filling near
the SP band [6]. Subsequent experiments have confirmed the coexistence of
multiple SP bands near $E_{\textrm{F}}$. $A\textrm{V}_{3}\textrm{Sb}_{5}$
exhibits various competing or intertwined orders including CDW,
superconductivity, pair density wave, nematicity, and time-reversal-symmetry
breaking [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. The relationship of
these ordered states to the unique band structure in
$A\textrm{V}_{3}\textrm{Sb}_{5}$ is a matter of intensive debate.

A key strategy to address this issue is to apply perturbations (e.g.,
pressure, elemental substitution, and carrier doping [19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]) that
presumably tune the band structure and the ordered states, to clarify their
relationship. Recent experiments on $\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$
revealed intriguing phase diagrams resulting from such perturbations. For
instance, isovalent substitution of V with Nb or Ta was shown to decrease the
CDW transition temperature ($T_{\textrm{CDW}}$) while increasing the
superconducting transition temperature ($T_{\textrm{c}}$) [26, 37, 38]. In
contrast, substituting V with Ti to dope hole carriers initially decreases
$T_{\textrm{c}}$ as $T_{\textrm{CDW}}$ is lowered, and then begins to increase
$T_{\textrm{c}}$, forming a two-dome-like superconducting state [24, 25]. More
recently, $\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ samples with Cr
substitution, introducing electron carriers, was reported to show reduction of
both CDW and superconductivity [Fig. 1(b)] [31, 32]. Given these remarkable
differences in the phase diagram, an experimental study on the underlying
electronic states would provide valuable information into the key electronic
states responsible for CDW and superconductivity. However, such studies are
scarce, especially regarding the Cr-substituted case. This limits
comprehensive understanding, such as whether CDW is commonly suppressed when
SPs are shifted away from $E_{\textrm{F}}$.

In this article, we report a systematic ARPES study of Cr-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$,
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$ ($x=0,0.25$, and
0.45). We found an orbital-selective band shift upon Cr substitution,
characterized by a downward shift of V-derived SP and Dirac-cone bands while
pinning the energy position of the Sb-derived electron bands at the $\Gamma$
point. We also found that the appearance of CDW gap and band doubling is well
correlated with a reduction of $T_{\textrm{CDW}}$. We discuss implications of
the present results in relation to the mechanism of CDW and superconductivity,
and further discuss the role of chemical substitution to the band diagram
through a comparison with Ti- and Nb-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$.

High-quality single crystals of
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$ were synthesized with
the self-flux method [32]. ARPES measurements were performed using a Scienta-
Omicron DA30 spectrometer at BL-28A in Photon Factory, KEK, with a micro-
focused beam spot of $10\times 12$ $\mu$m2 [43]. Circularly polarized photons
with photon energy of $h\nu=106$ eV were used to probe the $k_{z}=0$
($\Gamma$KM) plane [44]. Samples were cleaved in-situ at low temperatures and
the surface termination (Sb vs Cs) was determined by scanning surface area and
monitoring core-level intensity as detailed below. Energy resolution was set
to be 35 meV.

Figure 1:  (a) Crystal structure of
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$. The graphic was
created by VESTA [65]. (b) CDW and superconducting phase diagram of
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$ as a function of $x$
and temperature [32, 31]. (c) Side view of crystal with cleaved surface
terminated by either Sb or Cs atoms. (d) EDCs in the Cs $4d$ core-level region
for Cs- and Sb-terminated surfaces at $x=0.45$. (e) Spatial image of the core-
level intensity on a cleaved surface for the Cs $4d$ surface peak. Red star
indicates the location where ARPES measurements shown in Figs. 2 and 3 were
carried out. (f) EDCs in the V/Cr $3p$ core-level region.

Previous scanning tunneling microscopy (STM) and spatially-resolved ARPES
studies [10, 11, 12, 13, 14, 45, 46, 47, 48, 49, 18] suggested that the
cleaved surface of $A\textrm{V}_{3}\textrm{Sb}_{5}$ contains Sb and $A$
terminations [Fig. 1(c)]. The CDW properties are markedly different between
the two terminations; the $A$ termination reflects bulk-originated CDW whereas
the Sb termination shows no indication of a CDW gap. It is thus essential to
perform surface-termination-selective ARPES measurements also for
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$. As shown in Fig.
1(d), we distinguished Sb- and Cs-terminated surfaces in the Cr-substituted
samples based on the knowledge that the spectral weight of the surface-derived
Cs core levels is enhanced for the Cs-terminated surface (red curve) relative
to the Sb-terminated one (blue curve). We carefully selected the Cs-terminated
surface by scanning the whole area of cleaved surface [see representative
mapping image in Fig. 1(e)]. To elucidate the systematic evolution upon Cr
substitution, it is also essential to clarify whether the Cr atoms indeed
exist at the cleaved surface by ARPES. As shown in Fig. 1(f), after Cr
substitution, the Cr $3p$ core-level peak appears at the binding energy
($E_{\textrm{B}}$) of $\sim$42 eV, and its intensity for $x=0.45$ becomes
nearly comparable to that of the V $3p$ peak, consistent with the fact that
almost a half of V atoms are replaced with Cr atoms.

Figure 2:  (a)–(c) ARPES intensity of
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$ along the $\Gamma$KM
cut for $x=0$, 0.25, and 0.45, respectively, measured in the normal state
($T=120$ K for $x=0$, and $T=8$ K for $x=0.25$ and 0.45) for the Cs-terminated
surface. Solid curves are a guide for the eyes to trace the band dispersion.
Horizontal dashed lines and arrows highlight the energy difference of each
band among different $x$. (d) Comparison of EDC at the $\Gamma$ point between
$x=0$ and 0.45. (e) Comparison of MDC between $x=0$ and 0.45 along cuts shown
by black lines in (a) and (c) (red and black dots, respectively), together
with the results of numerical fittings (solid curves). (f) Energy shift of Sb
$5p_{z}$, SP1, SP2, and V/Cr $3d_{xz}/d_{yz}$ bands upon Cr substitution
relative to the band position for $x$ = 0, estimated by analyzing the peak
position of EDCs.

Next we present the evolution of the band structure upon Cr substitution.
Figures 2(a)–2(c) show a comparison of the band dispersion above
$T_{\textrm{CDW}}$ among pristine ($x=0$) and Cr-substituted ($x=0.25$ and
0.45) $\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$, measured along the
$\Gamma$KM cut of bulk hexagonal Brillouin zone. At $x=0$ [Fig. 2(a)], one can
immediately recognize several dispersive bands. For example, (i) highly
dispersive electron bands at the $\Gamma$ point (red curves), (ii) linearly
dispersive bands that intersect each other midway between the $\Gamma$ and K
points (blue solid and dashed lines), (iii) a holelike band forming a SP
(called SP1) at the M point slightly above $E_{\textrm{F}}$ which intersects
the M-centered electronlike band at the K point (green curves), and (iv)
another holelike band topped at the M point slightly below $E_{\textrm{F}}$
(blue curve) which forms another SP (called SP2). These bands are attributed
to the (i) Sb $5p_{z}$ and (ii-iv) V $3d$ kagome-lattice bands (namely, green,
solid blue, and dashed blue curves represent dominant contributions from V
$3d_{xy}/d_{x^{2}-y^{2}}$, V $3d_{xz}/d_{yz}$, and another V $3d_{xz}/d_{yz}$
orbitals, respectively) [12, 50, 51, 52, 11, 44, 53, 54, 55]. All these bands
are commonly observed in the Cr-substituted samples [Figs. 2(b) and 2(c)],
indicating that the 45% Cr doping does not qualitatively alter the overall
valence-band structure. However, a careful look at Figs. 2(a)–2(c) reveals a
quantitative difference in the band energies upon Cr substitution. For
example, the crossing point associated with the V $3d_{xz}/d_{yz}$ band (blue
circles) initially located at $E_{\textrm{B}}$ = 0.15 eV at $x=0$ moves
downward and reaches 0.25 eV at $x=0.45$. This trend is also seen for the
Dirac point formed by the SP1 band and an electron band (green circles) and
also for the $\Lambda$-shaped SP2 band (blue curve) although the total band
shift is smaller than that of the V $3d_{xz}/d_{yz}$ band (blue solid and
dashed lines). Interestingly, despite such a downward shift of the V-derived
bands, the bottom of the Sb-derived band at the $\Gamma$ point is stationary
to the $x$ variation (red curves and red dashed lines). This can be seen from
the EDCs at the $\Gamma$ point for $x=0$ and 0.45 in Fig. 2(d) which signifies
no apparent difference in the peak position. On the other hand, a direct
comparison of the MDCs along a k cut crossing the lower Dirac-cone-like band
[black solid lines in Figs. 2(a) and 2(c)] in Fig. 2(e) reveals a finite
difference in the peak position between $x=0$ and 0.45. The estimated energy
shift of each band [Fig. 2(f)] demonstrates a strong orbital dependence of the
total band shift upon $x$ variation; it is nearly zero for the Sb $5p_{z}$
band, $\sim$50 meV for the SP1 and SP2 bands, and $\sim$100 meV for the V
$3d_{xz}/d_{yz}$ bands. The observed downward shift of the V-derived bands
upon increasing $x$ would naturally be explained in terms of the electron
doping associated with the replacement of group-V element (V) with group-VI
element (Cr), which affects V orbitals more than Sb orbitals.

Figure 3:  (a)–(c) FS mapping for $x=0$, 0.25, and 0.45, respectively,
measured at $T=8$ K. (d)–(f) Enlarged view in the $k$ region enclosed by a
dashed rectangle in (a)–(c), respectively. (g) Comparison of EDCs at $T=8$ K
among $x=0$, 0.25 and 0.45, obtained at the $\textbf{k}_{\mathrm{F}}$ point of
the Sb $5p_{z}$ band. (h) Same as (g) but of the SP1 band.

Now that the overall band evolution upon Cr substitution is established, next
we present the influence of CDW to the electronic structure. Figures 3(a)–3(c)
compare the FS mappings at $k_{z}\sim 0$ plane for different $x$, measured at
$T=$ 8 K (the CDW state for $x=0$ and the normal state for $x=0.25$ and 0.45).
At $x=0$ [Fig. 3(a)], there are a small circular and a large hexagonal pocket
each centered at the $\Gamma$ point, together with a triangular pocket at the
K point [56, 57, 44, 58, 59, 60, 54, 55, 53] which are assigned to the Sb
$5p_{z}$, V $3d_{xz}/d_{yz}$, and V $3d_{xy}/d_{x^{2}-y^{2}}$ orbitals,
respectively (note that a circular pocket splits into two pockets [47, 49, 61,
48]). Upon increasing $x$, the hexagonal pocket gradually shrinks whereas the
circular pocket keeps the same volume, confirming the orbital-selective band
shift. An expanded view in Fig. 3(d) for $x=0$ signifies that the hexagonal
pocket shows a doubling as indicated by black arrows, likely due to the CDW-
induced band folding along the $k_{z}$ direction [47, 49, 62, 63, 64]
accompanied by the $2\times 2\times 2$ out-of-plane unit-cell doubling [45,
56]. Intriguingly, such a doubling is absent for $x=0.25$ [Fig. 3(e)] and
$x=0.45$ [Fig. 3(f)] in the non-CDW phase, supporting that the reported band
doubling [47, 49, 62, 63, 64] is indeed of CDW origin.

We further examined the influence of CDW to the low-energy excitations by
inspecting the EDCs. Figures 3(g) and 3(h) show EDCs at $T=8$ K compared among
$x=0$, 0.25 and 0.45, obtained at the Fermi-wave-vector
($\textbf{k}_{\mathrm{F}}$) points of the Sb-$5p_{z}$ band and the SP1 band,
respectively. One can immediately recognize in Fig. 3(g) a sharp peak near
$E_{\textrm{F}}$, being cut off by the Fermi-Dirac distribution function. The
leading-edge midpoint is always located at $E_{\textrm{F}}$, indicative of the
absence of CDW gap, consistent with previous reports [44, 60, 55]. On the
other hand, at the $\textbf{k}_{\mathrm{F}}$ point of the SP1 band [Fig. 3(h)]
for $x=0$, one can recognize a hump feature at $\sim$0.1 eV which is
attributed to the large CDW gap [44, 53]. In sharp contrast, the CDW gap is
absent for $x=0.25$ and 0.45, again consistent with the disappearance of CDW.

Figure 4:  (a) Schematic of the kagome lattice where V atoms are partly
substituted with other transition metals (Ti, Nb, or Cr). (b) Band diagrams
for Sb, SP1, and SP2 bands in pristine
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$. (c) Schematic FS of
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$. (d)–(f) Same as (b), but for Ti
(hole doped), Nb (isovalent), and Cr (electron doped) substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$. Lower panels show the CDW and
superconducting phase diagrams against $x$ and temperature.

To discuss the implications of the present results in relation to the
substitution-element dependence of the electronic phase diagram, we briefly
review previous studies on CDW in pristine, Nb-substituted, and Ti-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ (Fig. 4) [25, 26, 27]. As shown in
the schematic of normal-state band diagram in Fig. 4(b), pristine
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ has the Sb band at the $\Gamma$
point, and kagome-lattice bands with SP1 and SP2 slightly above and below
$E_{\textrm{F}}$, respectively, near the M point. Based on the largest CDW gap
opening ($\sim$70 meV) on the SP1 band below $T_{\textrm{CDW}}$ ($\sim$20 meV
on the SP2 band and $\sim$0 meV on the Sb band), electron scattering between
SP1 bands at different M points [as connected by cyan arrows in Fig. 4(c)] was
suggested to be crucial for the occurrence of CDW [44]. When V is substituted
with Ti, Nb, or Cr [Fig. 4(a)], CDW is monotonically suppressed for all cases
[24, 25, 26, 31, 32]. Regarding the band structure, SP1 shifts away from
$E_{\textrm{F}}$ for both Nb and Ti substitutions [Figs. 4(d) and 4(e)]. From
this observation, the suppression of CDW was explained by the weakening of
inter-SP1 scattering [25, 26, 27]. However, the present study indicates that
considering only SP1 is not sufficient to understand the phase diagram of Cr-
substituted $\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ because CDW is still
suppressed despite SP1 moving closer to $E_{\textrm{F}}$ by electron doping
due to Cr substitution [Fig. 4(f)]. To explain these contradictory behaviors,
it is suggested that the presence of both SP1 and SP2 in close proximity to
$E_{\textrm{F}}$ is the most favorable situation for CDW (it is noted that,
similarly to the consideration of only SP1, considering only SP2 does not
explain the behavior of Nb and Ti substitution that suppresses CDW while
moving SP2 closer to $E_{\textrm{F}}$). This can be understood qualitatively
if the peaks of DOS from SP1 and SP2 are energetically close to each other and
a total DOS has a peak in the middle between the two.

Regarding superconductivity, $T_{\textrm{c}}$ increases in Nb-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ within the CDW phase but decreases
in Ti- and Cr-substituted $\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ [24, 25,
26, 31, 32]. Previous ARPES measurements of Nb-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$ have suggested that the DOS near
the $\Gamma$ point increases due to a downward energy shift of the Sb band
[Fig. 4(e)], which may result in higher $T_{\textrm{c}}$ [27]. This
explanation, based on the Sb-band DOS, seems to universally explain the trend
of $T_{\textrm{c}}$ in the CDW phase of Ti- and Cr-substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$. Specifically, hole doping in the
Ti system would lead to a decrease in the DOS around the $\Gamma$ point due to
an upward energy shift of the Sb band [Fig. 4(d)], resulting in decrease in
$T_{\textrm{c}}$ until the end point of the CDW phase. In the Cr-substituted
case, the orbital-selective electron doping prevents a clear energy shift of
the Sb band [Fig. 4(f)]. This peculiar band evolution can explain why the Cr
substitution does not result in increase in $T_{\textrm{c}}$. Nevertheless,
there is a point that cannot be explained qualitatively by considering only
the Sb band. That is, in the Ti-substituted case, $T_{\textrm{c}}$ begins to
increase once the CDW disappears, whereas $T_{\textrm{c}}$ remains zero in the
Cr-substituted case. This dichotomy may originate from (i) differences in the
orbital of the SP band staying near $E_{\textrm{F}}$ ($d_{xz}/d_{yz}$ for SP2
in the Ti-substituted case and $d_{xy}/d_{x^{2}-y^{2}}$ for SP1 in the Cr-
substituted case) and/or (ii) stronger pair-breaking effects of the Cr atoms
introduced into the kagome lattice compared to the Ti atoms.

In conclusion, we reported a systematic ARPES study of
$\textrm{Cs(V}_{1-x}\textrm{Cr}_{x})_{3}\textrm{Sb}_{5}$. In the normal state
above $T_{\textrm{CDW}}$, we found an orbital-selective band shift upon Cr
substitution, characterized by the downward energy shift of the V-derived
bands in contrast to the stationary nature of the Sb-derived band. We also
found a band doubling due to the 3D CDW only for $x=0$, whereas the doubling
is absent in non-CDW phases of $x=0.25$ and 0.45. We have concluded that it is
necessary to take into account the presence of two SP bands and their
proximity to $E_{\textrm{F}}$ to universally understand the CDW and
superconducting phase diagram of chemically substituted
$\textrm{Cs}\textrm{V}_{3}\textrm{Sb}_{5}$. The present results lay a
foundation for comprehensively understanding the role of chemical substitution
to the phase diagram and relevant electronic states in kagome superconductors.

###### Acknowledgements.

This work was supported by JST-CREST (No. JPMJCR18T1), JST-PRESTO (No.
JPMJPR18L7), Grant-in-Aid for Scientific Research (JSPS KAKENHI Grant Numbers
JP21H04435, JP23KJ0099, and JP23K25812), KEK-PF (Proposal numbers: 2021S2-001,
2024S2-001). The work at Beijing Institute of Technology was supported by the
National Key R&D Program of China (Grant Nos. 2020YFA0308800, 2022YFA1403400),
the Natural Science Foundation of China (Grant No. 92065109), and the Beijing
Natural Science Foundation (Grant No. Z210006). S. Suzuki acknowledges support
from GP-MS at Tohoku University. T.K. acknowledges support from GP-Spin at
Tohoku University and JST-SPRING (No. JPMJSP2114). Z.W. thanks the Analysis &
Testing Center at BIT for assistance in facility support.

## References

  * Wang  _et al._ [2023] Y. Wang, H. Wu, G. T. McCandless, J. Y. Chan, and M. N. Ali, Nat. Rev. Phys. 5, 635 (2023).
  * Yu and Li [2012] S.-L. Yu and J.-X. Li, Phys. Rev. B 85, 144402 (2012).
  * Kiesel and Thomale [2012] M. L. Kiesel and R. Thomale, Phys. Rev. B 86, 121105(R) (2012).
  * Wang  _et al._ [2013] W.-S. Wang, Z.-Z. Li, Y.-Y. Xiang, and Q.-H. Wang, Phys. Rev. B 87, 115135 (2013).
  * Kiesel  _et al._ [2013] M. L. Kiesel, C. Platt, and R. Thomale, Phys. Rev. Lett. 110, 126405 (2013).
  * Ortiz  _et al._ [2019] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, E. Ertekin, T. M. McQueen, and E. S. Toberer, Phys. Rev. Mater. 3, 094407 (2019).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, Phys. Rev. Lett. 125, 247002 (2020).
  * Ortiz  _et al._ [2021] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. L. Teicher, R. Seshadri, and S. D. Wilson, Phys. Rev. Mater. 5, 034801 (2021).
  * Yin  _et al._ [2021] Q. Yin, Z. Tu, C. Gong, Y. Fu, S. Yan, and H. Lei, Chin. Phys. Lett. 38, 037403 (2021).
  * Chen  _et al._ [2021a] H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, S. Ma, S. Ni, H. Zhang, Q. Yin, C. Gong, Z. Tu, H. Lei, H. Tan, S. Zhou, C. Shen, X. Dong, B. Yan, Z. Wang, and H.-J. Gao, Nature 599, 222 (2021a).
  * Zhao  _et al._ [2021] H. Zhao, H. Li, B. R. Ortiz, S. M. L. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, and I. Zeljkovic, Nature 599, 216 (2021).
  * Jiang  _et al._ [2021] Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, J. Ruff, L. Kautzsch, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Wang, R. Thomale, T. Neupert, S. D. Wilson, and M. Z. Hasan, Nat. Mater. 20, 1353 (2021).
  * Shumiya  _et al._ [2021] N. Shumiya, M. S. Hossain, J.-X. Yin, Y.-X. Jiang, B. R. Ortiz, H. Liu, Y. Shi, Q. Yin, H. Lei, S. S. Zhang, G. Chang, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Guguchia, S. D. Wilson, and M. Z. Hasan, Phys. Rev. B 104, 035131 (2021).
  * Wang  _et al._ [2021] Z. Wang, Y.-X. Jiang, J.-X. Yin, Y. Li, G.-Y. Wang, H.-L. Huang, S. Shao, J. Liu, P. Zhu, N. Shumiya, M. S. Hossain, H. Liu, Y. Shi, J. Duan, X. Li, G. Chang, P. Dai, Z. Ye, G. Xu, Y. Wang, H. Zheng, J. Jia, M. Z. Hasan, and Y. Yao, Phys. Rev. B 104, 075148 (2021).
  * Nie  _et al._ [2022] L. Nie, K. Sun, W. Ma, D. Song, L. Zheng, Z. Liang, P. Wu, F. Yu, J. Li, M. Shan, D. Zhao, S. Li, B. Kang, Z. Wu, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Nature 604, 59 (2022).
  * Mielke III  _et al._ [2022] C. Mielke III, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Nature 602, 245 (2022).
  * Xu  _et al._ [2022] Y. Xu, Z. Ni, Y. Liu, B. R. Ortiz, Q. Deng, S. D. Wilson, B. Yan, L. Balents, and L. Wu, Nat. Phys. 18, 1470 (2022).
  * Jiang  _et al._ [2023] Z. Jiang, H. Ma, W. Xia, Z. Liu, Q. Xiao, Z. Liu, Y. Yang, J. Ding, Z. Huang, J. Liu, Y. Qiao, J. Liu, Y. Peng, S. Cho, Y. Guo, J. Liu, and D. Shen, Nano Lett. 23, 5625 (2023).
  * Du  _et al._ [2021] F. Du, S. Luo, B. R. Ortiz, Y. Chen, W. Duan, D. Zhang, X. Lu, S. D. Wilson, Y. Song, and H. Yuan, Phys. Rev. B 103, L220504 (2021).
  * Chen  _et al._ [2021b] K. Y. Chen, N. N. Wang, Q. W. Yin, Y. H. Gu, K. Jiang, Z. J. Tu, C. S. Gong, Y. Uwatoko, J. P. Sun, H. C. Lei, J. P. Hu, and J.-G. Cheng, Phys. Rev. Lett. 126, 247001 (2021b).
  * Zhang  _et al._ [2021] Z. Zhang, Z. Chen, Y. Zhou, Y. Yuan, S. Wang, J. Wang, H. Yang, C. An, L. Zhang, X. Zhu, Y. Zhou, X. Chen, J. Zhou, and Z. Yang, Phys. Rev. B 103, 224513 (2021).
  * Du  _et al._ [2022] F. Du, S. Luo, R. Li, B. R. Ortiz, Y. Chen, S. D. Wilson, Y. Song, and H. Yuan, Chin. Phys. B 31, 017404 (2022).
  * Song  _et al._ [2021] Y. Song, T. Ying, X. Chen, X. Han, X. Wu, A. P. Schnyder, Y. Huang, J. g. Guo, and X. Chen, Phys. Rev. Lett. 127, 237001 (2021).
  * Yang  _et al._ [2022] H. Yang, Z. Huang, Y. Zhang, Z. Zhao, J. Shi, H. Luo, L. Zhao, G. Qian, H. Tan, B. Hu, K. Zhu, Z. Lu, H. Zhang, J. Sun, J. Cheng, C. Shen, X. Lin, B. Yan, X. Zhou, Z. Wang, S. J. Pennycook, H. Chen, X. Dong, W. Zhou, and H.-J. Gao, Sci. Bull. 67, 2176 (2022).
  * Liu  _et al._ [2023] Y. Liu, Y. Wang, Y. Cai, Z. Hao, X.-M. Ma, L. Wang, C. Liu, J. Chen, L. Zhou, J. Wang, S. Wang, H. He, Y. Liu, S. Cui, B. Huang, J. Wang, C. Chen, and J.-W. Mei, Phys. Rev. Mater. 7, 064801 (2023).
  * Li  _et al._ [2022a] Y. Li, Q. Li, X. Fan, J. Liu, Q. Feng, M. Liu, C. Wang, J.-X. Yin, J. Duan, X. Li, Z. Wang, H.-H. Wen, and Y. Yao, Phys. Rev. B 105, L180507 (2022a).
  * Kato  _et al._ [2022a] T. Kato, Y. Li, K. Nakayama, Z. Wang, S. Souma, F. Matsui, M. Kitamura, K. Horiba, H. Kumigashira, T. Takahashi, Y. Yao, and T. Sato, Phys. Rev. Lett. 129, 206402 (2022a).
  * Zhou  _et al._ [2023] X. Zhou, Y. Li, Z. Liu, J. Hao, Y. Dai, Z. Wang, Y. Yao, and H.-H. Wen, Phys. Rev. B 107, 125124 (2023).
  * Oey  _et al._ [2022a] Y. M. Oey, B. R. Ortiz, F. Kaboudvand, J. Frassineti, E. Garcia, R. Cong, S. Sanna, V. F. Mitrović, R. Seshadri, and S. D. Wilson, Phys. Rev. Mater. 6, L041801 (2022a).
  * Oey  _et al._ [2022b] Y. M. Oey, F. Kaboudvand, B. R. Ortiz, R. Seshadri, and S. D. Wilson, Phys. Rev. Mater. 6, 074802 (2022b).
  * Ding  _et al._ [2022] G. Ding, H. Wo, Y. Gu, Y. Gu, and J. Zhao, Phys. Rev. B 106, 235151 (2022).
  * [32] S. Peng, Y. Han, Y. Li, J. Shen, Y. Miao, Y. Luo, L. Huai, Z. Ou, H. Li, Z. Xiang, Z. Liu, D. Shen, M. Hashimoto, D. Lu, Y. Yao, Z. Qiao, Z. Wang, and J. He, arXiv:2406.17769.
  * Zhong  _et al._ [2023a] Y. Zhong, S. Li, H. Liu, Y. Dong, K. Aido, Y. Arai, H. Li, W. Zhang, Y. Shi, Z. Wang, S. Shin, H. N. Lee, H. Miao, T. Kondo, and K. Okazaki, Nat. Commun. 14, 1945 (2023a).
  * Zhong  _et al._ [2023b] Y. Zhong, J. Liu, X. Wu, Z. Guguchia, J.-X. Yin, A. Mine, Y. Li, S. Najafzadeh, D. Das, C. Mielke III, R. Khasanov, H. Luetkens, T. Suzuki, K. Liu, X. Han, T. Kondo, J. Hu, S. Shin, Z. Wang, X. Shi, Y. Yao, and K. Okazaki, Nature 617, 488 (2023b).
  * Luo  _et al._ [2023] Y. Luo, Y. Han, J. Liu, H. Chen, Z. Huang, L. Huai, H. Li, B. Wang, J. Shen, S. Ding, Z. Li, S. Peng, Z. Wei, Y. Miao, X. Sun, Z. Ou, Z. Xiang, M. Hashimoto, D. Lu, Y. Yao, H. Yang, X. Chen, H.-J. Gao, Z. Qiao, Z. Wang, and J. He, Nat. Commun. 14, 3819 (2023).
  * Huai  _et al._ [2024] L. Huai, H. Li, Y. Han, Y. Luo, S. Peng, Z. Wei, J. Shen, B. Wang, Y. Miao, X. Sun, Z. Ou, B. Liu, X. Yu, Z. Xiang, M.-Q. Kuang, Z. Qiao, X. Chen, and J. He, npj Quantum Materials 9, 23 (2024).
  * Liu  _et al._ [2022] M. Liu, T. Han, X. Hu, Y. Tu, Z. Zhang, M. Long, X. Hou, Q. Mu, and L. Shan, Phys. Rev. B 106, L140501 (2022).
  * Liu  _et al._ [2024] J. Liu, Q. Li, Y. Li, X. Fan, J. Li, P. Zhu, H. Deng, J.-X. Yin, H. Yang, J. Li, H.-H. Wen, and Z. Wang, Sci. Rep. 14, 9580 (2024).
  * Qian  _et al._ [2021] T. Qian, M. H. Christensen, C. Hu, A. Saha, B. M. Andersen, R. M. Fernandes, T. Birol, and N. Ni, Phys. Rev. B 104, 144506 (2021).
  * Zhu  _et al._ [2024] K. J. Zhu, L. P. Nie, B. Lei, M. J. Wang, K. L. Sun, Y. Z. Deng, M. L. Tian, T. Wu, and X. H. Chen, Phys. Rev. Res. 6, 023295 (2024).
  * Nakayama  _et al._ [2022] K. Nakayama, Y. Li, T. Kato, M. Liu, Z. Wang, T. Takahashi, Y. Yao, and T. Sato, Phys. Rev. X 12, 011001 (2022).
  * Peng  _et al._ [2024] C. Peng, Y. Li, X. Chen, S. Dai, Z. Wu, C. Wu, Q. Wan, K. Zhao, R. Li, S. Mo, D. Qin, S. Yu, H. Zhong, S. Yuan, J. Guo, and N. Xu, Phys. Rev. B 109, 115415 (2024).
  * Kitamura  _et al._ [2022] M. Kitamura, S. Souma, A. Honma, D. Wakabayashi, H. Tanaka, A. Toyoshima, K. Amemiya, T. Kawakami, K. Sugawara, K. Nakayama, K. Yoshimatsu, H. Kumigashira, T. Sato, and K. Horiba, Rev. Sci. Instrum. 93, 033906 (2022).
  * Nakayama  _et al._ [2021] K. Nakayama, Y. Li, T. Kato, M. Liu, Z. Wang, T. Takahashi, Y. Yao, and T. Sato, Phys. Rev. B 104, L161112 (2021).
  * Liang  _et al._ [2021] Z. Liang, X. Hou, F. Zhang, W. Ma, P. Wu, Z. Zhang, F. Yu, J.-J. Ying, K. Jiang, L. Shan, Z. Wang, and X.-H. Chen, Phys. Rev. X 11, 031026 (2021).
  * Yu  _et al._ [2022] J. Yu, Z. Xu, K. Xiao, Y. Yuan, Q. Yin, Z. Hu, C. Gong, Y. Guo, Z. Tu, P. Tang, H. Lei, Q.-K. Xue, and W. Li, Nano Lett. 22, 918 (2022).
  * Kato  _et al._ [2022b] T. Kato, Y. Li, K. Nakayama, Z. Wang, S. Souma, M. Kitamura, K. Horiba, H. Kumigashira, T. Takahashi, and T. Sato, Phys. Rev. B 106, L121112 (2022b).
  * Luo  _et al._ [2022a] Y. Luo, S. Peng, S. M. L. Teicher, L. Huai, Y. Hu, Y. Han, B. R. Ortiz, Z. Liang, Z. Wei, J. Shen, Z. Ou, B. Wang, Y. Miao, M. Guo, M. Hashimoto, D. Lu, Z. Qiao, Z. Wang, S. D. Wilson, X. Chen, and J. He, Phys. Rev. B 105, L241111 (2022a).
  * Kato  _et al._ [2023] T. Kato, Y. Li, M. Liu, K. Nakayama, Z. Wang, S. Souma, M. Kitamura, K. Horiba, H. Kumigashira, T. Takahashi, Y. Yao, and T. Sato, Phys. Rev. B 107, 245143 (2023).
  * Tan  _et al._ [2021] H. Tan, Y. Liu, Z. Wang, and B. Yan, Phys. Rev. Lett. 127, 046401 (2021).
  * Wang  _et al._ [2022] C. Wang, S. Liu, H. Jeon, and J.-H. Cho, Phys. Rev. B 105, 045135 (2022).
  * Fu  _et al._ [2021] Y. Fu, N. Zhao, Z. Chen, Q. Yin, Z. Tu, C. Gong, C. Xi, T. Park, M. Ye, and L. Balents, Phys. Rev. B 104, 035142 (2021).
  * Kato  _et al._ [2022c] T. Kato, Y. Li, T. Kawakami, M. Liu, K. Nakayama, Z. Wang, A. Moriya, K. Tanaka, T. Takahashi, Y. Yao, and T. Sato, Commun. Mater. 3, 30 (2022c).
  * Hu  _et al._ [2022a] Y. Hu, X. Wu, B. R. Ortiz, S. Ju, X. Han, J. Ma, N. C. Plumb, M. Radovic, R. Thomale, S. D. Wilson, A. P. Schnyder, and M. Shi, Nat. Commun. 13, 2220 (2022a).
  * Kang  _et al._ [2022] M. Kang, S. Fang, J.-K. Kim, B. R. Ortiz, S. H. Ryu, J. Kim, J. Yoo, G. Sangiovanni, D. D. Sante, B.-G. Park, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, S. D. Wilson, J.-H. Park, and R. Comin, Nat. Phys. 18, 301 (2022).
  * Li  _et al._ [2021] H. Li, T. T. Zhang, T. Yilmaz, Y. Y. Pai, C. E. Marvinney, A. Said, Q. W. Yin, C. S. Gong, Z. J. Tu, E. Vescovo, C. S. Nelson, R. G. Moore, S. Murakami, H. C. Lei, H. N. Lee, B. J. Lawrie, and H. Miao, Phys. Rev. X 11, 031050 (2021).
  * Lou  _et al._ [2022] R. Lou, A. Fedorov, Q. Yin, A. Kuibarov, Z. Tu, C. Gong, E. F. Schwier, B. Buchner, H. Lei, and S. Borisenko, Phys. Rev. Lett. 128, 036402 (2022).
  * Liu  _et al._ [2021] Z. Liu, N. Zhao, Q. Yin, C. Gong, Z. Tu, M. Li, W. Song, Z. Liu, D. Shen, Y. Huang, K. Liu, H. Lei, and S. Wang, Phys. Rev. X 11, 041010 (2021).
  * Cho  _et al._ [2021] S. Cho, H. Ma, W. Xia, Y. Yang, Z. Liu, Z. Huang, Z. Jiang, X. Lu, J. Liu, Z. Liu, J. Li, J. Wang, Y. Liu, J. Jia, Y. Guo, J. Liu, and D. Shen, Phys. Rev. Lett. 127, 236401 (2021).
  * Luo  _et al._ [2022b] H. Luo, Q. Gao, H. Liu, Y. Gu, D. Wu, C. Yi, J. Jia, S. Wu, X. Luo, Y. Xu, L. Zhao, Q. Wang, H. Mao, G. Liu, Z. Zhu, Y. Shi, K. Jiang, J. Hu, Z. Xu, and X. J. Zhou, Nat. Commun. 13, 273 (2022b).
  * Cai  _et al._ [2024] Y. Cai, Y. Wang, Z. Hao, Y. Liu, X. Sui, Z. Liang, X.-M. Ma, F. Zhang, Z. Shen, C. Zhang, Z. Jiang, Y. Yang, W. Liu, Q. Jiang, Z. Liu, M. Ye, D. Shen, H. Gao, H. Xiao, Z. Liu, Z. Sun, Y. Liu, S. Cui, J. Chen, L. Wang, C. Liu, J. Lin, B. Huang, Z. Wang, X. Chen, J.-W. Mei, J. Wang, and C. Chen, Commun. Mater. 5, 31 (2024).
  * Kang  _et al._ [2023] M. Kang, S. Fang, J. Yoo, B. R. Ortiz, Y. M. Oey, J. Choi, S. H. Ryu, J. Kim, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, J. G. Checkelsky, S. D. Wilson, J.-H. Park, and R. Comin, Nat. Mater. 22, 186 (2023).
  * Hu  _et al._ [2022b] Y. Hu, X. Wu, B. R. Ortiz, X. Han, N. C. Plumb, S. D. Wilson, A. P. Schnyder, and M. Shi, Phys. Rev. B 106, L241106 (2022b).
  * Li  _et al._ [2022b] C. Li, X. Wu, H. Liu, C. Polley, Q. Guo, Y. Wang, X. Han, M. Dendzik, M. H. Berntsen, B. Thiagarajan, Y. Shi, A. P. Schnyder, and O. Tjernberg, Phys. Rev. Res. 4, 033072 (2022b).
  * Momma and Izumi [2011] K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).