We can rewrite the subproblem in (158) as follows:

$\displaystyle\begin{aligned} &\min\limits_{z,y}\Big{\\{}-\langle
b,\,y\rangle+\delta_{\cal
K}^{*}(-z)+\frac{\sigma}{2}\|z+A^{*}y-\widehat{c}^{k}\|^{2}-\langle\delta^{k},\,y\rangle\Big{\\}}\\\
=&\min\limits_{y}\left\\{-\langle
b,\,y\rangle-\langle\delta^{k},\,y\rangle+\min\limits_{z}\left\\{\delta^{*}_{{\cal
K}}(-z)+\frac{\sigma}{2}\|z+A^{*}y-\hat{c}^{k}\|^{2}\right\\}\right\\}\\\
=&\min\limits_{y}\left\\{-\langle
b,\,y\rangle-\langle\delta^{k},\,y\rangle+\sigma M_{\delta^{*}_{{\cal
K}}/\sigma}(A^{*}y-\hat{c}^{k})\right\\},\end{aligned}\quad$

where the last equality comes from the definition of Moreau-Yosida envelope,
and the equality holds for $z=-\mathop{\rm{Prox}}_{\delta^{*}_{{\cal
K}}/\sigma}(A^{*}y-\hat{c}^{k})$. By the property of Moreau-Yosida proximal
mapping, $\mathop{\rm{Prox}}_{\delta^{*}_{{\cal
K}}/\sigma}(A^{*}y-\hat{c}^{k})=(A^{*}y-\hat{c}^{k})-\mathop{\rm{Prox}}_{\sigma\delta_{{\cal
K}}}(\sigma(A^{*}y-\hat{c}^{k}))$, therefore, $z=\Pi_{{\cal
K}}(\sigma(A^{*}y-\hat{c}^{k}))-(A^{*}y-\hat{c}^{k})$, which gives (158). The
derivation of (166) is similar, we omit the details.

### Statements and Declarations

The research of Kim-Chuan Toh is supported by the Ministry of Education,
Singapore, under its Academic Research Fund Tier 3 grant call
(MOE-2019-T3-1-010). All authors certify that they have no affiliations with
or involvement in any organization or entity with any financial interest or
non-financial interest in the subject matter or materials discussed in this
manuscript.

### Data availability statement

The links of datasets from the literature that are used in this paper have
been included where they are mentioned. For the other simulated data, they are
available from the corresponding author on reasonable request.

## References

  * (1) Benders, J.F.: Partitioning procedures for solving mixed-variables programming problems. Numerische Mathematik. 4(3), 238-252 (1962)
  * (2) Berkelaar, A., Dert, C., Oldenkamp, B., Zhang, S.: A primal-dual decomposition-based interior point approach to two-stage stochastic linear programming. Operations Research. 50, 904-915 (2002)
  * (3) Birge, J.R.: Current trends in stochastic programming computation and applications. Technique Report 95-15, Department of Industrial and Operations Engineering, University of Michigan (1995)
  * (4) Birge, J.R., Louveaux, F.: Introduction to Stochastic Programming. Springer-Verlag (1997)
  * (5) Birge, J.R., Qi, L.: Computing block-angular Karmarkar projections with applications to stochastic programming. Management Science. 34, 1472-1479 (1988)
  * (6) Chatzipanagiotis, N., Dentcheva, D., Zavlanos, M.M.: An augmented Lagrangian method for distributed optimization. Mathematical Programming. 152, 405-434 (2015)
  * (7) Chen, L., Li, X.D., Sun, D.F., Toh, K.C.: On the equivalence of inexact proximal ALM and ADMM for a class of convex composite programming. Mathematical Programming. 185, 111-161 (2021)
  * (8) Chen, M., Mehrotra, S.: Self-concordance and decomposition based interior point methods for the two stage stochastic convex optimization problem. SIAM J. Optimization. 21, 1667-1687 (2011)
  * (9) Chen, L., Sun, D.F., Toh, K.C.: An efficient inexact symmetric Gauss-Seidel based majorized ADMM for high-dimensional convex composite conic programming. Math. Program. 161, 237-270 (2017)
  * (10) Dantzig, G.: Upper bounds, secondary constraints and block triangularity in linear programming. Econometrica. 23, 174-183 (1955)
  * (11) Fischetti, M., Ljubić, I., Sinnl, M.: Redesigning Benders decomposition for large-scale facility location. Management science. 63(7), 2146-2162 (2017)
  * (12) G$\stackrel{{\scriptstyle..}}{{\rm u}}$nl$\stackrel{{\scriptstyle..}}{{\rm u}}$k, O., Lee, J., Weismantel, R.: MINLP strengthening for separable convex quadratic transportation-cost UFL. Technical Report RC24213 (W0703-042), Research Division, IBM, Yorktown Heights, NY (2007)
  * (13) Kall, P., Wallace, S.W.: Stochastic Programming. John Wiley & Sons, NY (1994)
  * (14) Karmarkar, N.: A New Polynomial Time Algorithm for Linear Programming. Combinatorica. 4(4), 373–395 (1984)
  * (15) King, A.J.: Stochastic programming problems: examples from the literature. In: Ermoliev, Y., Wets, R. J-B. (eds.) Numerical Techniques for Stochastic Optimization, pp. 543–567. Springer-Verlag (1988)
  * (16) Kontogiorgis, S., Meyer, R.: A variable-penalty alternating directions method for convex optimization. Mathematical Programming. 83, 29-53 (1998)
  * (17) Lam, X. Y., Sun, D., Toh, K. C.: Semi-proximal Augmented Lagrangian-Based Decomposition Methods for Primal Block-Angular Convex Composite Quadratic Conic Programming Problems. Informs Journal on Optimization, 3(3), 254-277 (2021)
  * (18) Lau, K.K., Womersley, R.S.: Multistage quadratic stochastic programming. Journal of Computational and Applied Mathematics. 129, 105-138 (2001)
  * (19) Li, M., Wu, Z.: On the Convergence Rate of Inexact Majorized sGS ADMM with Indefinite Proximal Terms for Convex Composite Programming. Asia-Pacific Journal of Operational Research, 38(01), 2050035 (2021)
  * (20) Li, X.D., Sun, D.F., Toh, K.C.: A highly efficient semismooth Newton augmented Lagrangian method for solving Lasso problems. SIAM J. Optimization. 28, 433-458 (2018)
  * (21) Li, X.D., Sun, D.F., Toh, K.C.: A block symmetric Gauss-Seidel decomposition theorem for convex composite quadratic programming and its applications. Mathematical Programming. 175, 395-418 (2019)
  * (22) Linderoth, J., Shapiro, A., Wright, S.: The empirical behavior of sampling methods for stochastic programming. Annals of Operations Research. 142, 215-241 (2006)
  * (23) Lubin, M., Hall, J.A.J., Petra, C.G., Anitescu, M.: Parallel distributed-memory simplex for large-scale stochastic LP problems. Computational Optimization and Applications. 55, 571-596 (2013)
  * (24) Mayer, J.: Stochastic Linear Programming Algorithms: A Comparison based on a Model Management System. Gordon and Breach Science Publishers (1998)
  * (25) Mehrotra, S., Özevin, M.: Decomposition-based interior point methods for two-stage stochastic semidefinite programming. SIAM Journal on Optimization. 18, 206-222 (2007)
  * (26) Mehrotra, S., Özevin, M.: Decomposition based interior point methods for two-stage stochastic convex quadratic programs with recourse. Operations Research. 57, 964-974 (2009)
  * (27) Peters, R.J., Boskma, K., Kuper, H.A.E.: Stochastic programming in production planning: a case with non-simple recourse. Statistica Neerlandica. 31, 113-126 (1977)
  * (28) Prékopa, A.: Stochastic Programming. Kluwer Academic Publishers (1995)
  * (29) Rockafellar, R.T., Wets, R.: Scenarios and policy aggregation in optimization under uncertainty. Mathematics of Operations Research. 16, 119-147 (1991)
  * (30) Ruszczyński, A.: On convergence of an augmented Lagrangian decomposition method for sparse convex optimization. Math. Oper. Res. 20, 634-656 (1995)
  * (31) Ruszczyński, A.: Some advances in decomposition methods for stochastic linear programming. Ann. Oper. Res. 85, 153-172 (1999)
  * (32) Sun, J., Xu, H.L., Zhang, M.: A new interpretation of the progressive hedging algorithm for multistage stochastic minimization problems. Journal of Industrial $\&$ Management Optimization. 16, 1655-1662 (2020)
  * (33) Torres, J., Li, C., Apap, R., Grossmann, I.: A review on the performance of linear and mixed integer two-stage stochastic programming algorithms and software. Eprints for the optimization community (2019)
  * (34) Van Slyke, R.M., Wets, R.: L-Shaped linear programs with applications to optimal control and stochastic programming. SIAM Journal on Applied Mathematics. 17(4), 638-663 (1969)
  * (35) Wollenberg, T.: Two-Stage Stochastic Semidefinite Programming: Theory, Algorithms, and Application to AC Power Flow under Uncertainty. PhD Thesis, Universität Duisburg-Essen (2017)
  * (36) Zhang, N., Wu, J., Zhang, L.: A linearly convergent majorized ADMM with indefinite proximal terms for convex composite programming and its applications. Mathematics of Computation. 89, 1867-1894 (2020)
  * (37) Zhao, G.: Interior-point methods with decomposition for solving large-scale linear programs. J. Optimization Theory and Applications. 102, 169-192 (1999)
  * (38) Zhao, G.: A log-barrier method with Benders decomposition for solving two-stage stochastic programs. Mathematical Programming. 90, 507-536 (2001)
  * (39) Zverovich, V., Fábián, C.I., Ellison, E.F.E., Mitra, G.: A computational study of a solver system for processing two-stage stochastic LPs with enhanced Benders decomposition. Mathematical Programming Computation. 4, 211–238 (2012)


# A Non-Topological Approach to Understanding Weyl Semimetals

Antonio Levy1<EMAIL_ADDRESS>Albert F. Rigosi1 Francois Joint2 Gregory
S. Jenkins3 1National Institute of Standards and Technology, Gaithersburg,
Maryland 20899, USA  2Department of Physics, University of Maryland, College
Park, MD 20742, USA  3Laboratory for Physical Sciences, College Park, MD
20740, USA

###### Abstract

In this work, chiral anomalies and Drude enhancement in Weyl semimetals are
separately discussed from a semi-classical and quantum perspective, clarifying
the physics behind Weyl semimetals while avoiding explicit use of topological
concepts. The intent is to provide a bridge to these modern ideas for
educators, students, and scientists not in the field using the familiar
language of traditional solid-state physics at the graduate or advanced
undergraduate physics level.

##  I Introduction

Weyl fermions (defined below) have historically been of interest in answering
fundamental questions about the universe, particularly the observation of the
matter-antimatter imbalance.1 The family of elementary particles classified as
fermions, or particles of half-integer spin, are important in the Standard
Model that unifies three of the four known forces of nature. Within the model
are twenty-four families of fermions. Almost all of them are massive Dirac
fermions. Within the family of Dirac fermions lies a subset class known as
Weyl fermions, the set of fermions that are massless. Those well-versed in the
physics of the weak nuclear force will recall that those interactions are
stronger for both left-chiral matter and right-chiral antimatter than their
counterparts of opposite chirality (chirality summary shown in Fig. 1).

Weyl and 3D Dirac semimetals are the only systems in which signatures of Weyl
fermions have been observed. Although they are not fundamental particles, the
excitations in these material systems offer a unique playground to study Weyl
fermion physics like the chiral anomaly.2-6 Novel properties predicted for
Weyl semimetals (WSMs) like significantly reduced scattering, Drude
enhancement along applied magnetic fields, and long-lived spin-polarized
currents in the presence of magnetic fields could lead to myriad applications
in fields like spintronics and quantum computing.

Most introductions to WSMs and Weyl Fermions are mathematically intense,
making it difficult to develop an intuitive understanding of their physics.
There are a few works that try to explain the concepts behind Weyl fermions in
an educational context.7-9 This paper aims to provide a conceptual
introduction accessible to educators, students, and scientists who have some
understanding of traditional condensed matter physics.

##  II Background on Weyl Fermions

The concept of Weyl fermions originated from the fusion of two familiar
topics. The first is the energy-momentum relation from relativity expressed
as: $E^{2}=\left(pc\right)^{2}+\left(mc^{2}\right)^{2}$. The second is the
time-dependent Schrödinger equation from quantum mechanics expressed as:
$-\frac{\hslash^{2}\triangledown^{2}\psi}{\left(2m\right)}=i\hslash\frac{\partial\psi}{\partial
t}$.

To incorporate relativity into a quantum mechanical formula, Dirac treated
each quantity in the energy-momentum relation as an operator acting on a
wavefunction. Using $p=-i\hslash\triangledown$ and
$E=i\hslash\frac{\partial}{\partial t}$, the resulting equation is:

$\left(-\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}+\triangledown^{2}\right)\psi=\frac{m^{2}c^{2}}{\hslash^{2}}\psi$ (1)

This is known as the Klein-Gordon equation. It does not include spin and is
therefore applicable to zero-spin bosons.

Dirac transformed this equation by taking the square root of the operators,
which requires consideration of all three spatial dimensions and the time
dependence. The transformation involves 4 $\times$ 4 matrices, called the
Dirac Gamma matrices, that are directly related to the Pauli matrices
$\overrightarrow{\sigma}$ that describe half-spin fermions. The Dirac equation
can be written in matrix block form with a particle-hole (p-h) – or particle-
antiparticle – basis:

$\left[\begin{array}[]{cc}\left(\frac{\varepsilon}{c}-mc\right)&-\overrightarrow{p}\cdot\overrightarrow{\sigma}\\\
-\overrightarrow{p}\cdot\overrightarrow{\sigma}&\left(-\frac{\varepsilon}{c}+mc\right)\end{array}\right]\left[\begin{array}[]{cc}\psi_{p}\\\
\psi_{h}\end{array}\right]=0,$ (2)

where $\varepsilon$ is the energy of the particle. On the other hand, rather
than using a particle-hole basis, the Weyl equation uses left- and right-
chiral particles as the basis:

$\left[\begin{array}[]{cc}mc&\frac{\varepsilon}{c}-\overrightarrow{p}\cdot\overrightarrow{\sigma}\\\
\frac{\varepsilon}{c}+\overrightarrow{p}\cdot\overrightarrow{\sigma}&mc\end{array}\right]\left[\begin{array}[]{cc}\psi_{L}\\\
\psi_{R}\end{array}\right]=0,$ (3)

Figure 1: A left-handed chiral particle has a spin $\sigma$ that is
antiparallel to its momentum $p$. The spin and momentum are parallel for a
right-handed chiral particle.

Massive particles in the Weyl equation prevent pure chiral eigenstates from
emerging out of Equation (3). Pure chiral particles must be massless. The
eigenstates of massless Weyl particles only involve the momentum and a
parallel or antiparallel spin as shown in Figure 1.

In solid state physics, the dispersion relation is a mathematical description
of the available energies for an electron moving about in a material with a
given momentum. It is frequently plotted as energy versus momentum (also
referred to as k-space). Further details on the basics of solid state physics
can be found in standard textbooks.10 The relativistic energy-momentum
dispersion relation for massless particles is $E=pc$ . For particle-like
excitations (or quasiparticles) moving through solid materials, the velocity
can be a constant other than the speed of light while still obeying the Weyl
equation. Casting the momentum in terms of wavenumber
$\overrightarrow{p}=\hslash\overrightarrow{k}$ , the operative Weyl dispersion
relation becomes $E=\hslash kv$, where the $v$ is the Fermi velocity which
plays a role similar to $c$ in the relativistic equations. The zero-energy is
defined at zero momentum and is called the Weyl point, and its presence in
Weyl semimetals is responsible for the novel physics they are predicted to
exhibit.

##  III Chiral Anomalies

To appreciate how a chiral anomaly emerges, it will help to recall a few
concepts from statistical and solid-state physics. Within any material,
electrons and other quasiparticles only have access to a limited number of
momenta. Mapping out these allowable momenta creates what is commonly referred
to as k-space (closely related to the reciprocal lattice of a crystalline
solid). The highest-energy electrons, in some sense, define the bounds of
achievable momenta, and these bounds sketch out an object in k-space called a
Fermi surface. The Fermi surface for a Weyl particle is a sphere centered on
the zero-momentum Weyl point.

Due to quantum mechanics, the allowable energies at which electrons are
permitted to exist are discrete (and by extension, are also discrete in
k-space via Fourier transformation). These quantized energies make up the band
structure of a crystal. For the remainder of this paper, it is assumed that
the Fermi energy exists solely within the Weyl band and that other bands are
sufficiently separated in energy so as not to perturb the Weyl state. With
this last concept recalled, we are ready to describe the anomaly.

Consider the divergence theorem as applied to a charge Q, charge density
$\rho$ and electrical current density j:

$\frac{dQ}{dt}\mathrm{=}\iiint\left(\frac{\partial\rho}{\partial
t}+\overrightarrow{\triangledown}\cdot\overrightarrow{j}\right)d^{3}x$ (4)

Since charge is conserved, the integral extends over the entire volume of the
system, leaving both sides of Eq. (4) equal to zero. This expression indicates
that the total amounts of charge entering and leaving the system are equal.
Combining Eq. (4) with Eq. (3) for m = 0 under non-zero electromagnetic fields
(which is done by replacing $\overrightarrow{p}$ with
$\overrightarrow{p}-\frac{e\overrightarrow{A}}{c}$ ), the difference between
the number of right- and left-chiral particles, $n_{R}-n_{L}$ can be obtained
after considerable effort:

$\frac{d}{dt}\left(n_{R}-n_{L}\right)\mathrm{\propto}\iiint\overrightarrow{E}\cdot\overrightarrow{B}d^{3}x$
(5)

Equation (5) shows that if a population of Weyl fermions is subjected to non-
zero applied electric and magnetic fields, the chirality of the population
will change with time. This is equivalent to saying that Weyl fermions of a
single chirality will be annihilated and replaced with Weyl fermions of the
opposite chirality. This is the chiral anomaly.

The chiral anomaly in a crystal is best understood intuitively by considering
the effects of applied magnetic and electric fields on the Fermi surface with
linear dispersions illustrated in Fig. 2. The momentum direction of any
quasiparticle on a spherical Fermi surface is radially outward. In the
depicted generic WSM, there are two locations in k-space around which the band
structure obeys the Weyl equations. Each location hosts one species of chiral
particle where the spin is either radially inward opposite the momentum
(orange Fermi surface pocket on the left side of each subfigure) or purely
aligned radially outward with the momentum (green Fermi surface pocket on the
right side of each subfigure). As will be discussed later, including only a
single chiral Fermi pocket violates the conservation of energy.

So what happens to the Fermi surface if you apply a magnetic field? It will
distort, expanding or contracting by an amount proportional to the dot product
of the spin and the magnetic field (Zeeman effect). To intuitively understand
these distortions, consider the left-handed chiral (orange) Fermi surface
pocket in Figure 2(e-g) at the point where the spin of a quasiparticle is
purely antiparallel with applied magnetic field. The quasiparticle’s energy
will decrease with applied field. This decrease in energy, through the linear
energy dispersion, causes the momentum to decrease. The decrease in momentum
shifts this point of the Fermi surface toward the center of the Fermi pocket
sphere. Likewise, the opposite point, where the field and spin are parallel,
shifts to higher momentum away from the center. All the points where the spin
is perpendicular to the field do not change their energy or momentum.
Visualizing the smooth interpolation between these points results in a Fermi
surface distorted into an egg-like shape with the major axis aligned with
field as shown in Figure 2(e).

Similar effects on the Fermi sea occur from the charge of the quasiparticle
subjected to an applied electric field. In this case, all negatively charged
quasiparticles, regardless of their spin or location on the Fermi surface,
decrease their (vector) momentum along the direction of the electric field
causing a net momentum shift of the Fermi surface.

With either an applied electric or magnetic field, the changes in a singly-
chiral Fermi pocket induces a net momentum and therefore an electric current.
This current also transports a net spin current. However, when both
(oppositely) chiral spherical Fermi surface pockets are considered, the
application of a magnetic field (with no electric field) causes two equal and
oppositely-oriented distortions that negate the total current. In fact, this
can be understood as the reason that all Weyl semimetals must include pairs of
oppositely-chiral Fermi pockets; a static magnetic field should not be able to
generate a perpetual current in a material system with a finite scattering
rate like Weyl semimetals, as this would violate the conservation of energy.
Furthermore, in order not to violate conservation of energy, this current
would have to be superconducting. At high fields, there would be a strong
current that would flow consistently in one direction, which would persist
independent of the direction in which an external electric field was applied
to the material. This would mean that carriers could flow in the direction
opposite to the bias applied to the material, leading to negative power
dissipation in the material as a result of an applied external electric field
(provided that field was not large enough to fundamentally alter the band
structure).

To formulate an expression for the generated currents, the quasiparticles
removed and added at various momenta with applied fields must be properly
counted. Recall the density of states, which is the derivative of the number
of states with respect to energy in the system. The density of electronic
states (g, which we define as the number of states per unit energy at a
specific point in k-space) is energy dependent and therefore becomes dependent
on k-space location ($\theta$, defined as the polar angle between the k-vector
and the applied fields) with the application of fields.

Consider for simplicity the case with both fields pointing in the
$\mathrm{-}x$ direction (as in Fig. 2). Then, for the right-chiral location,
we track the difference between right- and left-chiral particles
$n_{R}(\theta)-n_{L}(\theta)$, which is proportional to integrated subtraction
of the density of states, expressed as
$g_{R(L)}(\varepsilon_{F})=g(\varepsilon_{F}+(-)a_{1}Bcos(\theta)+a_{2}Ecos(\theta))$
for the right and left pockets, where $a_{1}$ and $a_{2}$ are constants. The
difference in sign in front of the $a_{1}Bcos(\theta)$ is due to the opposite
spin of the two pockets, highlighting the antiparallel nature of the
chirality. Approximating for small fields, $g_{R(L)}(\varepsilon_{F})$ can be
rewritten:
$g_{R(L)}(\varepsilon_{F})=g(\varepsilon_{F})+(-)a_{1}Bcos(\theta)g^{\prime}(\varepsilon_{F}+a_{2}Ecos(\theta))$
or, equivalently, we could say:
$g_{R}(\varepsilon_{F})-g_{L}(\varepsilon_{F})=2a_{1}Bcos(\theta)g^{\prime}(\varepsilon_{F}+a_{2}Ecos(\theta)$.
The difference in the population of right- and left- chiral Weyl fermions is
given by:

$\displaystyle\int_{0}^{\pi}\left(g_{R}\left(\varepsilon_{F},\theta\right)-g_{L}~{}\left(\varepsilon_{F},\theta\right)\right)d\theta\approx$
(6) $\displaystyle
2a_{1}B\int_{0}^{\pi}\cos\left(\theta\right)g^{{}^{\prime}}\left(\varepsilon_{F}+a_{2}Ecos\left(\theta\right)\right)d\theta$

This subtraction is reduced to a derivative of the density of states for small
fields. With additional symmetry considerations, the bounds of the integral
are halved as long as we take into account $E\rightarrow-E$ for the respective
halves of the Fermi pockets:

$\displaystyle
2a_{1}B\int_{0}^{\pi}\cos\left(\theta\right)g^{{}^{\prime}}\left(\varepsilon_{F}+a_{2}Ecos\left(\theta\right)\right)d\theta\approx$
(7) $\displaystyle
2a_{1}B\int_{0}^{\frac{\pi}{2}}\left(g^{{}^{\prime}}\left(\varepsilon_{F}+a_{2}E\cos\left(\theta\right)\right)-g^{{}^{\prime}}\left(\varepsilon_{F}-a_{2}Ecos\left(\theta\right)\right)\right)$
$\displaystyle\cos\left(\theta\right)d\theta$ $\displaystyle\propto
2a_{1}a_{2}EB\int_{0}^{\frac{\pi}{2}}g^{{}^{\prime\prime}}\left(\varepsilon_{F}\right)\cos^{2}\left(\theta\right)d\theta{}$

where the approximation $g^{\prime}\left(\varepsilon_{F}\pm
a_{2}E\cos\left(\theta\right)\right)=g^{\prime}\left(\varepsilon_{F}\right)\pm
a_{2}E\cos\left(\theta\right)g^{\prime\prime}\left(\varepsilon_{F}\right)$ is
used to go from the second line of Eq. (7) to the third.

With these considerations, the origin of Eq. (5) becomes apparent. Since the
first time derivative of $\left(n_{R}-n_{L}\right)$ is a nonzero quantity
under the experimental conditions, a net chiral current, the $``$chiral
anomaly,$"$ emerges and begins pumping particles from the left-chiral Fermi
pocket to the right-chiral Fermi pocket leading to a net accumulation of
right-chiral particles. Furthermore, the total number of right-chiral
particles increases faster than the corresponding decrease seen in the left-
chiral particles. For static fields, this ever-growing chirality imbalance is
eventually counteracted by quasiparticle scattering between and within the two
chiral Fermi pockets. A chirality-imbalanced equilibrium value is established
that depends on both pumping rate and (current relaxation) scattering rate.

Figure 2: The chiral anomaly in Weyl semimetals is understood by distortions
and shifts of the Fermi surface caused by the application of electric and
magnetic fields. (a,b) With no applied field, the Fermi surface is two
spheres: one Weyl pocket is left-handed chiral (orange) with spins
perpendicular to the Fermi surface (pointing inward) and the other right-
handed chiral (green) with spins (red arrows) perpendicular to the Fermi
surface (pointing outward). (c) The $k_{x}\mathrm{-}k_{y}~{}$ planar cut
through the center of the Fermi pocket spheres results in two circles that
demarcate allowable quasiparticle momenta. (d) The energy of the
quasiparticles, with momenta taken as a $k_{x}\mathrm{-}$ line cut through the
centers of the circles, is given by the linear dispersion $E=\hslash k_{x}v$
(where velocity $v$ is constant). (g) Application of an external magnetic
field in the - $x$ direction raises or lowers the energy of the quasiparticle
depending on the direction of the spin. Changes in energy are connected to
changes in momenta through the linear energy dispersion resulting in (e,f) $k$
-space shifts and distortions of the Fermi surface. (j) The addition of an
applied electric field interacts with the (negative) charge of the
quasiparticle causing a vectoral increase of the momenta in the - $x$
direction that either increases or decreases the quasiparticle kinetic energy
depending on the initial direction of the momentum. (h,i) The combined effects
of the applied fields distort and shift the Fermi surface, generating currents
that, in the depicted case, are largest for the (green) right-handed chiral
carriers.

##  IV Semi-Classical Picture of Drude Enhancement

Another interesting phenomenon regarding WSMs is the existence of Drude
enhancement. To get comfortable with this topic, one must recall the Drude
model of AC conductivity from solid state physics:

$\sigma_{D}\left(\omega\right)=\frac{\omega_{D}^{2}}{4\pi\left(\frac{1}{\tau}-i\omega\right)}$
(8)

As $\omega\rightarrow 0$ (when $B=0$ ) the DC conductivity appears as
$\sigma_{DC}=\frac{\omega_{D}^{2}\tau}{4\pi}$ , where the term
$\omega_{D}^{2}$ is referred to as the Drude weight. The Drude weight and
conductivity of a material are properties that are generally well-described by
the equations of motion for the charge carriers.

For most metallic and semi-metallic systems in the semi-classical limit,
application of a magnetic field parallel to the driving electric field does
not affect current since there is no Lorentz force. In the case of Weyl and 3D
Dirac semimetals, such an arrangement of applied fields is expected to
increase the Drude weight and this was recently verified experimentally.11
Most explanations of this very unusual effect rely heavily on the mathematical
consequences of the Berry’s phase curvature near the Weyl points. Such
explanations, though rigorous in the semi-classical limit, do not provide an
intuitive explanation of the origin of such an effect.

Again, studying the Fermi surface of the Weyl state provides a non-topological
approach to understanding this effect. Since the Weyl pockets are represented
in k-space as asymmetrically distorted (egg-shaped) Fermi surface pockets
(Fig. 2 (h)), any applied magnetic field will have additional, distinct
effects on both populations of chiral particles, pushing more of the particles
into the portion of the Fermi surface that contributes to the net current, as
shown in Figs. 2(h) and 2(i). The Drude weight is proportional to the current
generated by the applied external electric field. Integrating over the Fermi
surface of a pair of Weyl points in the semi-classical small-field limit, with
$\sqrt{\frac{2\hslash v_{x}v_{y}}{l_{B}^{2}}}\ll\varepsilon_{F}$:

$\displaystyle\omega_{D}^{2}\propto\int
d\Omega[g_{+}(\varepsilon_{F},\overrightarrow{B})+g_{-}(\varepsilon_{F},\overrightarrow{B})]{v_{F}^{\parallel}}^{2}(\theta)$
(9) $\displaystyle\propto\int
d\Omega\left[g\left(\varepsilon_{F}\right)+\frac{1}{2}g^{{}^{\prime\prime}}\left(\varepsilon_{F}\right)\left(a_{1}B\cos\left(\theta\right)\right)^{2}\right]{v_{F}^{\parallel}}^{2}\left(\theta\right)$

where $v_{F}^{\parallel}\left(\theta\right)=v_{F}\cos\left(\theta\right)$ is
the Fermi velocity component parallel to the applied fields and
$l_{B}=\sqrt[]{\frac{\hslash c}{eB}}$ is the magnetic length. The second term
in (9) yields the $B^{2}$ Drude weight enhancement obtained using the Berry’s
phase curvature.12

##  V Quantum Mechanical Picture of the Chiral Anomaly and Drude Enhancement

For a more comprehensive picture of Drude enhancement, we turn to the quantum
limit of the Drude weight of a Weyl pocket in 3D. The energy dispersion of the
nth excited Landau level (LL) for $B\parallel x$ is given by:

$\varepsilon_{n}\left(k_{x}\right)=\sqrt[]{\frac{\left(2n+\left(1+s\right)\right)\hslash
v_{y}v_{z}}{l_{B}^{2}}+\hslash^{2}v_{x}^{2}k_{x}^{2}}$ (10)

where $s=\pm 1$ correspond to the two spin eigenstates. Note that this
approach can convert the 3D system into a series of 1D systems. Every LL
(except $n=0$ ) has an energy dispersion containing an effective mass.
Therefore, only the $n=0$ Landau level is chiral, making it responsible for
all of the nontrivial behavior associated with WSMs. Consider the chiral
anomaly. For all $n\neq 0$ LLs, the application of an electric field simply
shifts electrons within the LL. They do not change the population of the
pocket or even the LL itself.

Figure 3: The energy dispersion of the various Landau levels along the
magnetic field direction is illustrated. Left-handed particles are pumped from
one Weyl pocket to the right-handed pocket. The formation of Landau levels
makes the dispersion quasi one-dimensional, due to the discretization of
kinetic energy perpendicular to the magnetic field (i.e., in the y and z
directions).

Under an electric field parallel to the magnetic field (see Fig. 3), the
number of carriers in one of the two $n=0$ LLs will decrease while the number
in the other $n=0$ LL will increase. The dispersion in the two pockets’ $n=0$
LLs are given by $\varepsilon_{+}\left(k_{x}\right)=+\hslash v_{x}k_{x}$ and
$\varepsilon_{-}\left(k_{x}\right)=-\hslash v_{x}k_{x}$ for the left- and
right-chiral cases, respectively. Applying an electric field along the
x-direction will deplete the population of one $n=0$ LL while increasing the
population of the other $n=0$ LL by the same amount.

The transfer of charge from one pocket to the other occurs entirely through
the chiral $n=0$ LL of both pockets, since charge will be redistributed within
each excited LL under an electric field. The amount of charge transfer is
proportional to the electric field component parallel to the magnetic field
and the degeneracy of the $n=0$ LL, which is proportional to
$\frac{1}{l_{B}^{2}}$ due to Landau quantization perpendicular to the magnetic
field which, in our case, is the yz-plane. From this, the chiral anomaly will
persist at both low and high field strengths in WSMs.

The spectral weight in the chiral $n=0$ LL is independent of the Fermi energy.
From the quantum mechanical perspective, the increase in spectral weight with
field can be understood as transferring more and more carriers from slower
bands (where they are closer to the vertex) to the faster, massless $n=0$ LL.
This quantum mechanical perspective on Drude enhancement gives a more nuanced
understanding of the complexity of this effect without the need to rely on
topological descriptions.

##  VI Conclusions

It is rather uncommon to find approaches to the physics behind Weyl fermions
in the literature without involving a discussion on topology. Here, we explain
the concepts of chiral anomalies and Drude enhancement in the context of Weyl
semimetals intended for those not in the field. Non-topological approaches
offer a more intuitive conceptual framework for understanding these physical
systems. This understanding is useful when designing, performing, or analyzing
data from experiments.

## References

  * [1] 1\. H. R. Quinn, Physics Today 56, 30 (2003).
  * [2] 2\. H. Weyl, Z. Phys. 56, 330 (1929).
  * [3] 3\. S. Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C. C. Lee, and S. M. Huang, Science 349, 613-617 (2015).
  * [4] 4\. L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, M. Soljačić, Science 349, 622-624 (2015).
  * [5] 5\. A. A. Soluyanov, D. Gresch, Z. Wang, Q. Wu, M. Troyer, X. Dai, and B. A. Bernevig, Nature 527, 495-498, (2015).
  * [6] 6\. B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard et al. Phys. Rev. X 5, 031013, (2015).
  * [7] 7\. P. B. Pal, Am. J. Phys. 79, 485-498, (2011).
  * [8] 8\. I. Bowser, K. Kiers, and E. Mitchell, Am. J. Phys. 88, 769-783, (2020).
  * [9] 9\. L. Lerner, Eur. J. Phys. 17, 172-175, (1996).
  * [10] 10\. Ashcroft, N.W. and Mermin, N.D., 1976. Solid state physics int. ed. Holt-Saunders, Philadelphia, 16.
  * [11] 11\. A. L. Levy, A. B. Sushkov, F. Liu, B. Shen, N. Ni, H. D. Drew, and G. S. Jenkins, Phys. Rev. B 101, 125102, (2020).
  * [12] 12\. J. Hofmann, S. Das Sarma, Phys. Rev. B 93, 241402(R), (2016).


# Black holes, multiple propagation speeds and energy extraction

Vitor Cardoso  Niels Bohr International Academy, Niels Bohr Institute,
Blegdamsvej 17, 2100 Copenhagen, Denmark  CENTRA, Departamento de Física,
Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Avenida Rovisco
Pais 1, 1049-001 Lisboa, Portugal  Shinji Mukohyama  Center for Gravitational
Physics and Quantum Information, Yukawa Institute for Theoretical Physics,
Kyoto University, 606-8502, Kyoto, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), The University of Tokyo Institutes for
Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Naritaka Oshita  Center for Gravitational Physics and Quantum Information,
Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto,
Japan  The Hakubi Center for Advanced Research, Kyoto University, Yoshida
Ushinomiyacho, Sakyo-ku, Kyoto 606-8501, Japan  RIKEN iTHEMS, Wako, Saitama,
351-0198, Japan  Kazufumi Takahashi  Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, 606-8502, Kyoto, Japan

###### Abstract

The Standard Model of particle physics predicts the speed of light to be a
universal speed of propagation of massless carriers. However, other
possibilities exist—including Lorentz-violating theories—where different
fundamental fields travel at different speeds. Black holes are interesting
probes of such physics, as distinct fields would probe different horizons.
Here, we build an exact spacetime for two interacting scalar fields which have
different propagation speeds. One of these fields is able to probe the black
hole interior of the other, giving rise to energy extraction from the black
hole and a characteristic late-time relaxation. Our results provide further
stimulus to the search for extra degrees of freedom, black hole instability,
and extra ringdown modes in gravitational-wave events.

††preprint: YITP-24-44, IPMU24-0014, RIKEN-iTHEMS-Report-24

##  I Introduction

General Relativity is an extremely successful description of the gravitational
interaction, which passed numerous tests, ranging over different orders of
magnitude in length scale and field strength [1]. The advent of gravitational-
wave astronomy [2] has opened the way for unprecedented tests of the strong-
field regime, including new tests of Einstein’s theory and the strengthening
of the black hole (BH) paradigm [3, 4, 5, 6, 7, 8, 9].

BHs are a bizarre solution of the mathematical equations of General
Relativity, which nature seems to abide to [10, 11, 12, 6]. Their role in
fundamental physics is highlighted by two powerful results. The first concerns
BH uniqueness, the most general vacuum BH solution belongs to the Kerr family
[13, 14, 15]. The second is that BH interiors harbor the failure of the
underlying theory or setup from which the very notion of BHs arises [16, 17].

Figure 1: Schematic picture of how we could probe the doubled horizons due to
the difference in the speed of propagation of fields, $\psi$ and $\chi$. The
parameter $\epsilon$ controls the speed of the propagation of $\chi$. The
ringdown emission and energy extraction from a non-spinning BH are discussed
in Sec. IV.2 and in Sec. IV.3, respectively.

In the Standard Model of particle physics, together with General Relativity,
massless fields travel at the speed of light. Hence, spacetime horizons are
common to all interactions. However, this need not be the case. One can
consider for example Lorentz-violating theories where propagation speeds are
different and hence different interactions-carriers travel at different
speeds. In such setups, one field, say $\psi$, is able to probe the region
inside the horizon of another field $\chi$, possibly carrying important
information about BH interiors. We refer to Fig. 1.

For example, in the context of ghost condensation [18], a scalar field $\phi$
with a timelike gradient $\phi_{\mu}\coloneqq\partial_{\mu}\phi$ provides a
preferred time slicing and, as a result, spontaneously breaks Lorentz
invariance. A spherically symmetric BH solution in ghost condensation was
constructed in [19] by taking into account the effect of a higher derivative
term, which is present in the theory of ghost condensation (and which was
later dubbed scordatura [20]), to the stealth Schwarzschild solution found for
the first time also in [19]. (See also [21] for an extension to more generic
scalar-tensor theories with higher derivatives.) For any practical purposes at
astrophysical scales, the metric is well approximated by the Schwarzschild
geometry, and the gradient $\phi^{\mu}\coloneqq g^{\mu\nu}\phi_{\nu}$ is
tangent to a congruence of geodesics which is regular on the BH horizon. As we
shall see in the following section, coupling between $\phi^{\mu}$ and
derivative of other scalar fields can change the propagation speeds of the
latter fields and thus provides a concrete setup for the situation depicted in
Fig. 1.

An important effect of having different horizons for different fields is the
possibility of energy extraction: inside $\chi$’s horizon, it can have
negative energy which may, via coupling between $\chi$ and $\psi$, be created
together with positive energy of $\psi$. While the former (negative energy)
remains inside the horizon, the latter (positive energy) can be extracted to
the BH exterior [22, 23]. No concrete realization of this setup is known, nor
is there any evidence that energy extraction indeed occurs. In addition, even
if no energy extraction occurs, the probing of the horizon interior may affect
the dynamical behavior in the BH exterior, i.e., ringdown or the excitation of
quasinormal modes (QNMs). Unfortunately, nothing is known about possible
characteristic signatures. The goal of this work is to fill the gap. We
discuss and study a simple example of a theory with two interacting scalar
fields, which probe different spacetime geometries, where all of the above is
realized.

The rest of this paper is organized as follows. In Section II, we introduce an
example of a theory where two scalar fields couple to two different effective
geometries and have different propagation speeds. In Section III, we provide
details of the procedure to solve numerically the relevant equations on a
fixed background. In Section IV, we provide a number of important benchmark
tests. We also show two of our main results: i. the late-time relaxation
contains imprints of both horizons, modulating the ringdown phase, with
exciting observational consequences; ii. energy extraction occurs in our
setup, via scattering of wavepackets. Finally, we draw our conclusions in
Section V.

##  II Framework

We want to investigate a setup where different fields propagate at different
speeds, by coupling to different effective spacetime geometries. In
particular, we focus on two scalar fields $\psi$ and $\chi$ that propagate on
metrics $g_{\mu\nu}$ and $\tilde{g}_{\mu\nu}$, respectively, and that interact
with each other. For that, we consider a scalar field $\psi$ coupled to a
metric $g_{\mu\nu}$, and a scalar $\chi$ coupled to a disformally transformed
version of it, $\tilde{g}_{\mu\nu}$.

Our toy model that realizes such a setup is defined by the following action,

$\displaystyle I$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\biggl{(}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\psi\partial_{\nu}\psi-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+\frac{\epsilon}{2}(\phi^{\mu}\partial_{\mu}\chi)^{2}$
(1) $\displaystyle-$
$\displaystyle\frac{1}{2}(\nabla_{\mu}\phi^{\mu})^{2}(\epsilon_{11}\chi^{2}+2\epsilon_{12}\chi\psi+\epsilon_{22}\psi^{2})\biggr{)}\,,$

where $\phi^{\mu}=g^{\mu\nu}\partial_{\nu}\phi$ with $\phi$ being a scalar
field that satisfies $g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi=-1$, and
$\epsilon$ and $\epsilon_{11,12,22}$ are dimensionless constants. Note that
the $Z_{2}$ symmetry of $\phi$ is respected and the interaction between the
fields $\chi$ and $\psi$ is turned off at infinity if $\phi^{\mu}$ is chosen
to approach the asymptotic timelike Killing vector at infinity. This action
can be rewritten as

$\displaystyle I$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\biggl{(}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\psi\partial_{\nu}\psi-\frac{1}{2}\tilde{g}^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi$
(2) $\displaystyle-$
$\displaystyle\frac{1}{2}(\nabla_{\mu}\phi^{\mu})^{2}(\epsilon_{11}\chi^{2}+2\epsilon_{12}\chi\psi+\epsilon_{22}\psi^{2})\biggr{)}\,.$

Here, the inverse disformal metric has been defined as
$\tilde{g}^{\mu\nu}\coloneqq g^{\mu\nu}-\epsilon\phi^{\mu}\phi^{\nu}$, which
corresponds to the following disformal metric:

$\displaystyle\tilde{g}_{\mu\nu}=g_{\mu\nu}+\frac{\epsilon}{1+\epsilon}\phi_{\mu}\phi_{\nu}\,.$
(3)

As can be seen from Eq. (2), the effective metric governing the dynamics of
$\chi$ is the disformal metric $\tilde{g}_{\mu\nu}$, which differs from the
metric $g_{\mu\nu}$ governing the scalar $\psi$ when $\epsilon\neq 0$.
Therefore, in general the horizon of $\chi$ is different from that of $\psi$.

The model (2) is symmetric under the exchange of $\psi$ and $\chi$ in the
sense that the action is invariant under the following replacements:

$\displaystyle\psi\to(1+\epsilon)^{1/4}\chi\,,\quad\chi\to(1+\epsilon)^{1/4}\psi\,,$
$\displaystyle
g_{\mu\nu}\to\tilde{g}_{\mu\nu}\,,\quad\phi\to(1+\epsilon)^{-1/2}\phi\,,$
$\displaystyle\epsilon\to-\frac{\epsilon}{1+\epsilon}\,,\quad\epsilon_{ij}\to(1+\epsilon)^{-1}\epsilon_{ji}\,,$
(4)

provided that $\epsilon>-1$, where we have used the fact that
$\sqrt{-\tilde{g}}=(1+\epsilon)^{-1/2}\sqrt{-g}$ and that $\tilde{g}_{\mu\nu}$
and $\nabla_{\mu}\phi^{\mu}$ transform under (4) as

$\displaystyle\tilde{g}_{\mu\nu}\to
g_{\mu\nu}\,,\qquad\nabla_{\mu}\phi^{\mu}\to(1+\epsilon)^{1/2}\nabla_{\mu}\phi^{\mu}\,.$
(5)

Under the duality transformation (4), $\epsilon\in(-1,0)$ is mapped to
$\epsilon\in(0,\infty)$, and vice versa. Therefore, in what follows, we assume
$\epsilon>0$ without loss of generality. Note that the condition
$g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi=-1$ is preserved under the
above replacements.

For concreteness, consider a spherically symmetric, static spacetime,

$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}$ $\displaystyle=$ $\displaystyle-
Adt^{2}+\frac{dr^{2}}{B}+r^{2}d\Omega^{2}$ (6) $\displaystyle=$
$\displaystyle-d\tau^{2}+\left(1-A\right)d\rho^{2}+r^{2}d\Omega^{2}\,,$

where $d\Omega^{2}$ is the metric of the unit $2$-sphere, $A=A(r)$ and
$B=B(r)$. Here, $\tau$ and $\rho$ represent Gaussian-normal (or Lemaître)
coordinates defined by

$d\tau=dt+\sqrt{\frac{1-A}{AB}}dr\,,\qquad
d\rho=dt+\frac{dr}{\sqrt{AB(1-A)}}\,.$ (7)

Note that $\partial_{\rho}r=-\partial_{\tau}r=\sqrt{B(1-A)/A}$. The condition
$g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi=-1$ can be realized if we
choose

$\phi=\tau\,,$ (8)

which is precisely the typical scalar field profile for the stealth
Schwarzschild solution in the Gaussian-normal coordinate system [19, 24, 25].
The corresponding disformally transformed metric is

$\tilde{g}_{\mu\nu}dx^{\mu}dx^{\nu}=-\frac{d\tau^{2}}{1+\epsilon}+\left(1-A\right)d\rho^{2}+r^{2}d\Omega^{2}\,.$
(9)

Note also that we assume $\epsilon>-1$ so that the coordinate $\tau$ is
timelike with respect to not only $g_{\mu\nu}$ but also $\tilde{g}_{\mu\nu}$.

In what follows, we specialize to the only vacuum static solution of General
Relativity, i.e., the Schwarzschild background where

$\displaystyle A=B=1-\frac{r_{s}}{r}\eqqcolon f\,.$ (10)

We find in the standard Schwarzschild coordinates,

$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}$ $\displaystyle=$ $\displaystyle-
fdt^{2}+\frac{dr^{2}}{f}+r^{2}d\Omega^{2}\,,$ (11)
$\displaystyle\tilde{g}_{\mu\nu}dx^{\mu}dx^{\nu}$ $\displaystyle=$
$\displaystyle-\frac{1}{1+\epsilon}\left(1-(1+\epsilon)\frac{r_{s}}{r}\right)dt^{2}+r^{2}d\Omega^{2}$
(12) $\displaystyle+$ $\displaystyle
2\sqrt{\frac{r_{s}}{r}}f^{-1}\frac{\epsilon
dtdr}{1+\epsilon}+\left(1-\frac{r_{s}}{(1+\epsilon)r}\right)\frac{dr^{2}}{f^{2}}\,.$

It should be noted that the disformal metric $\tilde{g}_{\mu\nu}$ also
describes a Schwarzschild spacetime but with a different horizon (up to a
coordinate redefinition) [26, 27] (see also [28, 29]). More concretely, the
Schwarzschild horizon for the disformal metric is given by
$\tilde{r}_{s}=(1+\epsilon)r_{\rm s}$, with $r_{s}$ being the Schwarzschild
horizon for $g_{\mu\nu}$.*1*1*1We can use the following change of variables,
$dt\to\sqrt{1+\epsilon}\,d\tilde{t}+\frac{\epsilon\sqrt{1-f}}{f(-\epsilon+(1+\epsilon)f)}dr$
with $f=1-r_{s}/r$, to bring the line element for $\tilde{g}_{\mu\nu}$ to the
form
$\tilde{g}_{\mu\nu}dx^{\mu}dx^{\nu}=-\left(1-\frac{r_{s}(1+\epsilon)}{r}\right)d\tilde{t}^{2}+\left(1-\frac{r_{s}(1+\epsilon)}{r}\right)^{-1}dr^{2}+r^{2}d\Omega^{2}\,,$
(13) which is clearly a Schwarzschild spacetime. In our setup of $\epsilon>0$,
the horizon for $\chi$ lies outside that for $\psi$.

In other words, the simple model above is a realization of a theory where
different fields probe different geometries and different horizons. Although
this is just one example among many possibilities, it provides a concrete
theoretical framework in which one can explore the rich physics of a BH with
nested multi-horizons of the effective geometries. In the following, we will
explore some of the dynamics of the theory.

##  III Numerical scheme

We note that a minimal coupling of the scalar fields to the curvature allows
for the same vacuum solutions as General Relativity. We will study the
dynamics of $\psi$ and $\chi$ perturbatively in their amplitude, therefore
fixing the background geometry to be the Schwarzschild one. The equations of
motion for $\psi$ and $\chi$ are

$\displaystyle\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\psi)-\frac{9}{4}\frac{r_{s}}{r^{3}}(\epsilon_{12}\chi+\epsilon_{22}\psi)=0\,,$
(14)
$\displaystyle\frac{1}{\sqrt{-\tilde{g}}}\partial_{\mu}(\sqrt{-\tilde{g}}\tilde{g}^{\mu\nu}\partial_{\nu}\chi)-\frac{9}{4}\frac{r_{s}}{r^{3}}(\epsilon_{11}\chi+\epsilon_{12}\psi)=0\,.$
(15)

The simplest possible solution to the above systems would be a nontrivial,
static profile for the scalars. We have looked for regular, spherically
symmetric solutions of $\psi,\chi$ down to the inner horizon at $r=r_{s}$, and
we were unable to find any set of parameters for which this was possible. In
other words, we find no linearized hairy BHs, except the trivial $\psi=\chi=0$
everywhere (and a special solution, which requires that $\psi,\chi$ are
constants and only exist for a certain combination of coupling constants). We
thus resort to a numerical evolution of time-dependent configurations.

###  III.1 Setup

Given the spherical symmetry of the background, we expand the scalar fields
$\psi$ and $\chi$ in their multipolar components,

$\displaystyle\psi(\tau,\rho,\theta,\varphi)=\sum_{\ell m}\psi_{\ell
m}(\tau,\rho)Y_{\ell
m}(\theta,\varphi)\,,\qquad\chi(\tau,\rho,\theta,\varphi)=\sum_{\ell
m}\chi_{\ell m}(\tau,\rho)Y_{\ell m}(\theta,\varphi)\,,$ (16)

where the $Y_{\ell m}$ are spherical harmonics on the two-sphere. We will drop
the $\ell m$ subscripts onwards, with the understanding that we are always
discussing the multipolar components of the field. Note also that the
azimuthal number $m$ never plays a role, due to the spherical symmetry of the
background. We then numerically solve the two partial differential equations
(PDEs) (14) and (15), which reduce to the following equations after performing
the coordinate transformation introduced in Eq. (7):

$\displaystyle-\partial_{\tau}^{2}\psi+\frac{r}{r_{s}}\partial_{\rho}^{2}\psi-\frac{\ell(\ell+1)}{r^{2}}\psi+\frac{1}{2r^{2}}\sqrt{\frac{r}{r_{s}}}\left[3r_{s}\partial_{\tau}\psi+5r\partial_{\rho}\psi\right]-\frac{9}{4}\frac{r_{s}}{r^{3}}(\epsilon_{12}\chi+\epsilon_{22}\psi)=0\,,$
(17)
$\displaystyle-(1+\epsilon)\partial_{\tau}^{2}\chi+\frac{r}{r_{s}}\partial_{\rho}^{2}\chi-\frac{\ell(\ell+1)}{r^{2}}\chi+\frac{1}{2r^{2}}\sqrt{\frac{r}{r_{s}}}\left[3r_{s}(1+\epsilon)\partial_{\tau}\chi+5r\partial_{\rho}\chi\right]-\frac{9}{4}\frac{r_{s}}{r^{3}}(\epsilon_{11}\chi+\epsilon_{12}\psi)=0\,.$
(18)

Here, we have used the following relation to simplify the equations:

$r=\left\\{\frac{3}{2}r_{s}^{1/2}(\rho-\tau)\right\\}^{2/3}\,.$ (19)

The two fields $\psi$ and $\chi$ have two distinct horizons. We then solve the
equations while covering the interior of both horizons. To avoid numerical
instability that may arise near the singularity, we will use the following
coordinates $(T,r_{\ast})$ that nicely avoid it (see Figure 2):

Figure 2: Coordinates covering the interior of the two horizons and nicely
avoid the singularity. Each curve corresponds to a constant-$T$ hypersurface.

$\displaystyle d\tau$
$\displaystyle=dT+\frac{1}{f}\left(U\sqrt{\frac{r_{s}}{r}}+V\right)dr_{\ast}\,,$
(20) $\displaystyle d\rho$
$\displaystyle=dT+\frac{1}{f}\left(U\sqrt{\frac{r}{r_{s}}}+V\right)dr_{\ast}\,,$
(21)

where

$\displaystyle U$ $\displaystyle\coloneqq 1-r_{U}/r\,,$ (22) $\displaystyle V$
$\displaystyle\coloneqq-\sqrt{U^{2}-f}\,,$ (23) $\displaystyle r_{\ast}$
$\displaystyle\coloneqq r+r_{U}\log[(r-r_{U})/r_{s}]\,,$ (24)

and $r_{U}\,(\leq r_{s}/2)$ is a constant controlling how close constant-$T$
slices can approach the singularity.*2*2*2The relation between $(T,r_{\ast})$
and $(t,r)$ is given by  $\displaystyle dt=dT+\frac{V}{f}dr_{\ast}\,,\qquad
dr=Udr_{\ast}\,.$ Note that the coefficients of $dr_{\ast}$ in (20) and (21)
are regular at $r=r_{s}$. We then have the following relations:

$\displaystyle\partial_{\rho}$
$\displaystyle=-\frac{r_{s}/r+\sqrt{r_{s}/r}(V/U)}{f}\partial_{T}+\frac{1}{U}\sqrt{\frac{r_{s}}{r}}\partial_{r_{\ast}}\,,$
(25) $\displaystyle\partial_{\tau}$
$\displaystyle=\frac{1+\sqrt{r_{s}/r}(V/U)}{f}\partial_{T}-\frac{1}{U}\sqrt{\frac{r_{s}}{r}}\partial_{r_{\ast}}\,.$
(26)

In what follows, we will express all our results in units of $r_{s}$ and thus
set $r_{s}=1$. With this, we find the equations of motion

$\displaystyle\begin{split}&\left[\partial_{T}^{2}-f\partial_{r_{\ast}}^{2}+2V\partial_{T}\partial_{r_{\ast}}+U^{2}\frac{\ell(\ell+1)}{r^{2}}\right]\psi+\frac{1}{U^{2}}\left[fV^{\prime}-\frac{UV}{2r^{2}}\right]\left(\partial_{T}+V\partial_{r_{\ast}}\right)\psi\\\
&-\frac{2U}{r}\left[1-\frac{3}{4r}\right]\partial_{r_{\ast}}\psi+\frac{2UV}{r}\partial_{T}\psi+\frac{9}{4}\frac{U^{2}}{r^{3}}(\epsilon_{22}\psi+\epsilon_{12}\chi)=0\,,\end{split}$
(27) $\displaystyle\begin{split}&\left\\{\left(1+\epsilon
W^{2}\right)\partial_{T}^{2}-f\partial_{r_{\ast}}^{2}+2V\partial_{T}\partial_{r_{\ast}}+U^{2}\frac{\ell(\ell+1)}{r^{2}}\right\\}\chi\\\
&+\frac{1}{U^{2}}\left[fV^{\prime}-\frac{UV}{2r^{2}}\right]\left(\partial_{T}+V\partial_{r_{\ast}}\right)\chi-\frac{2U}{r}\left[1-\frac{3}{4r}\right]\partial_{r_{\ast}}\chi+\frac{2UV}{r}\partial_{T}\chi+\frac{9}{4}\frac{U^{2}}{r^{3}}(\epsilon_{11}\chi+\epsilon_{12}\psi)\\\
&+\epsilon\left[-2\sqrt{\frac{1}{r}}W\partial_{T}\partial_{r_{\ast}}-U\sqrt{\frac{1}{r}}\left(\frac{W}{U}\right)^{\prime}\partial_{T}+\frac{1}{r}\partial_{r_{\ast}}^{2}-\frac{2rU^{\prime}+U^{2}}{2r^{2}U}\partial_{r_{\ast}}\right.\\\
&\left.-\frac{3U}{2r}\sqrt{\frac{1}{r}}\left(W\partial_{T}-\sqrt{\frac{1}{r}}\partial_{r_{\ast}}\right)\right]\chi=0\,,\end{split}$
(28) $\displaystyle W$ $\displaystyle\coloneqq\frac{U+\sqrt{1/r}V}{f}\,,$ (29)

where a prime denotes differentiation with respect to $r_{\ast}$. We note that
the function $W(r)$ is regular at $r=1$. Indeed, we have

$\displaystyle W=1+\frac{r_{U}^{2}}{2(1-r_{U})}+{\cal O}(r-1)\,.$ (30)

Setting $r_{U}=1/2$, the function $W(r)$ reduces to

$W=\frac{2r+2\sqrt{r}+1}{2\sqrt{r}(\sqrt{r}+1)}\,,$ (31)

which is regular in the whole range of $r_{\ast}$. In the following, we take
$r_{U}=1/2$ and apply the formula (31) to the wave equation.

The two second-order PDEs of (27) and (28) are decomposed into four first-
order PDEs:

$\displaystyle\partial_{T}\psi$ $\displaystyle=\Pi_{\psi}\,,$ (32)
$\displaystyle\begin{split}\partial_{T}\Pi_{\psi}&=f\partial_{r_{\ast}}^{2}\psi-2V\partial_{r_{\ast}}\Pi_{\psi}-U^{2}\frac{\ell(\ell+1)}{r^{2}}\psi-\frac{1}{U^{2}}\left[fV^{\prime}-\frac{UV}{2r^{2}}\right]\left(\Pi_{\psi}+V\partial_{r_{\ast}}\psi\right)\\\
&+\frac{2U}{r}\left[1-\frac{3}{4r}\right]\partial_{r_{\ast}}\psi-\frac{2UV}{r}\Pi_{\psi}-\frac{9}{4}\frac{U^{2}}{r^{3}}(\epsilon_{22}\psi+\epsilon_{12}\chi)=0\,,\end{split}$
(33) $\displaystyle\partial_{T}\chi$ $\displaystyle=\Pi_{\chi}\,,$ (34)
$\displaystyle\begin{split}(1+\epsilon
W^{2})\partial_{T}\Pi_{\chi}&=f\partial_{r_{\ast}}^{2}\chi-2V\partial_{r_{\ast}}\Pi_{\chi}-U^{2}\frac{\ell(\ell+1)}{r^{2}}\chi-\frac{1}{U^{2}}\left[fV^{\prime}-\frac{UV}{2r^{2}}\right]\left(\Pi_{\chi}+V\partial_{r_{\ast}}\chi\right)\\\
&+\frac{2U}{r}\left[1-\frac{3}{4r}\right]\partial_{r_{\ast}}\chi-\frac{2UV}{r}\Pi_{\chi}-\frac{9}{4}\frac{U^{2}}{r^{3}}(\epsilon_{11}\chi+\epsilon_{12}\psi)-\epsilon\left[-2W\sqrt{\frac{1}{r}}\partial_{r_{\ast}}\Pi_{\chi}+W\frac{U^{\prime}}{U}\sqrt{\frac{1}{r}}\Pi_{\chi}\right.\\\
&\left.-W^{\prime}\sqrt{\frac{1}{r}}\Pi_{\chi}+\frac{1}{r}\partial_{r_{\ast}}^{2}\chi-\frac{2rU^{\prime}+U^{2}}{2r^{2}U}\partial_{r_{\ast}}\chi-\frac{3U}{2r}\sqrt{\frac{1}{r}}\left(W\Pi_{\chi}-\sqrt{\frac{1}{r}}\partial_{r_{\ast}}\chi\right)\right]=0\,.\end{split}$
(35)

We numerically solve the PDEs of (32)–(35) with the following initial
conditions:*3*3*3The initial conditions may involve partially outgoing modes
even when we set the wavepacket at a distant region. This is not problematic
in our computation of the reflectivity performed later as we read the time-
domain data inside the initial position of the wavepackets. That is, outgoing
modes do not contaminate the time-domain data we use.

$\displaystyle\psi(T=0,r_{\ast})$ $\displaystyle=\psi_{\rm
ini}(r_{*})\coloneqq
A_{\psi}\cos(\Omega_{\psi}r_{\ast}+\delta_{\psi})\exp\left[-\frac{(r_{\ast}-r^{\psi}_{\ast,s})^{2}}{\sigma_{\psi}^{2}}\right],$
(36) $\displaystyle\Pi_{\psi}(T=0,r_{\ast})$
$\displaystyle=\partial_{r_{*}}\psi_{\rm
ini}(r_{*})=-A_{\psi}\left(\Omega_{\psi}\sin(\Omega_{\psi}r_{\ast}+\delta_{\psi})+2\frac{(r_{\ast}-r^{\psi}_{\ast,s})}{\sigma_{\psi}^{2}}\cos(\Omega_{\psi}r_{\ast}+\delta_{\psi})\right)\exp\left[-\frac{(r_{\ast}-r^{\psi}_{\ast,s})^{2}}{\sigma_{\psi}^{2}}\right],$
(37) $\displaystyle\chi(T=0,r_{\ast})$ $\displaystyle=\chi_{\rm
ini}(r_{*})\coloneqq
A_{\chi}\cos(\Omega_{\chi}r_{\ast}+\delta_{\chi})\exp\left[-\frac{(r_{\ast}-r^{\chi}_{\ast,s})^{2}}{\sigma_{\chi}^{2}}\right],$
(38) $\displaystyle\Pi_{\chi}(T=0,r_{\ast})$
$\displaystyle=\partial_{r_{*}}\chi_{\rm
ini}(r_{*})=-A_{\chi}\left(\Omega_{\chi}\sin(\Omega_{\chi}r_{\ast}+\delta_{\chi})+2\frac{(r_{\ast}-r^{\chi}_{\ast,s})}{\sigma_{\chi}^{2}}\cos(\Omega_{\chi}r_{\ast}+\delta_{\chi})\right)\exp\left[-\frac{(r_{\ast}-r^{\chi}_{\ast,s})^{2}}{\sigma_{\chi}^{2}}\right],$
(39)

where $A_{\psi/\chi}$, $\Omega_{\psi/\chi}$, $\sigma_{\psi/\chi}$,
$\delta_{\psi/\chi}$, and $r^{\psi/\chi}_{\ast,s}$ are arbitrary parameters to
specify the shape, initial position, and initial velocity of the input
Gaussian wavepackets. We control numerical high-frequency unstable modes by
introducing the Kreiss-Oliger dissipation [30]. We confirmed that this does
not affect the benchmark tests we perform later. We also perform the
resolution test which is described in the Appendix.

###  III.2 Diagnostics for energy extraction

In this section, we explain how we discuss the energy extraction from a BH.
For this purpose, we obtain the time-domain waveform of $\psi$ and $\chi$ at
some fixed position $r_{\ast}=r_{\ast,o}\gg 1$. We set the center of the
initial Gaussian wavepacket at $r_{\ast}=r^{\psi/\chi}_{\ast,s}>r_{\ast,o}$,
so that we have an incoming wavepacket in the early times, and after a while,
we detect waves scattered by the double-horizon BH.

The energy flux $T^{\tau\rho}$ evaluated at $r_{\ast}=r_{\ast,o}$ is given by

$T^{\tau\rho}\sim-\left[\partial_{T}\psi\partial_{r_{\ast}}\psi+(1+\epsilon)\partial_{T}\chi\partial_{r_{\ast}}\chi\right],$
(40)

where we have approximated the expression by assuming $r_{\ast,o}\gg 1$ and
have omitted an overall factor which is irrelevant in the following
discussion. Moreover, since $\psi\simeq\psi(T\mp r_{\ast})$ and
$\chi\simeq\chi(T/\sqrt{1+\epsilon}\mp r_{\ast})$ at $r_{\ast}=r_{\ast,o}\gg
1$, the expression for $T^{\tau\rho}$ can be rewritten as

$\displaystyle
T^{\tau\rho}\sim\pm\left[(\partial_{T}\psi)^{2}+(1+\epsilon)^{3/2}(\partial_{T}\chi)^{2}\right],$
(41)

with the plus (minus) sign in front corresponding to the outgoing (incoming)
flux. At a distant region, the total energy flux is given by the sum of the
two energy fluxes for $\psi$ and $\chi$. To compute the energy flux of
incoming and scattered (outgoing) flux, we separate the energy flux into
incoming and outgoing modes in the time domain by introducing the truncation
parameters $T=T_{\rm cut1/cut2}$ as

$\displaystyle\tilde{\psi}_{\rm in}\coloneqq\int^{T_{\rm
cut1}}_{0}dTe^{i\omega T}(\partial_{T}\psi)\,,\qquad\tilde{\psi}_{\rm
out}\coloneqq\int^{T_{\rm end}}_{T_{\rm cut1}}dTe^{i\omega
T}(\partial_{T}\psi)\,,$ (42) $\displaystyle\tilde{\chi}_{\rm
in}\coloneqq\int_{0}^{T_{\rm cut2}}dTe^{i\omega
T}(\partial_{T}\chi)\,,\qquad\tilde{\chi}_{\rm out}\coloneqq\int_{T_{\rm
cut2}}^{T_{\rm end}}dTe^{i\omega T}(\partial_{T}\chi)\,,$ (43)

to obtain the spectral amplitude of each signal. Then the energy flux per
frequency is

$\displaystyle\frac{dE_{\rm in}}{d\omega}\sim\left[|\tilde{\psi}_{\rm
in}|^{2}+(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm in}|^{2}\right],$ (44)
$\displaystyle\frac{dE_{\rm out}}{d\omega}\sim\left[|\tilde{\psi}_{\rm
out}|^{2}+(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm out}|^{2}\right].$ (45)

Finally, we compute the net reflectivity ${\cal R}^{2}$ as

${\cal R}^{2}\coloneqq{\cal R}_{\psi}^{2}+{\cal
R}_{\chi}^{2}=\frac{|\tilde{\psi}_{\rm
out}|^{2}+(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm out}|^{2}}{|\tilde{\psi}_{\rm
in}|^{2}+(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm in}|^{2}}\,,$ (46)

where

${\cal R}_{\psi}^{2}\coloneqq\frac{|\tilde{\psi}_{\rm
out}|^{2}}{|\tilde{\psi}_{\rm in}|^{2}+(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm
in}|^{2}}\,,\qquad{\cal
R}_{\chi}^{2}\coloneqq\frac{(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm
out}|^{2}}{|\tilde{\psi}_{\rm in}|^{2}+(1+\epsilon)^{3/2}|\tilde{\chi}_{\rm
in}|^{2}}\,.$ (47)

When the net reflectivity exceeds unity for some range of $\omega$, we
conclude that the energy extraction from a BH occurs. We shall apply the above
methodology to our model in Sec. IV.3.

##  IV Results

###  IV.1 Benchmarks

Figure 3: Evolution of a Gaussian wavepacket of the fields $\psi$ and $\chi$,
governed by the dynamical equations (27) and (28). We set $\epsilon_{ij}=0$
and $\epsilon=0.5$. The ringdown of the field $\psi$ (black solid) is
consistent with the fundamental QNM frequency (50) (red dashed). The field
$\chi$, on the other hand, couples to an effective geometry with horizon
radius $r_{s}(1+\epsilon)$ and thus decays slower. We show in the right panel
that its decay is equivalent to that of $\psi$ with the appropriate
rescalings, $T\to T/\sqrt{1+\epsilon}$ and $\omega_{\ell mn}\to\omega_{\ell
mn}/(1+\epsilon)$. Throughout the discussion, we set units such that
$r_{s}=1$.  Figure 4: Energy spectrum for the field $\chi$, for
$\epsilon=0.1$, together with a least squares fit analysis. Red solid lines
indicate the value of $\text{Re}(\tilde{\omega}_{\rm Q})$ for
$\epsilon_{11}=\pm 2$ and $0$. We find excellent agreement with the modified
QNM $\omega_{\rm Q}$ derived in (51).

We start by a few benchmark tests on the equations of motion and our numerical
implementation. Throughout the paper, unless otherwise mentioned, we consider
the quadrupolar mode ($\ell=2$), $r_{\ast,o}=60$, and set the initial Gaussian
wavepacket with $\sigma_{\psi}=2=\sigma_{\chi}$,
$\Omega_{\psi}=1=\Omega_{\chi}$, and $\delta_{\psi}=0=\delta_{\chi}$. Also, we
set $A_{\psi}=1=A_{\chi}$ in Sections IV.1 and IV.2. To model a situation
where the fields are simultaneously excited near the horizons, we set its
position at $r^{\psi/\chi}_{\ast,s}=5$ unless otherwise noted. Throughout this
section, we turn off the coupling between the fields by setting
$\epsilon_{12}=0$.

####  IV.1.1 Ringdown

In Schwarzschild-like coordinates, with the coupling turned off
($\epsilon_{12}=0$) and with $r\psi(t,r)=e^{-i\omega t}\psi_{a}(r)$,
$r\chi(t,r)=e^{-i\omega t}\psi_{b}(r)$, we find

$\displaystyle
f^{2}\psi^{\prime\prime}_{a}+ff^{\prime}\psi^{\prime}_{a}+\left(\omega^{2}-V_{a}\right)\psi_{a}=0\,,$
(48) $\displaystyle
V_{a}=f\left(\frac{\ell(\ell+1)}{r^{2}}+\frac{r_{s}}{r^{3}}(1+9\epsilon_{22}/4)\right),$
(49)

and a similar (but more complicated) equation for $\psi_{b}$. This equation
has the same form as that of massless fields, and standard methods can be used
to calculate QNMs, for example a continued-fractions approach [31, 32]*4*4*4In
particular, we find the same continued fraction representation as Leaver [31],
with his $\epsilon\to-(1+9\epsilon_{22}/4)$.. In particular, for
$\epsilon_{22}=0$, one finds a classical result

$r_{s}\omega_{220}=0.967288-0.193518i\,,$ (50)

for the fundamental QNM frequency of a quadrupolar ($\ell=2$) mode, where a
QNM frequency, $\omega=\omega_{\ell mn}$, is labeled by angular modes
$(\ell,m)$ and the overtone number $n$.*5*5*5As mentioned earlier, the QNMs
are independent of the azimuthal number $m$, and we set $m=2$ just for
concreteness. For small $\epsilon_{22}$, we find with a continued-fraction
approach,

$r_{s}(\omega_{\rm
Q}-\omega_{220})=\frac{9\epsilon_{22}}{4}(0.05112-0.0032973i)\,,$ (51)

which agrees with the parametrized results of Ref. [33] (in their notation,
$\beta_{3}^{s}=9\epsilon_{22}/4$).

Our time-domain results are summarized in Fig. 3. We scatter a Gaussian
wavepacket in the spacetime, as described previously. After an immediate
response as a consequence of direct on-light-cone propagation, the field
configuration relaxes in a series of exponentially damped sinusoids. This
stage is known as ringdown stage and is dominated by the BH QNMs [32]. Figure
3 shows a clear ringdown waveform for both fields. Fitting our numerical data,
we recover the prediction (50) with very good accuracy when $\epsilon_{22}=0$
for the $\psi$ field.

The above result concerns the field $\psi$. It is also apparent that the field
$\chi$ decays slower. In fact, as we pointed out before in footnote *1, the
effective metric probed by the field $\chi$ is simply a rescaled version of
the Schwarzschild metric. In particular, $\chi$ is coupling to a geometry with
a Schwarzschild radius $r_{s}(1+\epsilon)$ and time coordinate rescaled by
$\sqrt{1+\epsilon}$. Taking the rescalings involved, we find that the mode
frequency $\tilde{\omega}$ of $\chi$ should be related to that of $\psi$ via

$\tilde{\omega}=\frac{\omega}{(1+\epsilon)^{3/2}}\,.$ (52)

The right panel of Fig. 3 compares the decay of $\chi$ against the prediction
(52). The agreement is excellent and further validated in Fig. 4, where we
extend the comparison to other values of $\epsilon_{11}$ and do it in the
frequency domain directly. We introduce the Lorentzian functions:

$\frac{A}{\omega-\tilde{\omega}_{\rm
Q}}+\frac{B}{\omega+\tilde{\omega}^{\ast}_{\rm Q}}\,,$ (53)

where $\tilde{\omega}_{\rm Q}\coloneqq\omega_{\rm Q}/(1+\epsilon)^{3/2}$. We
then fit the model (53) with the Fourier transform of $\chi$, denoted as
$\tilde{\chi}$ in Fig. 4, by using the least square method to determine the
fitting parameters $A$ and $B$. The two free parameters are real and are
relevant to the amplitude and phase. The agreement is excellent as the peak
(real part of QNM frequency) and broadness (relevant to the quality factor) of
$\tilde{\chi}$ fit very well with the Lorentzian function with
$\tilde{\omega}_{\rm Q}$.

####  IV.1.2 Tidal numbers

It is also straightforward to calculate the tidal response in the decoupling
limit. In the standard massless scalar in a Schwarzschild background, the
tidal Love numbers are zero. By contrast, in our setup they are nontrivial.
For a quadrupolar field for example, we find a “running” of the coupling, with
the regular solution at the horizon

$\displaystyle\psi_{a}=r^{3}\left[1+{\cal
O}\left(r^{-1}\right)\right]-\frac{1}{180}\cdot\frac{9\epsilon_{22}}{4}r^{-2}\left[\log
r+{\cal O}(1)\right]\left[1+{\cal O}\left(r^{-1}\right)\right]\,,$ (54)

a result which can also be read off from the parametrized study of Ref. [34].

####  IV.1.3 Instability

Figure 5: Evolution of a Gaussian wavepacket of the fields $\psi$ and $\chi$,
governed by the dynamical equations (27) and (28). We set $\epsilon=0.1$ and
$\epsilon_{11}=\epsilon_{22}=-7$. Both $\epsilon_{11}$ and $\epsilon_{22}$ are
mass-like parameters, and large negative values trigger a tachyonic-like
instability. The fields grow exponentially at late times with a rate
$r_{s}g_{\psi}=0.6281$ and $r_{s}g_{\chi}=0.4634$ as shown in the figures.
Those values are read by the least square method. The instability rate of
$\psi$ is in excellent agreement with $r_{s}g_{\psi}=0.6283$ and
$r_{s}g_{\chi}=0.4637$, obtained from a continued-fraction solution for the
unstable mode of the corresponding equation in the frequency domain, Eq. (56).

The results so far are perturbative in $\epsilon_{22}$ (or $\epsilon_{11}$).
For large negative values of this constant, we find an unstable mode of Eq.
(49). In fact, the effective potential $V_{a}$ becomes negative in some region
of $r$ and a sufficient (but not necessary) condition for an unstable mode of
$\psi$ to appear is that [35, 36]

$\int_{r_{s}}^{\infty}dr\frac{V_{a}}{f}<0\,,$ (55)

which amounts to the condition that $2\ell(\ell+1)+1+9\epsilon_{22}/4<0$. For
$\ell=2$ for example, it is sufficient that $\epsilon_{22}<-5.8$ for an
unstable mode to appear. A similar analysis reveals that a sufficient
condition for an unstable mode of $\chi$ to appear is given by
$\epsilon_{11}/(1+\epsilon)<-5.8$.

A continued-fraction solution yields accurate values for the unstable mode,
which corresponds to a purely imaginary value of $\omega$ with a positive
imaginary part. We define the growth rate for $\psi$ and $\chi$ as
$g_{\psi/\chi}\coloneqq-i\omega$. For $\epsilon=0.1$, $\ell=2$, and
$\epsilon_{22}=-7=\epsilon_{11}$, we find

$r_{s}g_{\psi}=0.6283\,,\qquad r_{s}g_{\chi}=0.4637\,.$ (56)

Our time-domain results are shown in Fig. 5, which shows a clear exponential
growth of both fields. The growth rate is in very good agreement with the
prediction (56).

####  IV.1.4 Greybody factors

Figure 6: The greybody factors for $\psi$ (black solid) and $\chi$ (grey
solid). We set $\epsilon_{11}=\epsilon_{22}=\epsilon_{12}=0$ and $\ell=2$. The
greybody factors obtained by the Heun function technique are shown with red
dashed lines. We set $T_{\rm cut1}=82$ and $T_{\rm cut2}=99$.

A measure of the permeability of the angular momentum barrier to incoming
waves is the transmission coefficients in a scattering experiment. Given
regular fields at the horizon $r=r_{s}$, a scattering experiment for field
$\psi$, say, consists on imposing the asymptotic behavior

$\displaystyle\tilde{\psi}=\left\\{\begin{array}[]{l}A_{\rm
out}(\omega)e^{i\omega x}+A_{\rm in}(\omega)e^{-i\omega x}\,,\quad
x\to\infty\,,\\\ A_{\rm T}e^{-i\omega x}\,,\qquad
x\to-\infty\,.\end{array}\right.$ (59)

Here, $x$ is the tortoise coordinate defined by $dr/dx=1-r_{s}/r$. A similar
expression holds for $\tilde{\chi}$, the Fourier transform of $\chi$. Each of
the amplitudes $A_{\rm out},A_{\rm in},A_{\rm T}$ is arbitrary, since the
problem is linear, but their ratio is fixed by boundary condition at the
horizon. As such, one defines the BH greybody factor $\Gamma_{\ell m}$ as

$\displaystyle\Gamma_{\ell m}\coloneqq 1-|A_{\rm out}(\omega)/A_{\rm
in}(\omega)|^{2}\,.$ (60)

This is an important quantity that characterizes classical BHs and also their
quantum spectrum, in particular for the emission rate of Hawking radiation
[37]. We computed the greybody factors for $\psi$ and $\chi$ with
$\epsilon_{11}=\epsilon_{22}=\epsilon_{12}=0$ and $\epsilon=0.5$. We compute
the energy flux of the injected Gaussian wavepacket and that of the reflected
waves by following the methodology presented in Sec. III. To obtain the
injected and reflected waveforms in the time domain, we here set
$r^{\psi/\chi}_{\ast,s}=80$ for both $\psi$ and $\chi$ and read the amplitudes
at $r_{\ast,o}=60$. Our results shown in Fig. 6 are in excellent agreement
with greybody factors computed with a different analytic technique, using the
Heun function [38].

###  IV.2 QNM excitation in ringdown

Figure 7: Evolution of a Gaussian wavepacket in a geometry with
$\epsilon_{11}=\epsilon_{22}=0$, $\epsilon_{12}=0.5$, and $\epsilon=0.5$. We
show the evolution of $\psi$ on the left panel and of $\chi$ on the right
panel. The least square fit is performed with the data in $T\geq 90$ and in
$T\geq 115$ for $\psi$ and $\chi$, respectively. We also show the QNM models
with a sum of the two fundamental modes $\omega_{220}$ and
$\tilde{\omega}_{220}$ (orange dashed), and with $\omega_{220}$ only (blue
thin solid).

In the benchmark tests performed in the previous section, we neglected the
coupling between the fields $\psi$ and $\chi$. This is a fundamental piece of
our setup, and we now discuss its impact on the dynamics of linearized fields.
It is important at this stage to reiterate that both fields propagate on a
Schwarzschild background, one with radius $r_{s}$, the other with
$r_{s}(1+\epsilon)$. It is also important to highlight that, although we focus
on two coupled scalars, ultimately we want to draw conclusions for gravity as
well. As such, observable gravitational waves can be associated either with
$g_{\mu\nu}$ or $\tilde{g}_{\mu\nu}$, and we may assume that either $\psi$ or
$\chi$ is in an observable sector, while the other is in a hidden sector, with
negligible couplings to the Standard Model of particle physics. As we now
show, there are exciting imprints of the hidden sector in the observable
sector, specifically in the relaxation stage. Our results of scattering
Gaussian wavepackets in this setup are summarized in Figs. 7–10.

Figure 7 shows the late-time relaxation of field $\psi$ (left) and of field
$\chi$ (right), for $\epsilon_{11}=\epsilon_{22}=0$, $\epsilon_{12}=0.5$, and
$\epsilon=0.5$. Notice that the fields are coupled now, via nonvanishing
$\epsilon_{12}$. After the direct signal (not shown), field $\psi$ decays in
two different stages, apparent from the figure. The first stage corresponds to
interaction with the light ring in “its own” spacetime, i.e., in geometry
$g_{\mu\nu}$ of Eq. (11), at $r=3r_{s}/2$. This explains the first
exponentially damped stage. However, the coupling to $\chi$ implies that the
field $\psi$ also has access to the exterior light ring, that of geometry
(12). In other words, the field $\psi$ is sensitive to the ringdown of field
$\chi$. The effective geometry of $\chi$ corresponds to a larger-mass BH,
hence rings at lower frequencies, as we showed already. Indeed, we show in
Fig. 7 also the result of fitting the signal with the two fundamental modes of
$\omega_{220}$ and
$\tilde{\omega}_{220}\coloneqq\omega_{220}/(1+\epsilon)^{3/2}$ in Eq. (52).
The agreement is excellent.

On the other hand, the field $\chi$ has no information on the inner light ring
(to leading order at least). In fact, for this choice of parameters, the light
ring for $\psi$ overlaps with the horizon for $\chi$. This explains why $\chi$
relaxes as a clean damped sinusoid corresponding to the fundamental
quadrupolar mode of the outer horizon (52). Such a pure ringdown is apparent
in the right panel of Fig. 7.

Figure 8: Waveform of $\chi$ from the scattering of a wavepacket in a
spacetime with $\epsilon=0.1$, $\epsilon_{12}=1$, and
$\epsilon_{11}=\epsilon_{22}=0$ (black solid line). The two light rings are
located outside the external horizon of $\chi$, and therefore both fields are
modulated by QNMs of each other. We fit with two fundamental modes for the
observable and hidden sectors, $\omega_{220}$ and $\tilde{\omega}_{220}$ [see
Eq. (52)], by using the least square fit (orange dashed). The modulation at
early times is well modeled by sum of the two fundamental modes. The right
panel is an enlarged version of the left one. Using only one fundamental mode
does not provide a good match (blue thin solid line).

For $\epsilon<0.5$ on the other hand, the light ring associated to $\psi$ lies
on the exterior of the $\chi$ horizon, and one would expect both fields to
carry imprints of the two light rings. Indeed, this feature is apparent in
Fig. 8, where we evolve both fields but now for $\epsilon=0.1$. The two
fundamental modes of $\psi$ and $\chi$ are excited in the ringdown of $\chi$
and lead to the modulation pattern as shown in Fig. 8.

Figure 9: Left Panel: Waveforms for the observable sector $\chi$, for
couplings $\epsilon_{12}=0.1,1$ (black thin and red thick lines,
respectively). We set $\epsilon_{11}=0=\epsilon_{22}$ and $\epsilon=0.1$ for
both cases. Right Panel: The mismatch between the ringdown waveform and our
QNM model consisting of four $\chi$ QNMs up to $n=3$. We perform the least
square fit with the different start times of fit, $T_{i}$.

What we have shown is that, given an “observable field” $\chi$ (i.e., one
interacting with our detectors), its relaxation properties can show imprints
of inner horizons for invisible fields $\psi$ with which it couples. To
quantify the effect, we compute the mismatch ${\cal M}$ between the ringdown
waveform and the superposition of $\chi$ QNMs. The mismatch ${\cal M}$ is
defined as

$\displaystyle{\cal M}$ $\displaystyle\coloneqq
1-\left|\frac{\braket{\chi_{\rm Q}}{\chi}}{\sqrt{\braket{\chi_{\rm
Q}}{\chi_{\rm Q}}\braket{\chi}{\chi}}}\right|,$ (61) $\displaystyle\chi_{\rm
Q}$ $\displaystyle\coloneqq\sum_{n=0}^{n_{\rm
max}}A_{n}\cos{(\tilde{\omega}_{n}T+\delta_{n})}\,,$ (62)

and we use up to the third overtone, i.e., $n_{\rm max}=3$. The inner product
between two functions $a(T)$ and $b(T)$ is given by

$\braket{a}{b}\coloneqq\int_{T_{i}}^{T_{f}}dT^{\prime}a(T^{\prime})b(T^{\prime})\,,$
(63)

with $T_{f}=200$ at which the amplitude of $\chi$ is well suppressed. The
fitting parameters $A_{n}$ and $\delta_{n}$ are obtained by using the least
square method. Results are shown in Fig. 9.

Figure 10: QNM amplitudes read by the least square fit with the three QNMs:
the hidden sector’s fundamental mode (red), $\omega_{220}$, the observable
sector’s fundamental mode and first tone (grey), $\tilde{\omega}_{220}$ and
$\tilde{\omega}_{221}$. The value of $\epsilon$ is set to $\epsilon=0.1$
(left) and to $\epsilon=0.2$ (right). The start time of fit is $T=76$ (left)
and $T=82$ (right). We set $\epsilon_{11}=\epsilon_{22}=0$ and
$\epsilon_{12}=1$.

We find that, for $\epsilon_{12}=0.1$, there is no modulation in ringdown and
the least value of the mismatch is of the order of ${\cal M}\sim{\cal
O}(10^{-5})$. On the other hand, a larger interaction with $\epsilon_{12}=1$
leads to the modulation in ringdown and the mismatch is ${\cal M}\sim{\cal
O}(10^{-2})$. It means that the hidden sector $\psi$ can be probed from the
ringdown of an observable sector and may affect the QNM measurement utilizing
the fitting analysis.

We also read the amplitude of each QNM. Concerning a possible over fitting, we
here use the three QNMs: the fundamental mode and the first overtone for
$\chi$, $\tilde{\omega}_{220}$ and $\tilde{\omega}_{221}$, and another
fundamental mode for $\psi$, $\omega_{220}$. We then find that the fundamental
mode for the hidden sector can be dominant as is shown in the left panel of
Figure 10. Indeed, when the two horizons are close to each other, say
$\epsilon=0.1$, the amplitude for the hidden sector’s fundamental mode is
larger than the QNM excitation of the observable sector $\chi$. For a larger
separation between the horizons, say $\epsilon=0.2$, the fundamental mode for
the hidden sector is less dominant as is shown in the right panel of Figure
10. The QNMs of $\psi$ and $\chi$ are more or less affected by the interaction
controlled by $\epsilon_{12}$. Nevertheless, this result shows that our fits
extract the amplitudes in a stable manner. It implies that at least for our
parameter choice, the value of QNMs are not significantly affected by the
interaction of $\epsilon_{12}=1$.

###  IV.3 Superradiance

Figure 11: Left Panel: Energy flux of $\psi$ measured at $r_{\ast}=60$. Solid
black and dashed grey indicate the incoming and outgoing energy fluxes,
respectively. The arrows indicate the three cutoff parameters of $T_{\rm
cut1/cut2}$ we used in the right panel. We set $\epsilon=0.5$ and
$\epsilon_{12}=2$, $\epsilon_{11}=\epsilon_{22}=0$. Center Panel: Energy flux
of $\chi$. All used parameters are the same as in the left panel. Right Panel:
Reflectivity ${\cal R}$ for the energy flux shown in the left and center
panels. For the different three lines, we change the cutoff times $T_{\rm
cut1}$ and $T_{\rm cut2}$ to split the flux into the ingoing and outgoing
fluxes in the time domain.

We now discuss possible energy extraction in our setup, an exciting
possibility of having one field probing the interior of another’s horizon. For
that, we compute the reflectivity ${\cal R}$ defined in Eq. (46). We first
decompose the energy flux into two sectors, ingoing and outgoing fluxes, by
separating the time-domain data at $T=T_{\rm cut1/cut2}$ (see Sec. III.2). We
first numerically detect the incoming waves, and afterwards, the outgoing flux
appears in the time-domain data. We then perform the Fourier analysis as
described in Sec. III.2.

We find superradiant amplification of waves, or in other words, energy
extraction out of a BH, as the reflectivity exceeds unity at lower
frequencies. The reflectivity we obtained is insensitive to the choice of the
value of $T_{\rm cut1/cut2}$. Our results are summarized in Figs. 11–17. Here,
we fix $r^{\psi}_{\ast,s}=80\sqrt{1+\epsilon}$ for $\psi$ and
$r^{\chi}_{\ast,s}=80$ for $\chi$ so that the two wavepackets arrive
simultaneously to the near horizon region.

Fluxes are shown in Fig. 11, where we can see that the pulse is initially
ingoing, after which the wavepackets interact with the geometry and get
scattered back to large spatial regions. A clear ringdown at late stages is
seen, whose features were already discussed in Sec. IV.2. However, the right
panel shows that reflectivity ${\cal R}$ [as defined in Eq. (46)] for this set
of parameters: it is larger than unity at small frequencies. The fields are
extracting energy away from the geometry, for reasons explained heuristically
in the introduction. It is interesting that such “superradiant” energy
extraction takes place at low frequencies only, something which is explained
by the fact that high frequency waves are hard to scatter back, and simply
fall onto the inner horizon.

Figure 12: Reflectivity ${\cal R}$ with $\epsilon_{11}=\epsilon_{22}=0$ and
$\epsilon=0.5$. We take $A_{\psi}=1=A_{\chi}$.  Figure 13: Reflectivity ${\cal
R}$ with $\epsilon_{11}=\epsilon_{22}=0$, $\epsilon_{12}=2$, and
$\epsilon=0.5$.

It is natural to expect that, since there is no energy extraction in the
decoupling limit $\epsilon_{12}=0$, energy extraction is larger at larger
couplings. This expectation is borne out in our numerical results, Fig. 12. We
can also see two strong suppressions at different frequencies in ${\cal R}$.
This is caused by the fact that the two light rings for $\psi$ and $\chi$ lead
to an abrupt decrease in the reflectivity at different frequencies as the two
light rings have different heights.

Let us now study the dependence of the superradiant amplification on the
initial data. Figure 13 shows the reflectivity ${\cal R}$ for different
combinations of the initial wave amplitudes, $A_{\psi}$ and $A_{\chi}$. The
other parameters of the wavepackets are set to the same ones used in the
ringdown analysis. We consider the following initial data:
$(A_{\psi},A_{\chi})=(1,1)$, $(A_{\psi},A_{\chi})=(1,0)$, and
$(A_{\psi},A_{\chi})=(0,1)$. We find superradiant amplification when $\psi$’s
wavepacket that can access the horizon interior is injected. On the other
hand, as is apparent in Fig. 13, we observe less significant energy extraction
when $A_{\chi}=0$ and $A_{\psi}=1$ with $\epsilon>0$. In a more realistic
setup, events that could excite significantly and simultaneously both the
observable (e.g., gravitons) and hidden sectors are BH mergers. Our setup
$A_{\psi}=1=A_{\chi}$ could be a mimicker of such scenarios, to see the
superradiant amplification that may be caused by strong gravity phenomena.

Given that the two wavepackets are simultaneously injected, interference
between $\psi$ and $\chi$ affects the superradiant amplification. Figure 14
shows the reflectivity ${\cal R}$ for different initial phases of $\psi$’s
wavepacket. We can infer that the energy extraction due to the multiple speeds
of propagation is sensitive not only to the interaction but also to the
details of the injected waves, such as interference effects.

Figure 14: Reflectivity for different phases of $\psi$’s wavepacket:
$\delta_{\psi}=0$, $\pi/2$, and $\pi$. We set $\epsilon_{11}=\epsilon_{22}=0$,
$\epsilon_{12}=3$, and $\epsilon=0.5$. Another phase $\delta_{\chi}$ is set to
zero and $A_{\psi}=1=A_{\chi}$.

Let us now “freeze” the interference and study instead the dependence on the
horizon size. To this end, we inject a single wavepacket of a field that can
access the interior of the outer horizon. That is, we consider
$(A_{\psi},A_{\chi})=(1,0)$. From our result shown in Fig. 15, one can read
that both the ${\cal R}_{\psi}$ and ${\cal R}_{\chi}$ are less than unity,
which means that both fields are important to yield superradiance. Also, at
lower frequencies, the superradiance is significant for a larger separation of
the two horizons, i.e., a larger value of $\epsilon$. This is reasonable as a
larger region to cause superradiance can efficiently accommodate lower
frequency modes relevant to superradiance. On the other hand, the frequency
range in which superradiance appears is narrower for larger values of
$\epsilon$. The height of the angular momentum barrier for $\chi$ is
suppressed and it cannot contribute to the reflectivity at higher frequencies
of $\omega\gtrsim 1/[r_{s}(1+\epsilon)]$.

Figure 15: Reflectivity ${\cal R}$ with $\epsilon_{11}=\epsilon_{22}=0$ and
$\epsilon_{12}=3$. Figure 16: Reflectivity ${\cal R}$ for different multipole
modes: $\ell=2$ (red thin solid), $\ell=3$ (black dashed), and $\ell=4$ (blue
thick solid). We set $\epsilon_{11}=0=\epsilon_{22}$, $\epsilon_{12}=3$, and
$\epsilon=0.5$.

The height of the angular momentum barrier depends on the multipole mode
$\ell$ and the mass of $\psi$ and $\chi$ as well. Indeed, we find that the
superradiant frequency is larger for a large value of $\ell$ as is shown in
Figure 16. Figure 17 shows that the superradiant frequency increases for
massive cases. For a higher multipole mode or massive modes, higher frequency
modes are scattered by high potential barriers, which increases the
superradiant frequency.

Figure 17: Reflectivity ${\cal R}$ (left) for different values of
$\epsilon_{11}$ with $\epsilon_{22}=0$ and (right) for different values of
$\epsilon_{22}$ with $\epsilon_{11}=0$. We set $\epsilon_{12}=3$ and
$\epsilon=0.5$.

##  V Conclusions

We have shown an explicit realization of a theory of two coupled scalars which
probe different effective geometries, and which propagate on distinct
spacetimes. Our results are very clear: one field can probe the BH interior of
the other via their coupling, and transport outwards important information on
the underlying theory. In particular, the relaxation of BHs in our setup
proceeds via two dominant QNMs, corresponding to the fundamental modes of the
two effective geometries. Also, our setup exhibits energy extraction out of a
non-spinning BH as expected from the following reasons.

In our system, the field $\chi$ that sees the outer horizon can indirectly
access the horizon interior via the interaction with $\psi$ that probes the
inner horizon. Given a field that can access the horizon interior where
apparent negative energy is available, it could trigger the Penrose process
even without the spin of the BH [22, 23, 39]. We indeed find superradiance in
our setup, by introducing multiple speeds of propagation in the BH background.
These are exciting results. The superradiance amplification factors are large
(i.e., without any tuning of parameters we find amplification of a few %). In
addition, new phenomena might be possible, such as instabilities, if this
system is enclosed in a cavity or if the fields are massive [23].

A classical result concerning ergoregions states that asymptotically flat,
horizonless spacetimes with ergoregions are dynamically unstable [40, 41, 23,
42, 43]. We can take inspiration from such a result. If we do away with the
horizon for $\psi$, for example by filling the interior of
$r<r_{s}(1+\alpha\epsilon)$, with $\alpha<1$ positive, then $\chi$ sees a
horizon of “size” $r_{s}(1+\epsilon)$ inside of which it carries negative
energies. However, there is no horizon for $\psi$, all it sees is a star.
Nevertheless, the field $\psi$ has access to an ergoregion (that of $\chi$,
through the coupling), and one may expect instabilities to develop, for any
sign of $\epsilon_{ij}$. A precise study of this phenomena is outside the
scope of this work.

We expect that our results would apply to any situation where two (or more)
coupled degrees of freedom have different propagation speeds. Such a situation
often happens in modified gravity models where additional field(s) are present
on top of the metric in the gravity sector. For instance, in the context of
scalar-tensor theories, black hole perturbations have been studied
extensively, where the scalar mode does not travel at the same speed as that
of the metric perturbations in general (see, e.g., Refs. [44, 45, 24, 25]).
This offers an interesting possibility that the energy extraction from a BH
and the characteristic late-time relaxation, which we have demonstrated for a
simple toy model, could be found in gravitational wave observations. We leave
this issue for future work.

###### Acknowledgements.

Vitor Cardoso is grateful for the warm hospitality at YITP. We acknowledge
support by VILLUM Foundation (grant no. VIL37766) and the DNRF Chair program
(grant no. DNRF162) by the Danish National Research Foundation. V.C. is a
Villum Investigator and a DNRF Chair. V.C. acknowledges financial support
provided under the European Union’s H2020 ERC Advanced Grant “Black holes:
gravitational engines of discovery” grant agreement no. Gravitas–101052587.
Views and opinions expressed are however those of the author only and do not
necessarily reflect those of the European Union or the European Research
Council. Neither the European Union nor the granting authority can be held
responsible for them. This project has received funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 101007855 and No 101131233. The work of
S.M. was supported in part by Japan Society for the Promotion of Science
Grants-in-Aid for Scientific Research No. 24K07017 and the World Premier
International Research Center Initiative (WPI), MEXT, Japan. N.O. is supported
by the Grant-in-Aid for Scientific Research (KAKENHI) project for FY2023
(23K13111) and by the Hakubi project at Kyoto University. K.T. was supported
by JSPS (Japan Society for the Promotion of Science) KAKENHI Grant Nos.
JP22KJ1646 and JP23K13101. *

##  Appendix A Resolution test and sanity check

When we solve the PDEs (32)–(35) numerically, we need to specify the time step
$\Delta T$ and the spatial grid size $\Delta r_{\ast}$. We here fix the ratio
between them as

$\lambda\coloneqq\frac{\Delta T}{\Delta r_{\ast}}=0.1\,,$ (64)

which is less than unity and satisfies the Courant condition. We then solve
the PDEs in the range of $-100\leq r_{\ast}\leq 600$ with $7500$ bins, i.e.,
$\Delta r_{\ast}=700/7500$. We perform the resolution test with three
different resolutions as is shown in Figure 18.

Figure 18: Resolution test with $\Delta r_{\ast}=700/6000$ (low), $700/7500$
(medium), and $700/9000$ (high).

We also check the superradiant amplification, discussed in Sec. IV.3, is
stable against the position of the initial wavepacket. To exclude the
interference effect on the superradiance reported in Figure 14, we inject
$\psi$’s wavepacket only (see Figure 19).

Figure 19: Reflectivity ${\cal R}$ for different initial positions of $\psi$’s
wavepacket: $r_{\ast,s}=80$, $125$, and $170$. We here set
$\epsilon_{11}=\epsilon_{22}=0$, $\epsilon_{12}=3$, $\epsilon=0.5$,
$A_{\psi}=1$ and $A_{\chi}=0$. We read the time-domain data at
$r_{\ast,o}=60$.

## References

  * Will [2014] C. M. Will,  _Living Rev. Rel._ 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. Lett._ 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Berti  _et al._ [2015] E. Berti _et al._ ,  _Class. Quant. Grav._ 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Yunes  _et al._ [2016] N. Yunes, K. Yagi, and F. Pretorius,  _Phys. Rev. D_ 94, 084002 (2016), arXiv:1603.08955 [gr-qc] .
  * Barack  _et al._ [2019] L. Barack _et al._ ,  _Class. Quant. Grav._ 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Cardoso and Pani [2019] V. Cardoso and P. Pani,  _Living Rev. Rel._ 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. D_ 100, 104036 (2019), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. D_ 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), arXiv:2112.06861 [gr-qc] .
  * Hawking and Israel [1987] S. W. Hawking and W. Israel, eds.,  _THREE HUNDRED YEARS OF GRAVITATION_ (Cambridge University Press, 1987).
  * Misner  _et al._ [1973] C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman, San Francisco, 1973).
  * Chandrasekhar [1985] S. Chandrasekhar,  _The mathematical theory of black holes_ (Oxford University Press, 1985).
  * Carter [1971] B. Carter,  _Phys. Rev. Lett._ 26, 331 (1971).
  * Robinson [1975] D. C. Robinson,  _Phys. Rev. Lett._ 34, 905 (1975).
  * Chrusciel  _et al._ [2012] P. T. Chrusciel, J. Lopes Costa, and M. Heusler,  _Living Rev. Rel._ 15, 7 (2012), arXiv:1205.6112 [gr-qc] .
  * Penrose [1965] R. Penrose,  _Phys. Rev. Lett._ 14, 57 (1965).
  * Penrose [1969] R. Penrose,  _Riv. Nuovo Cim._ 1, 252 (1969).
  * Arkani-Hamed  _et al._ [2004] N. Arkani-Hamed, H.-C. Cheng, M. A. Luty, and S. Mukohyama,  _JHEP_ 05, 074 (2004), arXiv:hep-th/0312099 .
  * Mukohyama [2005] S. Mukohyama,  _Phys. Rev. D_ 71, 104019 (2005), arXiv:hep-th/0502189 [hep-th] .
  * Motohashi and Mukohyama [2020] H. Motohashi and S. Mukohyama,  _JCAP_ 01, 030 (2020), arXiv:1912.00378 [gr-qc] .
  * De Felice  _et al._ [2023] A. De Felice, S. Mukohyama, and K. Takahashi,  _JCAP_ 03, 050 (2023), arXiv:2212.13031 [gr-qc] .
  * Eling  _et al._ [2007] C. Eling, B. Z. Foster, T. Jacobson, and A. C. Wall,  _Phys. Rev. D_ 75, 101502 (2007), arXiv:hep-th/0702124 .
  * Brito  _et al._ [2015] R. Brito, V. Cardoso, and P. Pani,  _Lect. Notes Phys._ 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
  * Khoury  _et al._ [2020] J. Khoury, M. Trodden, and S. S. Wong,  _JCAP_ 11, 044 (2020), arXiv:2007.01320 [astro-ph.CO] .
  * Takahashi and Motohashi [2021] K. Takahashi and H. Motohashi,  _JCAP_ 08, 013 (2021), arXiv:2106.07128 [gr-qc] .
  * Dubovsky and Sibiryakov [2006] S. L. Dubovsky and S. M. Sibiryakov,  _Phys. Lett. B_ 638, 509 (2006), arXiv:hep-th/0603158 .
  * Mukohyama [2009] S. Mukohyama,  _JHEP_ 09, 070 (2009), arXiv:0901.3595 [hep-th] .
  * Takahashi  _et al._ [2019] K. Takahashi, H. Motohashi, and M. Minamitsuji,  _Phys. Rev. D_ 100, 024041 (2019), arXiv:1904.03554 [gr-qc] .
  * Ben Achour  _et al._ [2020] J. Ben Achour, H. Liu, and S. Mukohyama,  _JCAP_ 02, 023 (2020), arXiv:1910.11017 [gr-qc] .
  * Kreiss and Oliger [1973] H. Kreiss and J. Oliger,  _Methods for the approximate solution of time dependent problems_, GARP Publications Series (SELBSTVERL.) FEBR, 1973).
  * Leaver [1985] E. W. Leaver,  _Proc. Roy. Soc. Lond. A_ 402, 285 (1985).
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets,  _Class. Quant. Grav._ 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Cardoso  _et al._ [2019] V. Cardoso, M. Kimura, A. Maselli, E. Berti, C. F. B. Macedo, and R. McManus,  _Phys. Rev. D_ 99, 104077 (2019), arXiv:1901.01265 [gr-qc] .
  * Katagiri  _et al._ [2024] T. Katagiri, T. Ikeda, and V. Cardoso,  _Phys. Rev. D_ 109, 044067 (2024), arXiv:2310.19705 [gr-qc] .
  * Buell and Shadwick [1995] W. F. Buell and B. A. Shadwick,  _American Journal of Physics_ 63, 256 (1995).
  * Cardoso  _et al._ [2009] V. Cardoso, M. Lemos, and M. Marques,  _Phys. Rev. D_ 80, 127502 (2009), arXiv:1001.0019 [gr-qc] .
  * Page [1976] D. N. Page,  _Phys. Rev. D_ 13, 198 (1976).
  * Gregory  _et al._ [2021] R. Gregory, I. G. Moss, N. Oshita, and S. Patrick,  _Class. Quant. Grav._ 38, 185005 (2021), arXiv:2103.09862 [gr-qc] .
  * Oshita  _et al._ [2021] N. Oshita, N. Afshordi, and S. Mukohyama,  _JCAP_ 05, 005 (2021), arXiv:2102.01741 [gr-qc] .
  * Friedman [1978] J. L. Friedman,  _Commun. Math. Phys._ 63, 243 (1978).
  * Sato and Maeda [1978] H. Sato and K.-i. Maeda,  _Prog. Theor. Phys._ 59, 1173 (1978).
  * Moschidis [2018] G. Moschidis,  _Commun. Math. Phys._ 358, 437 (2018), arXiv:1608.02035 [math.AP] .
  * Mukohyama [2000] S. Mukohyama,  _Phys. Rev. D_ 61, 124021 (2000), arXiv:gr-qc/9910013 .
  * Kobayashi  _et al._ [2014] T. Kobayashi, H. Motohashi, and T. Suyama,  _Phys. Rev. D_ 89, 084042 (2014), arXiv:1402.6740 [gr-qc] .
  * Franciolini  _et al._ [2019] G. Franciolini, L. Hui, R. Penco, L. Santoni, and E. Trincherini,  _JHEP_ 02, 127 (2019), arXiv:1810.07706 [hep-th] .