††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# Generalized Design Principles for Hydrodynamic Electron Transport in
Anisotropic Metals

Yaxian Wang  John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge, MA 02138, USA  Georgios Varnavides  John A.
Paulson School of Engineering and Applied Sciences, Harvard University,
Cambridge, MA 02138, USA  Department of Materials Science and Engineering,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Research
Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge,
MA 02139, USA  Polina Anikeeva  Department of Materials Science and
Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Research Laboratory of Electronics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Johannes Gooth  Max Planck Institute for Chemical
Physics of Solids, Nöthnitzer Strasse 40, 01187 Dresden, Germany  Claudia
Felser  Max Planck Institute for Chemical Physics of Solids, Nöthnitzer
Strasse 40, 01187 Dresden, Germany  Prineha Narang<EMAIL_ADDRESS>John A. Paulson School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA 02138, USA

###### Abstract

Interactions of charge carriers with lattice vibrations, or phonons, play a
critical role in unconventional electronic transport of metals and semimetals.
Recent observations of phonon-mediated collective electron flow in bulk
semimetals, termed electron hydrodynamics, present new opportunities in the
search for strong electron-electron interactions in high carrier density
materials. Here we present the general transport signatures of such a second-
order scattering mechanism, along with analytical limits at the Eliashberg
level of theory. We study electronic transport, using _ab initio_
calculations, in finite-size channels of semimetallic ZrSiS and TaAs2 with and
without topological band crossings, respectively. The order of magnitude
separation between momentum-relaxing and momentum-conserving scattering
length-scales across a wide temperature range make both of them promising
candidates for further experimental observation of electron hydrodynamics.
More generally, our calculations show that the hydrodynamic transport regime
can be realized in a much broader class of anisotropic metals and does not, to
first order, rely on the topological nature of the bands. Finally, we discuss
general design principles guiding future search for hydrodynamic candidates,
based on the analytical formulation and our _ab initio_ predictions. We find
that systems with strong electron-phonon interactions, reduced electronic
phase space, and suppressed phonon-phonon scattering at temperatures of
interest are likely to feature hydrodynamic electron transport. We predict
that layered and/or anisotropic semimetals composed of half-filled $d$-shells
and light group V/VI elements with lower crystal symmetry are promising
candidates to observe hydrodynamic phenomena in future.

##  I Introduction

Hydrodynamic electron transport, where charge carriers can flow collectively
akin to a classical fluid, has recently garnered significant attention as a
probe of strong electron interactions in conductors, with increasing
technological relevance as electronic devices approach the micro- and nano-
meter scale. For instance, in a hydrodynamic conductor resistive processes
occur predominantly at the boundaries, which alters the spatial distribution
of Joule heating and can thereby significantly impact thermal design. Further,
it has been demonstrated that in a narrow conducting channel, collective flow
can transfer charge more efficiently than the ballistic regime, thus achieving
“superballistic” transport [1]. Microscopically, this requires that the total
momentum of electrons is conserved, with the momentum-relaxing interactions of
electrons with impurities, lattice vibrations, or the device boundary being
significantly slower. The hydrodynamic transport regime has until recently
been inaccessible since at high temperatures the electron momentum is often
relaxed by lattice vibrations, while at low temperatures by extrinsic
scattering due to high impurity concentrations. In these systems charge
transport is governed by diffusive processes, where the electrons lose their
momentum after traveling an average “mean free path” distance ($l_{\rm mr}$).

To observe hydrodynamic effects, the momentum conserving length scale ($l_{\rm
mc}$) needs to dominate, necessitating electron-electron scattering to be
frequent enough. Since the strength of conventional Coulomb interactions
depends inversely on the carrier density [2], hydrodynamic transport is most
likely to occur in semiconductors or low carrier-density semimetals such as
graphene. This is supported by the first observations of hydrodynamic electron
flow in two dimensional electron gases (2DEG) (Al,Ga)As [3], and more recently
graphene [4, 5, 6]. However, potential applications and optimization of this
transport regime are hindered because in semiconductors sufficiently high
carrier concentrations are usually achieved by impurities, which in turn
dominate the scattering at low temperatures; while in bulk metals with high
carrier densities conventional Coulomb interactions are screened.

The growth of high quality single crystals with very low levels of impurities
has facilitated recent observations of hydrodynamic transport in a handful of
high carrier-density systems such as semimetals WP2 [7] and WTe2 [8], and
delafossite metals PdCoO2 and PtCoO2 [9, 10], where the observation could not
readily be explained using the Coulomb interaction [8, 11]. Instead, a
combination of theory and experiment shows an indirect electron-electron
interaction, mediated by a virtual phonon, dominates at intermediate
temperatures [8, 11]. This provides an opportunity to realize and optimize
hydrodynamic effects in a broader family of materials, ideally with high
carrier mobility and tunable sample quality. Despite these theoretical
advances in understanding the macroscopic observables of hydrodynamic
phenomena [12, 13, 14], general design principles guiding the discovery of
bulk hydrodynamic candidates remain elusive. In this _Article_ , we explore
the phonon-mediated electron-electron scattering mechanism in anisotropic
metals and semimetals with _ab initio_ calculations. We investigate the
possibility of hydrodynamic electron flow in semimetals ZrSiS and TaAs2, the
former of which has three-dimensional Dirac nodal lines while the later has no
Dirac or Weyl crossings. We find that both materials host strong phonon-
mediated electron-electron interactions and are promising candidates to
exhibit hydrodynamic behavior at relatively higher cryogenic temperatures than
previously reported systems. Inspired by the recent development of
temperature-dependent imaging techniques that can spatially resolve the
electron current profile via nitrogen-vacancy magnetometry [15, 16, 8], we
compute the current density profiles for various combinations of momentum
relaxing and momentum conserving length scales, and characterize the resulting
transport regimes at different temperatures to provide guidance for
experimentally-relevant hydrodynamic observables.

Further, in this _Article_ we investigate the analytical limits of the phonon-
mediated scattering lifetime, and provide empirical evidence to support these
limits in light of this work and previously explored semimetals [17, 18, 8].
At temperatures significantly lower than the system’s Debye temperature, when
anharmonic phonon-phonon interactions are significantly slower than electron-
phonon interactions, resulting in long-lived phonon-excitations which
preferentially transfer their momenta to the electronic system, the second-
order phonon-mediated electron lifetime is shorter than the first-order
electron-phonon lifetime. This effectively shifts the focus in the search for
hydrodynamic candidates to ultra-pure materials with long momentum-relaxing
mean free paths, significantly expanding the pool for future search. Finally,
we discuss the essential features leading to hydrodynamic electron transport
in metallic candidates: (a) strong electron-phonon coupling, typically found
in low-symmetry crystals composed of $d/p$ orbitals [19]; (b) suppressed
phonon-phonon interactions, e.g. realized through the “acoustic bunching”
effect in systems with relatively large atomic mass difference [20]; (c) high
Fermi velocity often correlated with highly dispersive electronic bands as
well as low levels of disorder. From these predictions and observations we
conclude that anisotropic quantum materials composed of transition metals and
group V/VI elements host multiple materials families suitable for exhibiting
and optimizing hydrodynamic electron flow.

##  II Overview of material properties in phonon-mediated hydrodynamic
candidates

At first glance, the observation of hydrodynamic electron transport in three
dimensional bulk conductors appears serendipitous in nature, since
conventional electron-electron Coulomb interactions are expected to be short-
ranged due to electron screening effects. More specifically, recent
temperature- and spatially-resolved measurements on WTe2 show evidence of non-
uniform flow at $\sim 20$ K [8]. Meanwhile, bulk transport signatures on WP2
feature strong deviations from the Wiedemann-Franz law [7], while the
temperature-dependence of phonon linewidths in WP2 cannot be readily explained
via lattice anharmonicity models and provides evidence of strong phonon-
electron coupling [18]. Similarly, the width-dependent resistivity in narrow
conducting channels of PdCoO2 indicated a $l_{\rm mc}/l_{\rm mr}\sim 0.1$
ratio meaning momentum-conserving mean free paths are as short as a few
microns [10], in contrast to conventional electron-electron interactions from
the Fermi-liquid theory estimating $l_{\rm mc}$ on the order of thousands of
microns [21]. In a series of recent theoretical predictions, phonon-mediated
electron-electron scattering is shown to be important for the non-diffusive
electron transport in these systems [17, 7, 18, 8, 11]. Since momentum-
relaxing scattering needs to be minimized to observe hydrodynamic effects, it
is important for material candidates to exhibit long momentum-relaxing
$l_{\mathrm{mr}}$. In this section we examine the metals and semimetals WP2,
WTe2, and PdCoO2 that exhibit signatures of hydrodynamics and the material
properties that enable this physics.

Figure 1: (a) Schematic energy dispersion diagram along a high symmetry
direction for a semimetal, with highly-dispersive $s/p$ electron-like (blue)
and “heavy” $d/f$ hole-like (red) bands crossing the Fermi energy. Due to the
separation of the Fermi pockets, lattice excitations with a large wave numbers
are needed for transition between $k_{s}$ and $k_{d}$. (b) Fermi surface
diagram of the semimetal shown in (a) across the first Brillouin zone, showing
intra-pocket (top-left) and inter-pocket (bottom) scattering events. The
square lattice Brillouin zone, with reciprocal lattice vector $G$, is shown in
solid lines. Adapted from Ref. 22. (c) Electron-phonon momentum relaxing
scattering event, where initial electronic state $\boldsymbol{k}$ is scattered
into final state $\boldsymbol{k^{\prime}}$, exchanging both energy and
momentum to the phonon mode
$\boldsymbol{q}=\boldsymbol{k^{\prime}}-\boldsymbol{k}$. The thermally
accessible states with energy $\sim k_{B}T$ around the Fermi energy are
indicated by the shaded area. (d) Electron-electron momentum conserving
scattering event, where two initial electronic states $\boldsymbol{k_{1}}$ and
$\boldsymbol{k_{2}}$ are scattered into final states
$\boldsymbol{k_{1}^{\prime}}$ and $\boldsymbol{k_{2}^{\prime}}$. The
scattering process conserves both energy and momentum. (e) A faceted Fermi
surface, on which a scattering event between wave vectors from
$\boldsymbol{k}$ to $\boldsymbol{k^{\prime}}$ may lead to negligible change in
the electron velocity such that the total momentum stays “quasi” conserved.
Figure 2: (a) Crystal lattice of ZrSiS, highlighting the layered structure.
(b-c) Electron-phonon lifetimes ($\tau_{\rm eph}$) at (b) 298 K and (c) 4 K
projected on the Fermi surface of ZrSiS. Similarly, while both electron and
hole pockets have long lived carriers, lowering temperature quenches the
scattering events on the hole pockets much more significantly.

Intuitively, delocalized electrons, for example from $s$ orbitals, usually
lead to large dispersion in band energy, that is, larger Fermi velocity,
$v_{F}$. However, they typically- result in symmetric bonding states often
found in isotropic lattices making them less sensitive to lattice vibrations.
To access the second-order phonon-mediated interactions, we seek strong
electron-phonon coupling across the Fermi surface (see Section IV), therefore
more localized electrons, for example hybridized $d$ orbitals, are
advantageous. By contrast, suppression of the first-order momentum-relaxing
scattering can be achieved in a number of ways [22, 23] (each individually
effective though a confluence would be beneficial): I) In semimetals with well
separated electron/hole Fermi pockets across the Brillouin zone. Crystals with
lower symmetry and complex orbital hybridization usually have their conduction
band minimum and valance band maximum shifted from the zone center, so their
Fermi surfaces are located away from each other. Here, at low temperatures,
phonon excitations lack the momentum to couple electrons from one pocket to
another, and only intra-band scattering is allowed (fig. 1a-b); II) At low
temperatures where the phonon phase-space is significantly reduced, momentum-
relaxing scattering is ineffective at relaxing the system’s momentum
appreciably (fig. 1c), while momentum-conserving electron-electron scattering
via the ‘instantaneous’ emission and re-absorption of phonon modes (fig. 1d)
is allowed to take place for larger $\boldsymbol{q}$; III) On Fermi surfaces
with non-trivial geometry where small scattering efficiency factor leads to
“quasi” conservation of electron momentum (fig. 1e). Specifically, at very low
temperatures, small phonon momenta $\boldsymbol{q}$ always result in a small
angle between the incident, $\boldsymbol{k}$, and scattered electron momentum,
$\boldsymbol{k^{\prime}}=\boldsymbol{k}-\boldsymbol{q}$. On a Fermi surface
with non-trivial geometry, the electron velocity $v_{F}$ can deviate
significantly from the wave vector and thus affect the transport-relevant
scattering efficiency [11, 24]. We can quantify this effect by taking into
account the scattering angle after scattering [17, 24]:

$1-\mathrm{cos}\theta=1-\frac{v_{n\boldsymbol{k}}\cdot
v_{m\boldsymbol{k+q}}}{|v_{n\boldsymbol{k}}||v_{m\boldsymbol{k+q}}|}.$ (1)

One such example is shown in fig. 1e, when the concave portions of the Fermi
surface, result in a strong velocity and mean free path distribution
anisotropy [25, 26], and can further reduce the momentum relaxing scattering
efficiency.

Recent first principles calculations have provided significant insights into
the state-resolved momentum relaxing lifetimes in these systems [17, 8, 11].
WP2 [17] and WTe2 [8] both show a similar trend in the carrier lifetime
distribution on their Fermi surfaces. At high temperatures, the momentum
relaxing lifetimes are distributed more evenly on electron/hole pockets
compared to extremely low temperatures, where the long lived electrons are
“focused” to specific spots on the hole pockets, which have open Fermi surface
shape and thus reduced electronic phase space. For PdCoO2, which has a faceted
open Fermi surface in a hexagonal shape, at high temperatures, long-lived
electrons span the entire Brillouin zone along both $k_{z}$ and $k_{x/y}$
[11]; while at 4 K, the long-lived electrons are located only at the corners
of the hexagonal Fermi surface at $k_{z}=0$. Taken together, the highly
dispersive hole bands give rise to large Fermi velocity and smaller phase
space for electron-phonon scattering. Yet, the strong electron-phonon coupling
arising from $d$-orbitals, leads to a strong second-order electron-electron
interaction mediated by a virtual phonon. These emerging common features
motivate us to search for hydrodynamic candidates in layered and/or
anisotropic structures with transition metals and group V/VI elements, and to
postulate generalized design principles. Meanwhile, various electronic
transport signatures are indicative of the, otherwise hard to experimentally
obtain, electron-phonon coupling strength. For example, large RRR (residual-
resistivity ratio) is an important sign of low-impurity samples, and large MR
(magnetoresistance) can indicate charge compensation, unconventional Fermi
surface topology and scattering mechanisms.

##  III Ab-initio Predictions of Electron Hydrodynamics in anisotropic metals

Following the design principles discussed in Sec. II, we investigate layered
compounds with transition metals that show unconventional transport
signatures, more specifically MR. This is because while hydrodynamic electron
flow has only been experimentally reported in a few bulk systems, anomalous MR
has been heavily studied in recent literature. Large positive (and mostly
linear) MR has been found in many semimetals including but not limited to WTe2
[27, 28], LaBi [29], NbIrTe4 [30] etc. In addition, large extremely
anisotropic MR is reported in Dirac nodal line semimetal ZrSiS [31, 32] and
ZrSi(Se,Te) [33]. Further, negative MR was found in topological semimetals
[34, 35, 36], a semimetal TaAs2 without Dirac dispersion [37], and at the
LaAlO3/SrTiO3 interface [38]. These observations suggest the pool of
candidates for hydrodynamic electron flow in bulk materials may be much larger
than previously thought. In this section we present calculations of ZrSiS and
TaAs2 as case studies to evaluate these design principles.

Figure 3: (a) Temperature dependence of momentum relaxing electron-phonon
($l_{\rm mr}$) and momentum conserving electron-electron ($l_{\rm mc}$) mean
free paths in ZrSiS calculated from _ab initio_. The impurity mean free path
is taken as 1 $\mu{\rm m}$ from Uykur _et al._ [39]. (b) Normalized $j_{x}$
curvature phase diagram. Overlaid lines show the trajectories with the
decreasing temperature with different channel widths.

ZrSiS crystallizes in the PbFCl-type tetragonal $P4/nmm$ structure (fig. 2a),
and has been shown to host a three-dimensional Dirac line node and feature a
chiral anomaly and extremely large non-saturating MR [40, 41, 42, 31, 43, 32].
Non-saturating MR results from charge compensation of its high mobility
carriers with a ratio of $\sim 0.94$ electrons to holes, similar to WTe2.
Moreover, the tube-shaped Fermi surface oriented along the crystallographic
$\hat{c}$ direction leads to open electron orbits under magnetic fields, which
can give rise to strong anisotropic MR [32, 43]. These signatures indicate the
possibility hydrodynamic electron flow in this material that we investigate
next.

Figure 2b-c shows the Fermi surface of ZrSiS, composed of four degenerate
electron and hole pockets due to the mirror symmetry in the $ab$ plane all in
open shape, in agreement with previous work [40, 32, 41, 31]. The behavior of
momentum-relaxing electron-phonon lifetimes ($\tau_{\rm eph}$) at 298 K and 4
K shows substantial similarity with those in other systems mentioned above,
such as WTe2 and WP2. While at 298 K longer-lived carriers appear evenly
distributed between electron and hole pockets, at 4 K the only long-lived
carriers are located on hole pockets at the zone boundary. Using the formalism
developed in Sec. IV, the phonon-mediated electron-electron lifetime
($\tau_{\rm ee(ph)}$) is much shorter than the momentum relaxing electron-
phonon lifetimes across all temperatures (fig. 3a). Therefore, the crossover
between the dominant scattering mechanism is determined by the impurity
scattering mean free path ($l_{\mathrm{imp}}$). High quality ZrSiS samples
have been shown to have an impurity mean free path $l_{\mathrm{imp}}>1\mu{\rm
m}$ [39, 44]. Using $l_{\mathrm{imp}}=1\mu{\rm m}$, momentum relaxing events
appear to dominate below $\sim 20$ K (fig. 3a), allowing for the possibility
of observing hydrodynamic flow at higher temperatures. Further, the expected
current density ($j_{x}$) profiles with various combinations of $l_{\rm mc}$
and $l_{\rm mr}$ in a narrow conducting channel are computed from solving the
spatially-resolved Boltzmann transport equation (BTE), discussed in our prior
work [8, 11]. The contour plot of the $j_{x}$ curvature gives a general metric
to characterize different transport regime limits. Using the $l_{\rm mc}$ and
$l_{\rm mr}$ values obtained from _ab initio_ calculations, we can bridge the
microscopic scattering mechanisms with these observables and examine how
strong the hydrodynamic effect is in samples with different dimensions (fig.
3b). Our results suggest hydrodynamic electron flow can be realized in narrow,
sub-micron devices. Alternatively, if sample quality is further improved to
$l_{\mathrm{imp}}\approx 5\mu{\rm m}$, the hydrodynamic window can be expanded
to lower temperatures and wider devices. Compared to that of WTe2,
$l_{\mathrm{mc}}$ is relatively small in ZrSiS, due to a larger carrier
concentration ($\sim 9\times 10^{19}{\rm cm}^{-3}$ in ZrSiS [31] versus $\sim
2\times 10^{19}{\rm cm}^{-3}$ in WTe2 [8]).

While hydrodynamic signatures in the topological semimetals we considered so
far (ZrSiS, WP2, and WTe2) are similar, it is intriguing to address the
necessity of topologically protected states in the context of hydrodynamic
transport. The role of topology in Dirac-Weyl semimetals that exhibit
signatures of hydrodynamics has sparked intense debates in the field. To
address this open question, we turn to TaAs2 as a case study for a semimetal
without nontrivial band crossings. TaAs2 crystallizes in the monoclinic
$C12/m1$ structure with two chemical sites for As atoms, one forming Ta-As
planes with Ta atoms and the other bridging the interlayer coupling (fig. 4a).
As a non magnetic material, TaAs2 has been shown to have giant MR from high
mobility compensated charge carriers [45, 37, 46, 47, 48]. However, the band
crossings are gapped in the presence of spin orbit coupling, making it an
interesting case to study hydrodynamic behavior in topologically trivial
semimetals. Owing to a small effective mass $m^{*}\sim 0.3m_{0}$ [45, 47], the
carrier mean free path is estimated to be $\sim 10\mu{\rm m}$ at 2 K,
significantly longer than that of WTe2 and ZrSiS. Moreover, in a high quality
TaAs2 crystal with carrier density $\approx 2.8\times 10^{18}$ cm-3, the MR is
nearly one magnitude larger than that of WTe2 [46], and the carrier mobility
is up to $\approx 1.2\times 10^{5}$ cm2/V/s [37]. More importantly, the
quadratic field dependence of its MR obeys the semiclassical model,
representative of conventional (semi)metals.

Figure 4: (a) Crystal lattice of TaAs2, highlighting the monoclinic structure.
(b-c) Electron-phonon lifetimes ($\tau_{\rm eph}$) at 298 K and (d-e) 4 K
projected on the Fermi surface of TaAs2. Different views are presented due to
the low symmetry Brillouin zone and the complexity of the electron and hole
pockets. The features are highlighted that hole pockets at low temperatures
feature much longer momentum relaxing lifetimes.

These observations motivate us to examine in further detail the temperature
dependent mean free paths for both momentum relaxing and momentum conserving
events, shown in fig. 5a. Similar to ZrSiS, the phonon mediated electron-
electron mean free path is shorter than the first-order electron-phonon mean
free path for all temperature ranges studied. Combined with a long impurity
mean free path, this provides a wide temperature window, as well as a wide
range of channel widths in which to expect hydrodynamic flow (fig. 5b).
Compared to ZrSiS (Fig. 3b) and WTe2 [8], hydrodynamic flow in TaAs2 appears
more pronounced, calling for experimental verification.

Investigating the momentum relaxing lifetime distribution on the Fermi surface
of TaAs2 (fig. 4b-e), we find similar features to those discussed before.
TaAs2 has the valence band maximum located at the L point while the conduction
band minimum is shifted between $\Gamma$ and L, leading to the coexistence of
a few pairs of electron/hole Fermi surfaces located at the boundary of the
Brillouin zone. The electron band forms a closed Fermi surface while the hole
band forms an open Fermi surface [45]. At high temperatures, the distribution
of $\tau_{\rm eph}$ lifetimes is narrow and fairly evenly distributed between
the electron and hole pockets. At low temperatures, however, the electron
lifetimes are considerably shorter than those of hole carriers, consistent
with other semimetals discussed in this work. Despite attempts to classify
TaAs2 as a new kind of topological material with $\mathbb{Z}_{2}$ invariant
(0;111) from density functional theory calculations [37], the absence of
linear Dirac or Weyl band dispersion indicates that topologically protected
band crossings are not a salient ingredient for observing hydrodynamic
electron flow. The possibility of observing hydrodynamic flow in TaAs2 is
intriguing, suggesting that strong phonon-mediated electron-electron
interactions can be found in topologically trivial metals.

##  IV Electron-Phonon Coupling Interactions and the analytical limits

Our ab initio calculations suggest that, in materials with strong electron-
phonon interactions, $l_{\mathrm{mc}}<l_{\mathrm{mr}}$ across all temperatures
of interest. In this section, we seek to understand this observation by
comparing the second-order phonon-mediated electron-electron interaction
analytically to the first-order electron-phonon interaction. To first order in
the atomic displacements, the coupled electron-phonon system is described by
the Hamiltonian [49]:

$\displaystyle\hat{H}=\sum_{n\boldsymbol{k}}\epsilon_{n\boldsymbol{k}}\hat{c}^{\dagger}_{n\boldsymbol{k}}\hat{c}_{n\boldsymbol{k}}+\sum_{\boldsymbol{q}\nu}\hbar\omega_{\boldsymbol{q}\nu}\left(\hat{a}^{\dagger}_{\boldsymbol{q}\nu}\hat{a}_{\boldsymbol{q}\nu}+1/2\right)$
$\displaystyle\;+N^{-1/2}\sum_{\boldsymbol{k}\boldsymbol{q}mn\nu}g_{mn\nu}(\boldsymbol{k},\boldsymbol{q})\hat{c}_{m\boldsymbol{k+q}}^{\dagger}\hat{c}_{n\boldsymbol{k}}\left(\hat{a}_{\boldsymbol{q}\nu}+\hat{a}^{\dagger}_{-\boldsymbol{q}\nu}\right),$
(2)

where $\epsilon_{n\boldsymbol{k}}$ is the single-particle energy for an
electron with crystal momentum $\boldsymbol{k}$ in band $n$,
$\omega_{\boldsymbol{q}\nu}$ is the frequency for a phonon with crystal
momentum $\boldsymbol{q}$ and polarization $\nu$,
$\hat{c}^{\dagger}_{n\boldsymbol{k}}$ and $\hat{c}_{n\boldsymbol{k}}$
($\hat{a}^{\dagger}_{\boldsymbol{q}\nu}$ and $\hat{a}_{\boldsymbol{q}\nu}$)
are the fermionic (bosonic) creation and annihilation operators respectively,
$g_{mn\nu}(\boldsymbol{k},\boldsymbol{q})$ is the matrix element coupling
electrons and phonons, and $N$ is the number of unit cells in the phonon
supercell.

Using Fermi’s golden rule to first and second order respectively, we can
arrive at expressions for the inverse lifetimes of the electronic system due
to the coupling with phonons [17]:

$\displaystyle\frac{1}{\tau_{n\boldsymbol{k}}^{\mathrm{eph}}}$
$\displaystyle=\frac{2\pi}{\hbar}\sum_{m\nu\pm}\int\frac{d\boldsymbol{q}}{\Omega_{\mathrm{BZ}}}\delta\left(\epsilon_{m\boldsymbol{k+q}}-\epsilon_{n\boldsymbol{k}}\mp\hbar\omega_{\boldsymbol{q}\nu}\right)$
$\displaystyle\times\left[n_{\boldsymbol{q}\nu}+\frac{1}{2}\mp\left(\frac{1}{2}-f_{m\boldsymbol{k+q}}\right)\right]\left|g_{mn\nu}(\boldsymbol{k},\boldsymbol{q})\right|^{2}$
(3) $\displaystyle\frac{1}{\tau_{n\boldsymbol{k}}^{\mathrm{ee(ph)}}}$
$\displaystyle=\frac{2\pi}{\hbar}\sum_{mpr}\int\frac{d\boldsymbol{q}}{\Omega_{\mathrm{BZ}}}\int\frac{d\boldsymbol{k^{\prime}}}{\Omega_{\mathrm{BZ}}}\left|M_{nmpr}(\boldsymbol{k},\boldsymbol{k^{\prime}},\boldsymbol{q})\right|^{2}$
$\displaystyle\times\left[f_{p\boldsymbol{k^{\prime}}}f_{m\boldsymbol{k+q}}+f_{r\boldsymbol{k^{\prime}+q}}\left(1-f_{p\boldsymbol{k^{\prime}}}-f_{m\boldsymbol{k+q}}\right)\right]$
$\displaystyle\times\delta\left(\epsilon_{n\boldsymbol{k}}+\epsilon_{r\boldsymbol{k^{\prime}+q}}-\epsilon_{p\boldsymbol{k^{\prime}}}-\epsilon_{m\boldsymbol{k+q}}\right),$
(4)

where $f$ and $n$ are the Fermi-Dirac and Bose-Einstein equilibrium
distribution functions of electrons and phonons respectively, and
$\Omega_{\mathrm{BZ}}$ is the volume of the first Brillouin zone. Equations 3
and 4 describe the momentum-relaxing electron-phonon
$\left(\tau_{n\boldsymbol{k}}^{\mathrm{e-ph}}\right)$ and momentum-conserving
phonon-mediated electron-electron
$\left(\tau_{n\boldsymbol{k}}^{\mathrm{ee(ph)}}\right)$ interactions at first
and second order respectively, with

$\displaystyle
M_{nmpr}(\boldsymbol{k},\boldsymbol{k^{\prime}},\boldsymbol{q})=\sum_{\nu}\frac{g_{mn\nu}(\boldsymbol{k},\boldsymbol{q})^{*}g_{pr\nu}(\boldsymbol{k^{\prime}},\boldsymbol{q})}{\hbar\omega_{\boldsymbol{q\nu}}+\epsilon_{n\boldsymbol{k}}-\epsilon_{m\boldsymbol{k+q}}+i\eta}.$
(5)

Since we typically only consider interactions near the Fermi level, we
introduce Fermi-surface averaged versions of eqs. 3 and 4 at the Eliashberg
level of theory [50, 17]:

$\displaystyle\tau_{\rm eph}^{-1}$
$\displaystyle=\frac{\pi\beta}{2g(\epsilon_{F})}\sum_{\nu}\int\frac{d\boldsymbol{q}}{\Omega_{\mathrm{BZ}}}G_{\boldsymbol{q}\nu}\frac{\omega_{\boldsymbol{q}\nu}}{\sinh^{2}(\hbar\beta\omega_{\boldsymbol{q}\nu}/2)}$
(6) $\displaystyle\tau_{\rm ee(ph)}^{-1}$ $\displaystyle=\frac{2\pi}{\hbar
g(\epsilon_{F})}\sum_{\nu}\int\frac{d\boldsymbol{q}}{\Omega_{\mathrm{BZ}}}G_{\boldsymbol{q}\nu}^{2}\gamma(\hbar\beta\bar{\omega}_{\boldsymbol{q}\nu}).$
(7)

Here, $\beta=1/k_{B}T$ is the dimensionless inverse temperature,
$g(\epsilon_{F})$ is the density of states at the Fermi level per unit cell,
and $G_{\boldsymbol{q}\nu}$ is the dimensionless Fermi-surface integrated
electron-phonon coupling strength:

$\displaystyle
G_{\boldsymbol{q}\nu}=\sum_{nm}\int\frac{g_{s}d\boldsymbol{k}}{\Omega_{\mathrm{BZ}}}\left|g_{mn\nu}(\boldsymbol{k},\boldsymbol{q})\right|^{2}\delta(\epsilon_{n\boldsymbol{k}}-\epsilon_{F})\delta(\epsilon_{m\boldsymbol{k+q}}-\epsilon_{F}),$
(8)

and $\gamma(x)$ is the complex-valued integral, evaluated at the complex
phonon frequency:

$\displaystyle\gamma(x)$
$\displaystyle\equiv\int_{-\infty}^{\infty}dy\frac{1}{4}\frac{y^{2}}{\sinh^{2}(y/2)|x-y|^{2}}$
(9) $\displaystyle\bar{\omega}_{\boldsymbol{q}\nu}$
$\displaystyle=\omega_{\boldsymbol{q}\nu}+i\tau_{\boldsymbol{q}\nu}^{-1}/2=\omega_{\boldsymbol{q}\nu}+i\pi\omega_{\boldsymbol{q}\nu}G_{\boldsymbol{q}\nu}.$

Equations 6 and 7 can now directly be compared, by seeking a suitable function
approximation to $\gamma(x)$. To this end, we express the integrand of eq. 9
as a power series near $\hbar\beta\omega_{\boldsymbol{q}\nu}$ using a (0,2)
Padé approximant:

$\displaystyle\gamma(\hbar\beta\bar{\omega}_{\boldsymbol{q}\nu})$
$\displaystyle\approx\int_{-\infty}^{\infty}dy\left(\hbar\beta\omega_{\boldsymbol{q}\nu}\right)^{2}/$
$\displaystyle\left\\{\cosh(\hbar\beta\omega_{\boldsymbol{q}\nu})\left[2y^{2}+6\left(\pi
G_{\boldsymbol{q}\nu}y\right)^{2}+\left(\pi
G_{\boldsymbol{q}\nu}y\hbar\beta\omega_{\boldsymbol{q}\nu}\right)^{2}-2y\hbar\beta\omega_{\boldsymbol{q}\nu}\left(2+8\left(G_{\boldsymbol{q}\nu}\pi\right)^{2}+\left(G_{\boldsymbol{q}\nu}\pi\hbar\beta\omega_{\boldsymbol{q}\nu}\right)^{2}\right)\right.\right.$
$\displaystyle\qquad\left.\left.+\left(\hbar\beta\omega_{\boldsymbol{q}\nu}\right)^{2}\left(2+12\left(G_{\boldsymbol{q}\nu}\pi\right)^{2}+\left(G_{\boldsymbol{q}\nu}\pi\hbar\beta\omega_{\boldsymbol{q}\nu}\right)^{2}\right)\right]-2\left(6\pi^{2}G_{\boldsymbol{q}\nu}^{2}+1\right)\hbar\beta\omega_{\boldsymbol{q}\nu}^{2}-2\left(3\pi^{2}G_{\boldsymbol{q}\nu}^{2}+1\right)y^{2}\right.$
$\displaystyle\qquad\left.+4\left(1+4\pi^{2}G_{\boldsymbol{q}\nu}^{2}\right)y\hbar\beta\omega_{\boldsymbol{q}\nu}-2\pi^{2}G_{\boldsymbol{q}\nu}^{2}\hbar\beta\omega_{\boldsymbol{q}\nu}(\hbar\beta\omega_{\boldsymbol{q}\nu}-y)(3\hbar\beta\omega_{\boldsymbol{q}\nu}-2y)\sinh(\hbar\beta\omega_{\boldsymbol{q}\nu})\right\\}$
(10)

Figure 5:  (a) Temperature dependence of momentum relaxing electron-phonon
($l_{\rm mr}$) and momentum conserving electron-electron ($l_{\rm mc}$) mean
free paths in TaAs2 calculated from _ab initio_. The impurity mean free path
is estimated as 10 $\mu{\rm m}$ from earlier works [45, 37]. (b) Normalized
$j_{x}$ curvature phase diagram. Overlaid lines show the trajectories with the
decreasing temperature with different channel width.

Applying the fundamental theorem of calculus to the antiderivative of eq. 10
and simplifying the expression we obtain:

$\displaystyle\gamma(\hbar\beta\bar{\omega}_{\boldsymbol{q}\nu})$
$\displaystyle\approx\frac{\hbar\beta\omega_{\boldsymbol{q}\nu}}{2G_{\boldsymbol{q}\nu}\sinh^{\frac{3}{2}}\left(\frac{\hbar\beta\omega_{\boldsymbol{q}\nu}}{2}\right)\sqrt{\pi^{2}G_{\boldsymbol{q}\nu}^{2}\left(\left((\hbar\beta\omega_{\boldsymbol{q}\nu})^{2}+8\right)\sinh\left(\frac{\hbar\beta\omega_{\boldsymbol{q}\nu}}{2}\right)-4\hbar\beta\omega_{\boldsymbol{q}\nu}\cosh\left(\frac{\hbar\beta\omega_{\boldsymbol{q}\nu}}{2}\right)\right)+4\sinh\left(\frac{\hbar\beta\omega_{\boldsymbol{q}\nu}}{2}\right)}}$
(11)

We then directly compare eqs. 6 and 7 and as such we simplify the expression
under the square root of eq. 11 by taking its power series around
$G_{\boldsymbol{q}\nu}=0$ to first order to obtain:

$\displaystyle\gamma(\hbar\beta\bar{\omega}_{\boldsymbol{q}\nu})$
$\displaystyle\approx\frac{\hbar\beta\omega_{\boldsymbol{q}\nu}}{4G_{\boldsymbol{q}\nu}\sinh^{2}(\hbar\beta\omega_{\boldsymbol{q}\nu}/2)}$
(12)

Substituting eq. 12 into eq. 7, we see that the two rates are identically
equal. As such, we conclude that any additional (positive) scattering terms
arising from terms higher than linear order in $G_{\boldsymbol{q}\nu}$ can
only decrease the phonon-mediated electron-electron lifetime, showing that
$\tau_{\rm ee(ph)}$ is strictly smaller than $\tau_{\rm e-ph}$.

Equations 11 and 12 are most valid for dimensionless electron-phonon couplings
$G_{\boldsymbol{q}\nu}\ll 1$, which is true for non-superconducting materials
systems. Further, the phonon-mediated electron-electron interaction proceeds
via the exchange of a ‘virtual’ phonon, i.e. the phonon emitted (absorbed) by
a pair of electrons is assumed to be instantaneously absorbed (emitted) by a
different pair of electrons. In practice, this means that other phonon
scattering mechanisms, such as anharmonic phonon-phonon scattering, must be
‘quenched’ out, that is, occur at a much slower rate. This assumption is
justified at low temperatures, but breaks down at temperatures approaching the
material’s Debye temperature. At those elevated temperatures, one can still
use eq. 11 by augmenting $G_{\boldsymbol{q}\nu}$ with a competing coupling
$G^{pp}_{\boldsymbol{q}\nu}$ capturing anharmonic phonon-phonon interactions.
Integrating eq. 7 numerically using this competing coupling, we find that
$\tau_{\rm ee(ph)}$ naturally increases, and can overtake $\tau_{\rm eph}$ at
high temperatures. In both ZrSiS and TaAs2, (as well as WP2 and WTe2), there
is no energy separation between the acoustic branches and the lowest-energy
optical modes. These lower-energy optical modes involve large displacements of
the heavy Zr/Ta atoms. The heavy transition metal atoms will result in a
narrow bandwidth of the acoustic branches and thus a smaller phase space and
long phonon-phonon lifetimes, as discussed in high thermal-conductivity system
cubic boron arsenide [20].

This analytical limit allows us to revisit and further establish our design
principles discussed in Sec.section II. To obtain a low enough
$l_{\mathrm{mc}}$, we need a material with non-vanishing density of states at
the Fermi level, and a large electron-phonon matrix element. While the former
observation encourages us to seek for short $\tau_{\rm ee(ph)}$ in metals, the
latter indicates that the electron potential is strongly sensitive to lattice
perturbations, i.e. phonons. This is the case when the atomic orbital mixing
has low symmetry (i.e. the band extrema are located off high symmetry points
in the Brillouin zone), where the electron wavefunction, or equivalently the
charge distribution in real space, is highly anisotropic. To maintain a long
enough $l_{\rm mr}$, large Fermi velocity facilitated by dispersive electronic
bands can be combined with long electron-phonon lifetime maintained by reduced
electronic phase space at low temperatures. Indeed, the electronic structures
of both ZrSiS and TaAs2 explored in this work, as well as WTe2 and WP2
reported earlier, display the above features. There, the Fermi surfaces arise
from a combination of $d$ orbitals from the transition metal and $p$ orbitals
from the metalloids. These material-specific insights are poised to inform
exploration hydrodynamic electron flow in various condensed matter systems.

##  V Conclusions and Outlook

This _Article_ identifies material signatures for hydrodynamic electron flow
in semimetals from first principles, and predicts two candidates with
prominent hydrodynamic effects across a wide temperature range that warrant
experimental exploration. We obtain the momentum-relaxing electron-phonon mean
free paths $l_{\rm mr}$ and the momentum-conserving phonon mediated electron-
electron mean free paths $l_{\rm mc}$ in ZrSiS and TaAs2, and show that the
indirect phonon-mediated electron-electron interaction could dominate over a
wide temperature range, facilitating hydrodynamic behavior. Through inspection
of the momentum relaxing lifetimes on their Fermi surfaces, we find that at
low temperatures, hole pockets with open Fermi surface shape feature much
longer lifetimes than electron pockets. These observations suggest that
topological Dirac/Weyl bands are not indispensable in the search for promising
candidates for hydrodynamic flow, but highly dispersive bands are beneficial
in realizing long $l_{\rm mr}$. We review these signatures in light of
previously studied hydrodynamic candidates and discuss the analytical limits
of the first- and second-order electron-phonon interactions to find that the
mechanism is more broadly applicable. At the Eliashberg level of theory, we
then propose general principles for experimental discovery of hydrodynamic
behavior in anisotropic metals, such as low symmetry crystals with $d/p$
atomic orbital mixing, suppressed phonon-phonon scattering, and reduced
electronic phase space. Future work distinguishing electron-electron Umklapp
scattering events and computing the electron viscosity tensor from first
principles is still needed to provide a full picture of electron hydrodynamics
in quantum materials.

##  VI Acknowledgement

The authors acknowledge fruitful discussions with Adam Jermyn, Uri Vool, Doug
Bonn, Ken Burch, and Amir Yacoby. This work (Y.W., G.V., P.N.) was primarily
supported by the Quantum Science Center (QSC), a National Quantum Information
Science Research Center of the U.S. Department of Energy (DOE). P.N. is a
Moore Inventor Fellow and gratefully acknowledges support through Grant No.
GBMF8048 from the Gordon and Betty Moore Foundation. This research used
resources of the Oak Ridge Leadership Computing Facility, which is a DOE
Office of Science User Facility supported under Contract DE-AC05-00OR22725 as
well as the resources of the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

## References

  * Bandurin  _et al._ [2018] D. A. Bandurin, A. V. Shytov, L. S. Levitov, R. K. Kumar, A. I. Berdyugin, M. B. Shalom, I. V. Grigorieva, A. K. Geim, and G. Falkovich, Nature communications 9, 1 (2018).
  * Kotov  _et al._ [2012] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012).
  * de Jong and Molenkamp [1995a] M. J. M. de Jong and L. W. Molenkamp, Phys. Rev. B 51, 13389 (1995a).
  * Bandurin  _et al._ [2016] D. A. Bandurin, I. Torre, R. K. Kumar, M. B. Shalom, A. Tomadin, A. Principi, G. H. Auton, E. Khestanova, K. S. Novoselov, I. V. Grigorieva, L. A. Ponomarenko, A. K. Geim, and M. Polini, Science 351, 1055 (2016).
  * Crossno  _et al._ [2016] J. Crossno, J. K. Shi, K. Wang, X. Liu, A. Harzheim, A. Lucas, S. Sachdev, P. Kim, T. Taniguchi, K. Watanabe, T. A. Ohki, and K. C. Fong, Science 351, 1058 (2016).
  * Kumar  _et al._ [2017] R. K. Kumar, D. A. Bandurin, F. M. D. Pellegrino, Y. Cao, A. Principi, H. Guo, G. H. Auton, M. B. Shalom, L. A. Ponomarenko, G. Falkovich, K. Watanabe, T. Taniguchi, I. V. Grigorieva, L. S. Levitov, M. Polini, and A. K. Geim, Nature Physics 13, 1182 (2017).
  * Gooth  _et al._ [2018] J. Gooth, F. Menges, N. Kumar, V. Su$\mathrm{\beta}$, C. Shekhar, Y. Sun, U. Drechsler, R. Zierold, C. Felser, and B. Gotsmann, Nature Communications 9 (2018).
  * Vool  _et al._ [2020] U. Vool, A. Hamo, G. Varnavides, Y. Wang, T. X. Zhou, N. Kumar, Y. Dovzhenko, Z. Qiu, C. A. Garcia, A. T. Pierce, _et al._ , arXiv preprint arXiv:2009.04477 (2020).
  * Nandi  _et al._ [2018] N. Nandi, T. Scaffidi, P. Kushwaha, S. Khim, M. E. Barber, V. Sunko, F. Mazzola, P. D. C. King, H. Rosner, P. J. W. Moll, M. König, J. E. Moore, S. Hartnoll, and A. P. Mackenzie, npj Quantum Materials 3 (2018).
  * Moll  _et al._ [2016] P. J. Moll, P. Kushwaha, N. Nandi, B. Schmidt, and A. P. Mackenzie, Science 351, 1061 (2016).
  * Varnavides  _et al._ [2021] G. Varnavides, Y. Wang, P. J. Moll, P. Anikeeva, and P. Narang, arXiv preprint arXiv:2106.00697 (2021).
  * Levitov and Falkovich [2016] L. Levitov and G. Falkovich, Nature Physics 12, 672 (2016).
  * de Jong and Molenkamp [1995b] M. J. M. de Jong and L. W. Molenkamp, Physical Review B 51, 13389 (1995b).
  * Lucas and Fong [2018] A. Lucas and K. C. Fong, Journal of Physics: Condensed Matter 30, 053001 (2018).
  * Ku  _et al._ [2020] M. J. H. Ku, T. X. Zhou, Q. Li, Y. J. Shin, J. K. Shi, C. Burch, L. E. Anderson, A. T. Pierce, Y. Xie, A. Hamo, U. Vool, H. Zhang, F. Casola, T. Taniguchi, K. Watanabe, M. M. Fogler, P. Kim, A. Yacoby, and R. L. Walsworth, Nature 583, 537 (2020).
  * Jenkins  _et al._ [2020] A. Jenkins, S. Baumann, H. Zhou, S. A. Meynell, D. Yang, K. Watanabe, T. Taniguchi, A. Lucas, A. F. Young, and A. C. B. Jayich, arXiv:2002.05065 (2020).
  * Coulter  _et al._ [2018] J. Coulter, R. Sundararaman, and P. Narang, Phys. Rev. B 98, 115130 (2018).
  * Osterhoudt  _et al._ [2021] G. B. Osterhoudt, Y. Wang, C. A. C. Garcia, V. M. Plisson, J. Gooth, C. Felser, P. Narang, and K. S. Burch, Phys. Rev. X 11, 011017 (2021).
  * Butler [1977] W. Butler, Physical Review B 15, 5267 (1977).
  * Lindsay  _et al._ [2013] L. Lindsay, D. A. Broido, and T. L. Reinecke, Phys. Rev. Lett. 111, 025901 (2013).
  * Giuliani and Quinn [1982] G. F. Giuliani and J. J. Quinn, Phys. Rev. B 26, 4421 (1982).
  * Kaveh and Wiser [1984] M. Kaveh and N. Wiser, Advances in physics 33, 257 (1984).
  * Wilson [1938] A. H. Wilson, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 167, 580 (1938).
  * van Delft  _et al._ [2021] M. R. van Delft, Y. Wang, C. Putzke, J. Oswald, G. Varnavides, C. A. Garcia, C. Guo, H. Schmid, V. Süss, H. Borrmann, _et al._ , Nat. Commun. 12, 1 (2021).
  * Ong [1991] N. P. Ong, Physical Review B 43, 193 (1991).
  * Wang and Narang [2020] Y. Wang and P. Narang, Phys. Rev. B 102, 125122 (2020).
  * Ali  _et al._ [2014] M. N. Ali, J. Xiong, S. Flynn, J. Tao, Q. D. Gibson, L. M. Schoop, T. Liang, N. Haldolaarachchige, M. Hirschberger, N. P. Ong, and R. J. Cava, Nature 514, 205 (2014).
  * Pletikosić  _et al._ [2014] I. Pletikosić, M. N. Ali, A. V. Fedorov, R. J. Cava, and T. Valla, Phys. Rev. Lett. 113, 216601 (2014).
  * Kumar  _et al._ [2016] N. Kumar, C. Shekhar, S.-C. Wu, I. Leermakers, O. Young, U. Zeitler, B. Yan, and C. Felser, Phys. Rev. B 93, 241106 (2016).
  * Zhou  _et al._ [2019] W. Zhou, B. Li, C. Q. Xu, M. R. Delft, Y. G. Chen, X. C. Fan, B. Qian, N. E. Hussey, and X. Xu, Advanced Electronic Materials 5, 1900250 (2019).
  * Ali  _et al._ [2016] M. N. Ali, L. M. Schoop, C. Garg, J. M. Lippmann, E. Lara, B. Lotsch, and S. S. P. Parkin, Science Advances 2, e1601742 (2016).
  * Lv  _et al._ [2016] Y.-Y. Lv, B.-B. Zhang, X. Li, S.-H. Yao, Y. B. Chen, J. Zhou, S.-T. Zhang, M.-H. Lu, and Y.-F. Chen, Applied Physics Letters 108, 244101 (2016).
  * Chiu  _et al._ [2019] Y.-C. Chiu, K.-W. Chen, R. Schönemann, V. L. Quito, S. Sur, Q. Zhou, D. Graf, E. Kampert, T. Förster, K. Yang, G. T. McCandless, J. Y. Chan, R. E. Baumbach, M. D. Johannes, and L. Balicas, Phys. Rev. B 100, 125112 (2019).
  * Zhang  _et al._ [2016] C.-L. Zhang, S.-Y. Xu, I. Belopolski, Z. Yuan, Z. Lin, B. Tong, G. Bian, N. Alidoust, C.-C. Lee, S.-M. Huang, T.-R. Chang, G. Chang, C.-H. Hsu, H.-T. Jeng, M. Neupane, D. S. Sanchez, H. Zheng, J. Wang, H. Lin, C. Zhang, H.-Z. Lu, S.-Q. Shen, T. Neupert, M. Z. Hasan, and S. Jia, Nature Communications 7 (2016).
  * Arnold  _et al._ [2016] F. Arnold, C. Shekhar, S.-C. Wu, Y. Sun, R. D. dos Reis, N. Kumar, M. Naumann, M. O. Ajeesh, M. Schmidt, A. G. Grushin, J. H. Bardarson, M. Baenitz, D. Sokolov, H. Borrmann, M. Nicklas, C. Felser, E. Hassinger, and B. Yan, Nature Communications 7 (2016).
  * Li  _et al._ [2016] H. Li, H. He, H.-Z. Lu, H. Zhang, H. Liu, R. Ma, Z. Fan, S.-Q. Shen, and J. Wang, Nature Communications 7 (2016), 10.1038/ncomms10301.
  * Luo  _et al._ [2016] Y. Luo, R. D. McDonald, P. F. S. Rosa, B. Scott, N. Wakeham, N. J. Ghimire, E. D. Bauer, J. D. Thompson, and F. Ronning, Scientific Reports 6 (2016), 10.1038/srep27294.
  * Diez  _et al._ [2015] M. Diez, A. M. R. V. L. Monteiro, G. Mattoni, E. Cobanera, T. Hyart, E. Mulazimoglu, N. Bovenzi, C. W. J. Beenakker, and A. D. Caviglia, Phys. Rev. Lett. 115, 016803 (2015).
  * Uykur  _et al._ [2019] E. Uykur, L. Z. Maulana, L. M. Schoop, B. V. Lotsch, M. Dressel, and A. V. Pronin, Phys. Rev. Research 1, 032015 (2019).
  * Schoop  _et al._ [2016] L. M. Schoop, M. N. Ali, C. Straßer, A. Topp, A. Varykhalov, D. Marchenko, V. Duppel, S. S. P. Parkin, B. V. Lotsch, and C. R. Ast, Nature Communications 7 (2016).
  * Fu  _et al._ [2019] B.-B. Fu, C.-J. Yi, T.-T. Zhang, M. Caputo, J.-Z. Ma, X. Gao, B. Q. Lv, L.-Y. Kong, Y.-B. Huang, P. Richard, M. Shi, V. N. Strocov, C. Fang, H.-M. Weng, Y.-G. Shi, T. Qian, and H. Ding, Science Advances 5, eaau6459 (2019).
  * Hu  _et al._ [2017] J. Hu, Z. Tang, J. Liu, Y. Zhu, J. Wei, and Z. Mao, Phys. Rev. B 96, 045127 (2017).
  * Singha  _et al._ [2017] R. Singha, A. K. Pariari, B. Satpati, and P. Mandal, Proceedings of the National Academy of Sciences 114, 2468 (2017).
  * Schilling  _et al._ [2017] M. Schilling, L. Schoop, B. Lotsch, M. Dressel, and A. Pronin, Physical Review Letters 119 (2017).
  * Butcher  _et al._ [2019] T. A. Butcher, J. Hornung, T. Förster, M. Uhlarz, J. Klotz, I. Sheikin, J. Wosnitza, and D. Kaczorowski, Phys. Rev. B 99, 245112 (2019).
  * Wu  _et al._ [2016] D. Wu, J. Liao, W. Yi, X. Wang, P. Li, H. Weng, Y. Shi, Y. Li, J. Luo, X. Dai, and Z. Fang, Applied Physics Letters 108, 042105 (2016).
  * Wang  _et al._ [2016] Y.-Y. Wang, Q.-H. Yu, P.-J. Guo, K. Liu, and T.-L. Xia, Phys. Rev. B 94, 041103 (2016).
  * Yuan  _et al._ [2016] Z. Yuan, H. Lu, Y. Liu, J. Wang, and S. Jia, Phys. Rev. B 93, 184405 (2016).
  * Giustino [2017] F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
  * Allen [1978] P. B. Allen, Phys. Rev. B 17, 3725 (1978).


# Scalar field perturbation around a rotating hairy black hole: quasinormal
modes, quasibound states and superradiant instability

Yun-He Lei<EMAIL_ADDRESS>Center for Gravitation and Cosmology, College
of Physical Science and Technology, Yangzhou University, Yangzhou, 225009,
China  Zhen-Hao Yang<EMAIL_ADDRESS>Center for Gravitation and
Cosmology, College of Physical Science and Technology, Yangzhou University,
Yangzhou, 225009, China  Xiao-Mei Kuang<EMAIL_ADDRESS>(corresponding
author) Center for Gravitation and Cosmology, College of Physical Science and
Technology, Yangzhou University, Yangzhou, 225009, China

###### Abstract

We consider the quasinormal modes, quasibound states and superradiant
instability of a rotating hairy black hole, which possesses a Horndeski hair
as deviation from Kerr black hole, under the perturbation of massive scalar
field. With the use of the matrix method, we mainly calculate the
eigenfrequencies related to those modes of the perturbation. Under the
perturbation of the massless scalar field, the Horndeski hair and spin
parameter have significant influences on the quasinormal frequency, but its
imaginary part is always finite negative and no unstable mode is found. Under
the perturbation of the massive scalar field, we focus on the eigenfrequencies
of quasibound states and find the modes of which the imaginary part of
eigenfrequencies is positive, indicating that the black hole undergoes
superradiant instability. Then we scan the parameters and figure out a diagram
in the space of Horndeski hair and spin parameters to distinguish the rotating
hairy black hole with superradiant instability from the stable one.

###### Contents

  1. I Introduction
  2. II Background geometry and setup of the scalar field perturbation
  3. III Results of the perturbing modes and their eigenfrequencies
    1. III.1 Quasinormal modes
    2. III.2 Quasibound states and superradiant instabilities
  4. IV Conclusions

##  I Introduction

The stability of black hole solution under certain perturbations is an over-
lasting crucial problem in black hole physics. Quasinormal modes (QNMs) with
specific eigenfrequencies (Quasinormal frequencies, or QNFs) play a remarkable
role in this scenario. QNFs of a field perturbation on the black hole are
infinite discrete spectrum of complex frequencies, of which the real part
determines the oscillation timescale of the QNMs, while the imaginary part
determines their exponential decaying timescale. It is one of the most
important characteristics of black hole geometry, and the interest of QNMs in
more fundamental physics can be referred into the reviews Nollert (1999);
Berti et al. (2009); Konoplya and Zhidenko (2011). Early in the 1950s, Regge
and Wheeler studied the stability of black holes in general relativity (GR) by
analyzing QNFs of gravitational perturbation in Schwarzschild black hole, and
they found that the perturbation field oscillates and decays over time such
that the Schwarzschild black hole is stable Regge and Wheeler (1957). Later,
various external field perturbations in the Schwarzschild black hole were
explored and further verified the stability of the Schwarzschild black hole in
GR Konoplya (2003); Cardoso et al. (2004); Cho (2003); Jing (2005) .

However, rotating neutral black holes(i.e. Kerr black holes) in GR may develop
instabilities under certain conditions. In fact, if an ingoing bosonic wave is
scattered off a Kerr black hole with the horizon angular velocity
$\Omega_{H}$, the scattered wave will be amplified once the superradiance
condition $0<\omega<m\Omega_{H}$ is satisfied, where $\omega$ is the frequency
of the wave and $m$ is the azimuthal number Starobinsky (1973). The extra
energy of the scattered wave is obtained from the rotational energy of the
black hole, so this is a wave analogue of the Penrose process Penrose and
Floyd (1971); Misner (1972). In this process, if there exists a hypothetic
reflective “mirror” that makes the amplified wave be scattered back and forth
between the “mirror” and the black hole, the amplifying proceeds continuously
and will finally cause the black hole background to suffer from superradiant
instability Brito et al. (2015). Such a system was also dubbed the “black hole
bomb” Press and Teukolsky (1972). The Kerr black hole was found to be stable
under massless scalar, electromagnetic or gravitational perturbations
according to the analysis of those QNMs Press and Teukolsky (1973); Teukolsky
and Press (1974). However, once the perturbing scalar field is massive, Kerr
black hole was found to undergo superradiant instability because the mass term
of the perturbing scalar field acts as the natural mirror Damour et al.
(1976); Zouros and Eardley (1979); Detweiler (1980); Li (2012); Dolan (2007);
Cardoso and Yoshida (2005); Hod and Hod (2010); Damour and Ruffini (1976).
Moreover, it was addressed in Brito et al. (2015); Dolan (2007) that such
superradiatively unstable mode may be induced from the quasi-bound state (QBS)
of the massive scalar field, which decays exponentially far away from the
black hole with certain eigenfrequency.

In GR, the eigenfrequencies for QNM and QBS around a Kerr black hole are only
characterized by three parameters: the mass and angular momentum of the black
hole and the mass of external field. Though recent observations of
gravitational wave generated from binary compact objects Abbott et al. (2016,
2019) and black hole shadow Akiyama et al. (2019, 2022) provide us with
remarkable chances to test GR in the strong field regime, the uncertainties in
the data leave some space for alternative theories of gravity and non-Kerr
metric of black holes. In addition, despite the great success, there are still
some challenging problems for GR, such as the explanation of the Universe
expansion history, the large scale structure and the understanding of quantum
gravity Riess et al. (1998); Weinberg (1989); Carroll (2001). Therefore, a
more general theory of gravity is required from both observational and
theoretical aspects. Plenty of modified gravitational theories have been
proposed Clifton et al. (2012); Bakopoulos (2020); Corelli (2020), among which
the usual way is to modify the action of GR. The scalar tensor theories, which
contains a scalar field in addition to the metric field are considered as the
simplest nontrivial extensions of GR. The most general scalar tensor theory
with second order equations of motion for both the metric and the scalar field
in four dimensions is the Horndeski theory, which is free of the Ostrogradski
instabilities Horndeski (1974); Kobayashi (2019) and is found to be equivalent
to the generalized Galileon theory Deffayet et al. (2011); Kobayashi et al.
(2011) derived in a different way. The observational constraints or bounds on
Horndeski gravity have been extensively explored in Bellini et al. (2016);
Bhattacharya and Chakraborty (2017); Kreisch and Komatsu (2018); Hou and Gong
(2018); Spurio Mancini et al. (2019); Allahyari et al. (2020). Moreover,
Horndeski gravity attracts lots of attention in both cosmological and
astrophysical communities because it has significant consequences in
describing the accelerated expansion and other interesting features (please
see Kobayashi (2019) for a review), as well as holographic applications Feng
et al. (2015); Kuang and Papantonopoulos (2016); Jiang et al. (2017); Baggioli
and Li (2017); Feng and Liu (2019); Wang et al. (2019); Zhang et al. (2023);
Bravo-Gaete and Santos (2022); Bravo-Gaete et al. (2022) and references
therein.

In particular, the Horndeski framework has also been used to test the no hair
theorem of black holes, which states that the isolated black hole in GR is
merely characterized by three parameters: mass, electric charge and angular
momentum. Due to the extra scalar field in the action, it is natural to figure
out if there are black hole solutions with scalar hair, and this makes sense
because similar to GR, the Horndeski theory has diffeomorphism invariance and
second order field equations. The no hair theorem for the shift symmetric
Horndeski theory Hui and Nicolis (2013) has proved that a static and
spherically symmetric black hole cannot sustain a nontrivial scalar field with
vanishing boundary condition at infinity. But later it turns out that there is
a loophole in the demonstration and the counter hairy black hole solution was
constructed in Sotiriou and Zhou (2014a, b). The no hair theorem has also been
considered in Babichev et al. (2017) for the shift symmetric Horndeski and
beyond Horndeski theories, which also admit the hairy static and spherically
symmetric black hole solutions with nontrivial scalar profile. Up to date,
many hairy black holes have been constructed and analyzed in Horndeski theory,
including the radially dependent Rinaldi (2012); Cisterna and Erices (2014);
Feng et al. (2015) and time–dependent scalar hair Babichev and Charmousis
(2014); Babichev et al. (2018); Ben Achour and Liu (2019), among which the
black holes with linear time–dependent scalar hair is unstable under
perturbation Khoury et al. (2020), hence ruling out the possibility of this
type of hairy solutions. Fortunately, a spherically symmetric hairy black hole
solution in shift symmetric Horndeski theory has recently been constructed in
Bergliaffa et al. (2021), which was verified to be stable under various
massless external field perturbations by us Yang et al. (2023). Some
observational investigations have been carried out for this static hairy
Horndeski black hole, for example, the strong gravitational lensing Kumar et
al. (2022) and photon rings in the black hole image Wang et al. (2023), which
shew that the Horndeski hair has a significant influence on the related
observable. Especially, with the use of the revised Newman–Janis algorithm
Azreg-Aïnou (2014a, b), the authors of Walia et al. (2022) soon constructed
its axially symmetric counterpart, known as rotating black hole with Horndeski
hair. The thermodynamic properties and weak gravitational lensing Walia et al.
(2022), black hole shadow Afrin and Ghosh (2022) and superradiant energy
extraction Jha et al. (2022) have been extensively studied in this rotating
hairy black hole.

In this paper, we will study the QNMs and QBSs for a massive perturbing scalar
field on the rotating hairy black hole in Horndeski theory. Our study
motivates from two aspects. One is that the computation of eigenfrequencies of
QNMs for scalar field perturbation is a first step to analyze the dynamical
stability of the rotating hairy black hole. The other is that the interesting
QBSs could be the potential source of the superradiatively unstable modes as
we aforementioned. It is noticed that a brief investigation of superradiant
instabilities in the rotating hairy metric was performed in Jha et al. (2022)
using approximate analytical methods. Our work expands the results of Jha et
al. (2022) by analyzing in detail the eigenfrequency of these unstable modes
as a function of the black hole parameters. We will explore how the hair
parameter influences this instability apart from the rotation. To this end, we
will use the matrix method to numerically determine the frequency spectrums of
QNMs, QBSs and superradiatively unstable modes.

This paper is organized as follows. In Sec.II we present the background
geometry of the rotating hairy Horndeski black hole, derive the master
equations of the scalar field perturbation and analyze the corresponding
boundary conditions. The eigenfrequencies of QNMs, QBSs and superradiant
instabilities of the perturbation are explored in Sec.III, in which we mainly
analyze the effect of the hair parameter and figure out a ‘phase’ diagram in
black hole parameter space. Sec.IV contains our conclusions.

##  II Background geometry and setup of the scalar field perturbation

The action for the dubbed quartic Horndeski gravity reads as Babichev et al.
(2017)

$\displaystyle S=\int
d^{4}x\sqrt{-g}\big{[}Q_{2}+Q_{3}\Box\phi+Q_{4}R+Q_{4,\chi}\left((\Box\phi)^{2}-(\nabla^{\mu}\nabla^{\nu}\phi)(\nabla_{\mu}\nabla_{\nu}\phi)\right)+Q_{5}G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi$
$\displaystyle-\frac{1}{6}Q_{5,\chi}\left((\Box\phi)^{3}-3(\Box\phi)(\nabla^{\mu}\nabla^{\nu}\phi)(\nabla_{\mu}\nabla_{\nu}\phi)+2(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\nu}\nabla^{\gamma}\phi)(\nabla_{\gamma}\nabla^{\mu}\phi)\right)\big{]},$
(1)

where $\chi=-\partial^{\mu}\phi\partial_{\mu}\phi/2$ is the canonical kinetic
term, $Q_{i}~{}(i=2,3,4,5)$ are arbitrary functions of $\chi$ and
$Q_{i,\chi}\equiv\partial Q_{i}/\partial\chi$, $R$ is the Ricci scalar and
$G_{\mu\nu}$ is the Einstein tensor. A static hairy black hole in a specific
quartic Horndeski theory, saying that $Q_{5}$ in the above action vanishes,
has been constructed in Bergliaffa et al. (2021) with the metric

$\displaystyle
ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2})~{}~{}~{}\mathrm{with}~{}~{}~{}f(r)=1-\frac{2M}{r}+\frac{h}{r}\ln\left(\frac{r}{2M}\right).$
(2)

Here, $M$ and $h$ are the parameters related to the black hole mass and scalar
hair. The metric reduces to Schwarzschild case as $h\to 0$, and it is
asymptotically flat. The hairy black hole described by the above metric was
found to be stable under various external field perturbations by computing
their QNFs and time evolutions Yang et al. (2023). Besides, some observable
phenomena affected by the Horndeski hair have also been explored in Kumar et
al. (2022); Wang et al. (2023). Very recently, the axial symmetric counterpart
of metric (2) was constructed, and in the Boyer-Lindquist coordinates the
metric of the corresponding rotating hairy black hole reads as Walia et al.
(2022)

$ds^{2}=-\left[1-\frac{2\tilde{M}(r)r}{\Sigma}\right]dt^{2}+\frac{\Sigma}{\Delta}dr^{2}+\Sigma
d\theta^{2}+\frac{(r^{2}+a^{2})^{2}-a^{2}\Delta\sin^{2}\theta}{\Sigma}\sin^{2}\theta
d\varphi^{2}-\frac{4a\tilde{M}(r)r}{\Sigma}\sin^{2}\theta dtd\varphi,$ (3)

with the functions

$\Delta=r^{2}-2\tilde{M}(r)r+a^{2},\quad\Sigma=r^{2}+a^{2}\cos^{2}\theta,\quad\tilde{M}(r)=M-\frac{h}{2}\ln\left(\frac{r}{2M}\right),$
(4)

where $a$ is the spin parameter. The axially symmetric spacetime is also
asymptotically flat and reduces to the Kerr spacetime as $h\to 0$. The
horizons of the rotating hairy black hole are determined by the real positive
roots of $\Delta=0$, which only exist in a certain parameter region of
$(a/M,h/M)$ shown in FIG.1. In detail, in the gray parameter region,
$\Delta=0$ has two real positive roots $r_{\pm}$, in which the larger root
$r_{+}$ indicates the event horizon while the smaller one $r_{-}$ denotes the
Cauchy horizon. In the white parameter region, $\Delta=0$ has no real root
meaning there is no horizon such that the metric (3) describes a naked
singularity. And the parameters on the black curve that separates these two
regions make the rotating hairy black hole become extremal with $r_{-}=r_{+}$.
The horizons $r_{\pm}$ as functions of the black hole parameters are
explicitly shown in FIG.2. More interesting properties, including the Komar
conserved quantities, thermodynamics, light deflection, shadow, energy
emission and energy extraction, of this rotating hairy black hole can be found
in Walia et al. (2022); Afrin and Ghosh (2022); Jha et al. (2022).

Figure 1: The $(a/M,h/M)$ parameter regions in which the metric (3) describes
a black hole (gray region) or a naked singularity (white region), the
parameters on the black curve that separates these two regions correspond to a
extreme black hole.

Figure 2:  The horizons $r_{\pm}$ as functions of the black hole parameters.
The orange part indicates the behavior of event horizon $r_{+}$ while the blue
part describe the behavior of Cauchy horizon $r_{-}$.

Here in order to analyze the dynamical (in)stability of the rotating hairy
black hole (3), we, as the first step, will introduce the massive scalar field
as the probe and study the eigenfrequencies of its master equation. The
dynamic of a perturbing scalar field around the rotating hairy black hole is
governed by the Klein–Gordon equation

$(\Box-\mu^{2})\Phi=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\Phi)-\mu^{2}\Phi=0,$
(5)

where $\mu$ is the mass of the scalar field. Then with the decomposition

$\Phi(t,r,\theta,\varphi)=R(r)S(\theta)e^{-i\omega t+im\varphi},$ (6)

in which $\omega$ is the frequency and $m$ is the azimuthal number, we can
separate Eq.(5), yielding the radial master equation

$\Delta\frac{d}{dr}\left[\Delta\frac{d}{dr}R(r)\right]+[\omega^{2}(a^{2}+r^{2})^{2}-\Delta(a^{2}\omega^{2}+\mu^{2}r^{2}+\lambda)-4a\omega
mr\tilde{M}(r)+a^{2}m^{2}]R(r)=0,$ (7)

and the angular master equation

$\frac{d}{du}\left[(1-u^{2})\frac{d}{du}S(u)\right]+\left[a^{2}u^{2}(\omega^{2}-\mu^{2})-\frac{m^{2}}{1-u^{2}}+\lambda\right]S(u)=0,$
(8)

where $u=\cos\theta$ and $\lambda$ is the separation constant.

Then, we shall analyze the physically relevant boundary conditions of Eq.(7)
and Eq.(8). The boundary conditions of $S(u)$ can be directly obtained from
its asymptotic solutions near the boundaries $u\to\pm 1$, which are

$S(u)\sim\begin{cases}(1-u)^{\pm\frac{m}{2}},\quad u\to 1\\\
(1+u)^{\pm\frac{m}{2}},\quad u\to-1\end{cases}.$ (9)

For the radial equation, by introducing the tortoise coordinate

$r_{*}=\int\frac{r^{2}+a^{2}}{\Delta}dr,$ (10)

and the function transformation

$K(r)=\sqrt{r^{2}+a^{2}}R(r),$ (11)

we can obtain the equation

$\frac{d^{2}}{dr_{*}^{2}}K(r)+V(r)K(r)=0,$ (12)

with

$V(r)=\left(\omega-\frac{ma}{r^{2}+a^{2}}\right)^{2}-\frac{\Delta}{(r^{2}+a^{2})^{2}}\left[a^{2}\omega^{2}-2am\omega+\mu^{2}r^{2}+\lambda+\frac{(a^{2}-2r^{2})\Delta+r(a^{2}+r^{2})d\Delta/dr}{(r^{2}+a^{2})^{2}}\right].$
(13)

It is straightforward to get the asymptotic behavior of $V(r)$ as

$V(r\to
r_{+})=(\omega-m\Omega_{H})^{2}~{}~{}~{}\text{and}~{}~{}~{}V(r\to\infty)=\omega^{2}-\mu^{2},$
(14)

with $\Omega_{H}=\frac{a}{r_{+}^{2}+a^{2}}$ the horizon angular velocity,
based on which we can easily obtain the asymptotic solutions as the boundary
conditions

$K(r\to r_{+})\sim e^{\pm
i\bar{\omega}r_{*}}~{}~{}~{}\text{and}~{}~{}~{}K(r\to\infty)\sim e^{\pm
ikr_{*}}.$ (15)

with $\bar{\omega}=\omega-m\Omega_{H}$ and $k=\sqrt{\omega^{2}-\mu^{2}}$.
Reminding that the perturbation cannot come out from the event horizon, we
only consider the case with the minus sign in the first formula of Eq.(15),
i.e., the ingoing boundary condition near the event horizon. For the second
formula of Eq.(15), we take the plus (minus) sign and ensure $Re(k)>0$
($Im(k)<0$) to get the purely outgoing (exponentially decaying) boundary
condition at the spatial infinity for the QNMs (QBSs) Dolan (2007). The
eigenfrequencies $\omega$ obtained from the master equations are discrete
complex numbers, of which the real part describes the oscillation while the
imaginary part can be employed to judge the (in)stability of the black hole
under the perturbations. QNMs are known to have negative $Im(\omega)$, which
means that the perturbation decays and the lifetime is evaluated by
$\tau=|Im(\omega)|^{-1}$. But the QBSs may become the source of the unstable
modes with positive $Im(\omega)$, implying an exponentially growing amplitude
due to the superradiant instability. Therefore, we shall solve the master
equations to obtain the eigenfrequencies, and analyze the possible modes of
the perturbation in the rotating hairy black hole.

Since it’s not possible to solve the master equations analytically, so we
shall solve them numerically with the use of the matrix method Lin and Qian
(2017); Lin et al. (2017); Lin and Qian (2019) and extract the
eigenfrequencies of the QNMs and QBSs. For the convenience to employ this
method, we shall make some transformations to Eq.(7) and Eq.(8) beforehand.
Firstly, by rewriting the boundary conditions of the radial equation Eq.(15)
into the form of function $R$ via Eq.(11), and expanding them at the
corresponding asymptotic region, we can obtain

$R(r\to
r_{+})\sim(r-r_{+})^{-i\sigma}~{}~{}~{}\text{and}~{}~{}~{}R(r\to\infty)\sim
e^{\pm ikr\mp\frac{1}{2}ihk(\ln\frac{r}{2M})^{2}}r^{-1\pm 2iMk},$ (16)

where
$\sigma=\frac{r_{+}(r_{+}^{2}+a^{2})}{r_{+}(r_{+}+h)-a^{2}}\bar{\omega}$, and
the upper (lower) signs correspond to the plus (minus) sign in the second
formula of Eq.(15). Then, in order to factor out the asymptotic behaviors of
the radial and angular functions, we redefine

$\begin{split}R(r)&=(r-r_{+})^{-i\sigma}r^{i\sigma}e^{\pm
ikr\mp\frac{1}{2}ihk(\ln\frac{r}{2M})^{2}}r^{-1\pm 2iMk}\tilde{R}(r),\\\
S(u)&=(1-u)^{\frac{m}{2}}(1+u)^{\frac{m}{2}}\tilde{S}(u),\end{split}$ (17)

so that $R(r)$ and $S(u)$ will automatically satisfy the aforementioned
boundary conditions. Next, we introduce the new coordinates

$x=\begin{cases}1-(\frac{r_{+}}{r})^{1/3}\quad\text{for QNMs}\\\
1-\frac{r_{+}}{r}\quad\text{for
QBSs}\end{cases}~{}~{}~{}\text{and}~{}~{}~{}v=\frac{1+u}{2},$ (18)

to bring the integration domains into $x\in[0,1]$ and $v\in[0,1]$ to replace
$r\in[r_{+},\infty)$ and $u\in[-1,1]$, respectively. Finally, by further
performing the transformation

$\chi_{R}(x)=\begin{cases}x(1-x)\tilde{R}(x)\quad\text{for QNMs}\\\
x(1-x)^{6}\tilde{R}(x)\quad\text{for
QBSs}\end{cases}~{}~{}~{}\text{and}~{}~{}~{}\chi_{S}(v)=v(1-v)\tilde{S}(v),$
(19)

we obtain the equations of $\chi_{R}$ and $\chi_{S}$, which we can directly
solve by using the matrix method, as

$\displaystyle\mathscr{C}_{R2}(x)\chi_{R}^{\prime\prime}(x)+\mathscr{C}_{R1}(x)\chi_{R}^{\prime}(x)+\mathscr{C}_{R0}(x)\chi_{R}(x)=0,$
(20)
$\displaystyle\mathscr{C}_{S2}(v)\chi_{S}^{\prime\prime}(v)+\mathscr{C}_{S1}(v)\chi_{S}^{\prime}(v)+\mathscr{C}_{S0}(v)\chi_{S}(v)=0,$
(21)

where the coefficient functions $\mathscr{C}_{Ri}$ and $\mathscr{C}_{Si}$
($i=0\sim 2$) are functions of $\omega$ and $\lambda$. Moreover, the
corresponding boundary conditions are then converted into the homogeneous form

$\displaystyle\chi_{R}(0)=\chi_{R}(1)=0,$ (22)
$\displaystyle\chi_{S}(0)=\chi_{S}(1)=0.$ (23)

The last step is to follow the matrix method algorithm to extract the
eigenfrequencies of QNMs and QBSs by numerically solving the algebraic
equations

$det[\mathscr{M}_{R}(\omega,\lambda)]=0,\quad\text{and}\quad
det[\mathscr{M}_{S}(\omega,\lambda)]=0,$ (24)

where the expressions of the matrix
$\mathscr{M}_{R~{}\text{or}~{}S}(\omega,\lambda)$ are determined by the
coefficient functions in the equations (20) and (21). Here we prefer not to
repeat the steps of the matrix method algorithm, and readers can refer to Lin
and Qian (2017); Lin et al. (2017); Lin and Qian (2019) for more details.

##  III Results of the perturbing modes and their eigenfrequencies

In this section, we shall analyze the results of QNMs, QBSs and superradiant
instabilities of the perturbing scalar field on the rotating hairy black hole
background (3). It is noted that all the physical quantities in this section
are rescaled by the mass parameter $M$ into be dimensionless.

###  III.1 Quasinormal modes

We first consider the dynamic behavior of a massless ($\mu M=0$) scalar field
perturbing around the rotating hairy Horndeski black hole. We mainly focus on
the effects of spin parameter $a/M$ and the Horndeski hair parameter $h/M$ on
the fundamental QNMs with overtone number $N=0$, so we select samples of the
two parameters to show the dependence of the fundamental QNFs on them. The
fundamental QNFs with the mode $\ell=m=0$ are shown in FIG.3. With fixed $a/M$
(see the dashed curves), as the hair parameter $h/M$ decreases, the real part
of the QNFs decreases while the imaginary part increases. This means that
comparing to the Kerr black hole, the Horndeski hair makes the massless scalar
field perturbation oscillate slower and decay slower. With fixed $h/M$ (see
the solid curves), the effect of the spin parameter on QNFs could depend on
the value of $h/M$. In detail, for weaker Horndeski hair, i.e., small $\mid
h/M\mid$, as $a/M$ increases, the real part of QNFs first increases and then
decreases which is similar to that found in Kerr black hole Dolan (2007), but
for large enough $\mid h/M\mid$ it decreases monotonically. Meanwhile, the
imaginary part first increases with the increasing of $a/M$ and then decreases
slightly as the black hole approaches the extreme for all the selected values
of $h/M$, but it tends to a finite negative value for the extremal case. This
implies that for the rotating hairy black hole with small $\mid h/M\mid$, the
faster spin makes the corresponding perturbation oscillate first faster and
then slower as the black hole becomes near extremal, similar to that occurs in
Kerr black hole. However, in the rotating hairy black hole with large enough
$\mid h/M\mid$, the oscillation always becomes slower as $a/M$ increases, this
is because the strong suppression effect of $h/M$ on the real part could
balance the enhancement effect of $a/M$ and even dominates when the black hole
closes to extreme.

Figure 3: The variation of the real and imaginary parts of eigenfrequencies
for fundamental $\ell=m=0$ QNM of the massless scalar field on the rotating
hairy Horndeski black hole in terms of the hair (spin) parameter $h/M$ ($a/M$)
with several selected spin (hair) parameters. For each curve with given $a/M$
($h/M$), the value of $h/M$ ($a/M$) starts from $0$ (Kerr case) [(static
case)] and then extends toward the negative (positive) direction until the
black hole becomes near extreme.

Figure 4:  The fundamental QNFs for the modes $\ell=1$, with $m=-1$ (left
panel), $m=0$ (middle panel) and $m=1$ (right panel). The setting of black
hole parameters are the same as those in FIG.3.

To compare the influences of different angular quantum number $\ell$ and
azimuthal number $m$ on QNMs, we compute the fundamental QNFs for $\ell=1$.
The results are depicted in FIG. 4, where we exhibit how the QNFs of the modes
$m=-1$ (left panel), $m=0$ (middle panel) and $m=1$ (right panel) are affected
by $h/M$ and $a/M$ respectively, with the same selected values of parameters
as those in $\ell=0$ case. We can extract the following features. (i) As $h/M$
decreases, the real (imaginary) part of the fundamental QNFs decreases
(increases) with any selected $a/M$ for all the $\ell=1$ modes, except for the
mode $\ell=m=1$ with large $a/M$ of which the real part first decreases then
increases with the decreasing of $h/M$. (ii) The effect of varying $a/M$ can
also be read off. For the modes $m=0$ and $m=1$, both the real and the
imaginary parts of the QNFs increase as the spin becomes faster for the
selected $h/M$. For the mode $m=-1$, as the rotation increases, the real part
decreases monotonically for any selected $h/M$, but the imaginary part
increases for small $\mid h/M\mid$ and decreases for large $\mid h/M\mid$ on
the whole. (iii) According to our study, all the imaginary part of
eigenfrequencies are finite negative value and we do not find any unstable
mode.

###  III.2 Quasibound states and superradiant instabilities

We then turn on the mass of the scalar field and study the existence of QBSs
and their eigenfrequencies. As we mentioned in Sec.I, the mass term may
trigger the superradiant instability that we will show in this subsection. To
be specific, we fix the mass of the scalar field $\mu M=0.4$ as it was done in
Siqueira and Richartz (2022), and compute the fundamental eigenfrequencies of
QBSs for the modes $\ell=0$ and $\ell=1$ for the same selected parameters as
those in the QNMs. The results are depicted in FIG.5 and FIG.6, from which we
can find some common features of the modes. (i) The real part of the
frequencies exhibits a typical feature of the QBSs, i.e., they are always
smaller than the the mass of the scalar field. (ii) Similar to the QNFs, the
strength of the Horndeski hair also influences the effect of the spin
parameter on the QBS frequency, which could be different for different modes.
Here we will pay more attention to the effect of Horndeski hair on the
eigenfrequencies rather than starchily setting out the effect of rotation.

Figure 5:  The fundamental eigenfrequencies of the quasibound states in the
mode $\ell=m=0$ of perturbing scalar field with the mass $\mu M=0.4$. The
setting of black hole parameters are the same as those in QNMs.

Figure 6:  The fundamental eigenfrequencies of the quasibound states in the
modes $\ell=1$ of perturbing scalar field with the mass $\mu M=0.4$. From left
to right, we have $m=-1$, $m=0$ and $m=1$, respectively. Again, the setting of
black hole parameters are the same as those in previous studies.

Figure 7:  The fundamental eigenfrequency of the QBSs in the mode $\ell=m=1$
of massive ($\mu M=0.4$) scalar field as a function of $h/M$ with several
selected $a/M$. In the plots, each curve for fixed $a/M$ extends toward the
negative $h/M$ direction until the black hole becomes near extreme. In the
left plot, there exists positive values of the imaginary part, indicating the
superradiant instabilities. In the right plot, the solid curves denote the
numerical results of the real part. The red dots denote the eigenfrequencies
of the bound state or scalar clouds for which the imaginary part vanishes, so
the parameters are extracted from the intersections between each curve and the
dashed line in the left plot. The real part of these bound states marked by
red dots are directly evaluated by $\omega_{sc}=m\Omega_{H}$.

In FIG.5 with $\ell=0$, the imaginary part of the eigenfrequency is always
negative, indicating that the mode is stable. With a fixed spin parameter, the
stronger Horndeski hair suppresses both the real and imaginary parts of the
eigenfrequency, which is different from the one that happens in QNFs.
Moreover, it is obvious that by comparing FIG.5 and FIG.3, the real part of
the QBSs of the massive scalar field with $\mu M=0.4$ is larger than those of
the QNMs of the massless scalar field, indicating faster oscillation for the
former. In FIG.6 with $\ell=1$, with the selected $a/M$, the real and the
imaginary parts of QBSs from all modes always decrease with the decreasing of
$h/M$, however, the inset of the right panel shows the exception for the mode
$\ell=m=1$ with $a/M=0.8$, where the imaginary part first decreases slightly
then increases and eventually becomes positive. These QBSs with positive
$Im(\omega)$ means that the amplitude of these states will grow exponentially
over time, leading to the superradiant instabilities. Moreover, note that the
superradiant instabilities only exist in the $m>0$ states, as the result of
the superradiance condition $0<\omega<m\Omega_{H}$.

Figure 8:  The parameter regions for the stable rotating hairy Horndeski black
hole and the one with superradiant instability under the perturbation of
massive ($\mu M=0.4$) scalar field, of the fundamental $\ell=m=1$ QBSs. The
black curve corresponds to the extreme black hole as in FIG. 1, and the red
curve corresponds to the onset of superradiant instability, or bound state or
scalar cloud with the eigenfrequency $\omega_{sc}=m\Omega_{H}$ which is purely
real.

Figure 9: The fundamental eigenfrequency of the QBSs in the mode $\ell=m=1$ as
a function of $h/M$ for selected masses of the scalar field. Here we fix
$a/M=0.86$. In the left plot for imaginary part, we see that the onset of
superradiant instability is very sensitive to the scalar mass. Again, the red
dots in the right plot denote the eigenfrequencies of bound state or scalar
cloud, which we evaluate via $\omega_{sc}=m\Omega_{H}$.

To clearly show how the parameters of hairy black hole affect the superradiant
instabilities, we further plot the eigenfrequencies as a function of the
Horndeski hair for several large spin parameters in FIG. 7. In the left plot,
the dashed line with $Im(M\omega)=0$ separates the decaying modes with
$Im(M\omega)<0$ from the superradiatively unstable modes with $Im(M\omega)>0$.
In fact, the modes with $Im(M\omega)=0$ are bound states which could form as
the scalar cloud around the rotating hairy black hole. As addressed in Hod
(2012), the eigenfrequency of the scalar cloud is purely real
$\omega_{sc}=m\Omega_{H}$, which means that the scalar cloud is synchronized
with the rotation of the hairy black hole such that there is no energy flux
towards the black hole Hod (2012); Herdeiro and Radu (2014); Benone et al.
(2014). In the right plot, our numerical real part matches well with
$\omega_{sc}=m\Omega_{H}$ marked by the red dots for the scalar cloud. In a
word, the Horndeski hair significantly affects the superradiant instability
and will make the onset of superradiant instability need a smaller spin
parameter. We collect the data for the scalar cloud and superradiatively
unstable modes into the $(a/M,h/M)$ space, depicted in FIG.8 where we rescale
the axes, and the red curve denotes the onset of instability or the black hole
with a scalar cloud. The rotating hairy Horndeski black hole with the
parameters in the light gray region is stable under the perturbation of the
scalar field with mass $\mu M=0.4$, while the black hole with the parameters
in the dark gray region could undergo superradiant instability. It is shown
that comparing to the Kerr black hole $(h=0)$, the superradiant instability
can occur for slower rotating hairy Horndeski black hole but the parameter
region of the existence of superradiant instability is narrower.

Before closing this section, we shall study the effect of the scalar field
mass $\mu M$ on the QBSs and the superradiant instabilities. As addressed in
Dolan (2007), the superradiant instability of Kerr black hole occurs when the
mass is in the range of $0<\mu M<\mu_{c}M$, where $\mu_{c}M$ depends on the
spin of black hole and we name it as the upper bound of the scalar mass
corresponding to $Im(M\omega)=0$. Here we fix the spin parameter $a/M=0.86$
and calculate the eigenfrequencies of the fundamental $\ell=m=1$ mode as a
function of the Horndeski hair. Due to the technical problems, it is difficult
for us to scan all the values of $\mu M$, so we only give the results with
selected values of $\mu M$ in FIG. 9 . The figure shows that in the selected
parameter region, as $\mu M$ increases, the imaginary (real) part of the
eigenfrequencies decreases (increases). Moreover, for the rotating hairy
Horndeski black hole with a smaller $h/M$, the upper bound of $\mu M$ for
superradiant instability increases. This suggests that comparing to the Kerr
black hole, the mass of the scalar cloud for the rotating hairy Horndeski
black hole is enhanced. In addition, we evaluate the eigenfrequencies of the
scalar cloud via $\omega_{sc}=m\Omega_{H}$ marked by the red dots in the right
plot, which are coincident with the numerical results. It is shown that the
scalar cloud with larger mass oscillates faster in the rotating hairy
Horndeski black hole.

##  IV Conclusions

Neutral black holes in general relativity are known to be characterized by two
hairs: the mass and spinning parameter. However, theories of gravity beyond GR
usually allow additional hairs, which in general modify various properties of
the spacetime. In this paper, with the use of matrix method, we investigated
the eigenfrequencies associated with the perturbation of a massive scalar
field around a rotating hairy black hole in Horndeski theory.

Under the perturbation of massless scalar field, we found that comparing to
the Kerr black hole with the same fixed spin parameter, the Horndeski hair
could make the perturbation mode $\ell=m=0$ oscillate slower and decay slower,
corresponding to smaller real part and larger imaginary part of the QNFs. The
effect of the spin parameter on the QNFs could depend on the fixed value of
$h/M$. For small $\mid h/M\mid$, the real part of QNFs first increases and
then decreases with the increase of $a/M$, which is similar to that found in
Kerr black hole. While for large enough $\mid h/M\mid$, it decreases
monotonically, which is understandable because the strong suppression effect
of $\mid h/M\mid$ on the real part could balance the enhancement effect of
$a/M$ and even dominates when the black hole closes to the extreme. Meanwhile,
for all the selected values of $h/M$, the imaginary part first increases with
the increasing of $a/M$ but then decreases slightly as the black hole closes
to the extreme. In addition, QNFs with $\ell=1$ and $m=-1,0,1$ were also
calculated. The Horndeski hair also has a significant influence on those QNFs,
though the rules may be different from those for the mode $\ell=m=0$.
Nevertheless, all the QNFs for various modes have finite negative imaginary
part, which in some sense suggests that the rotating hairy Horndeski black
hole is stable under the perturbation of a massless scalar field.

Then we turned on the mass of the scalar field perturbation, and mainly
checked the eigenfrequencies of the quasibound state. For the scalar mass $\mu
M=0.4$, we found significant effects of the Horndeski hair and spin parameter
on the eigenfrequencies of the quasibound states with $\ell=0$ and $\ell=1$.
Especially, the stronger Horndeski hair always suppresses both the real and
imaginary parts of the eigenfrequency for all the modes, except for the mode
$\ell=m=1$ with $a/M=0.8$ of which the imaginary part first decreases slightly
and then increases to be positive. This could be associated with the
superradiant instability of the black hole under the perturbed mode
$\ell=m=1$. Afterward, we scanned the parameter space to fix the onset of
superradiant instability, also known as bound state or scalar cloud, which
satisfies $Im(M\omega)=0$ and $Re(M\omega)=m\Omega_{H}$. Our results indicate
that due to the existence of the Horndeski hair, the rotating hairy black hole
with a smaller spin than Kerr black hole could undergo superradiant
instability. Finally, we discussed the effect of mass of the scalar field on
the quasibound state and superradiant instability. Our study in the selected
parameter region, shew that the upper bound of the scalar field mass for the
superradiant instability increases when the value $h/M$ is smaller. This
implies that comparing to the Kerr black hole, the range of the scalar field
mass for superradiant instability of the rotating hairy Horndeski black hole
may be extended.

To conclude, this study provided the preliminary step to analyze the
(in)stability of the rotating hairy black hole in Horndeski theory, which we
believe deserves to extend into the perturbations of fields with higher spin.
Our findings on superradiant instability expand the approximate analytical
results of Jha et al. (2022). Moreover, our results rich the phenomenal
properties introduced by the Horndeski hair, which we hope to give some
insight in the test of the no-hair theorem of black holes as well as the
observational applications of hairy black holes.

###### Acknowledgements.

This work is partly supported by Natural Science Foundation of China under
Grants No. 12375054, Natural Science Foundation of Jiangsu Province under
Grant No.BK20211601 and the Postgraduate Research & Practice Innovation
Program of Jiangsu Province under Grant No. KYCX23_3501.

## References

  * Nollert (1999) H.-P. Nollert, Class. Quant. Grav. 16, R159 (1999).
  * Berti et al. (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), eprint 0905.2975.
  * Konoplya and Zhidenko (2011) R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011), eprint 1102.4014.
  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957).
  * Konoplya (2003) R. A. Konoplya, Phys. Rev. D 68, 024018 (2003), eprint gr-qc/0303052.
  * Cardoso et al. (2004) V. Cardoso, J. P. S. Lemos, and S. Yoshida, Phys. Rev. D 69, 044004 (2004), eprint gr-qc/0309112.
  * Cho (2003) H. T. Cho, Phys. Rev. D 68, 024003 (2003), eprint gr-qc/0303078.
  * Jing (2005) J.-l. Jing, Phys. Rev. D 71, 124006 (2005), eprint gr-qc/0502023.
  * Starobinsky (1973) A. A. Starobinsky, Sov. Phys. JETP 37, 28 (1973).
  * Penrose and Floyd (1971) R. Penrose and R. M. Floyd, Nature 229, 177 (1971).
  * Misner (1972) C. W. Misner, Phys. Rev. Lett. 28, 994 (1972).
  * Brito et al. (2015) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015), eprint 1501.06570.
  * Press and Teukolsky (1972) W. H. Press and S. A. Teukolsky, Nature 238, 211 (1972).
  * Press and Teukolsky (1973) W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973).
  * Teukolsky and Press (1974) S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974).
  * Damour et al. (1976) T. Damour, N. Deruelle, and R. Ruffini, Lett. Nuovo Cim. 15, 257 (1976).
  * Zouros and Eardley (1979) T. J. M. Zouros and D. M. Eardley, Annals Phys. 118, 139 (1979).
  * Detweiler (1980) S. L. Detweiler, Phys. Rev. D 22, 2323 (1980).
  * Li (2012) R. Li, Phys. Lett. B 714, 337 (2012), eprint 1205.3929.
  * Dolan (2007) S. R. Dolan, Phys. Rev. D 76, 084001 (2007), eprint 0705.2880.
  * Cardoso and Yoshida (2005) V. Cardoso and S. Yoshida, JHEP 07, 009 (2005), eprint hep-th/0502206.
  * Hod and Hod (2010) S. Hod and O. Hod, Phys. Rev. D 81, 061502 (2010), eprint 0910.0734.
  * Damour and Ruffini (1976) T. Damour and R. Ruffini, Phys. Rev. D 14, 332 (1976).
  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), eprint 1602.03837.
  * Abbott et al. (2019) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), eprint 1811.12907.
  * Akiyama et al. (2019) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 875, L2 (2019), eprint 1906.11239.
  * Akiyama et al. (2022) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022).
  * Riess et al. (1998) A. G. Riess et al. (Supernova Search Team), Astron. J. 116, 1009 (1998), eprint astro-ph/9805201.
  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * Carroll (2001) S. M. Carroll, Living Rev. Rel. 4, 1 (2001), eprint astro-ph/0004075.
  * Clifton et al. (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rept. 513, 1 (2012), eprint 1106.2476.
  * Bakopoulos (2020) A. Bakopoulos, Ph.D. thesis, Ioannina U. (2020), eprint 2010.13189.
  * Corelli (2020) F. Corelli, Master’s thesis, Rome U. (2020), eprint 2112.12048.
  * Horndeski (1974) G. W. Horndeski, Int. J. Theor. Phys. 10, 363 (1974).
  * Kobayashi (2019) T. Kobayashi, Rept. Prog. Phys. 82, 086901 (2019), eprint 1901.07183.
  * Deffayet et al. (2011) C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, Phys. Rev. D 84, 064039 (2011), eprint 1103.3260.
  * Kobayashi et al. (2011) T. Kobayashi, M. Yamaguchi, and J. Yokoyama, Prog. Theor. Phys. 126, 511 (2011), eprint 1105.5723.
  * Bellini et al. (2016) E. Bellini, A. J. Cuesta, R. Jimenez, and L. Verde, JCAP 02, 053 (2016), [Erratum: JCAP 06, E01 (2016)], eprint 1509.07816.
  * Bhattacharya and Chakraborty (2017) S. Bhattacharya and S. Chakraborty, Phys. Rev. D 95, 044037 (2017), eprint 1607.03693.
  * Kreisch and Komatsu (2018) C. D. Kreisch and E. Komatsu, JCAP 12, 030 (2018), eprint 1712.02710.
  * Hou and Gong (2018) S. Hou and Y. Gong, Eur. Phys. J. C 78, 247 (2018), eprint 1711.05034.
  * Spurio Mancini et al. (2019) A. Spurio Mancini, F. Köhlinger, B. Joachimi, V. Pettorino, B. M. Schäfer, R. Reischke, E. van Uitert, S. Brieden, M. Archidiacono, and J. Lesgourgues, Mon. Not. Roy. Astron. Soc. 490, 2155 (2019), eprint 1901.03686.
  * Allahyari et al. (2020) A. Allahyari, M. A. Gorji, and S. Mukohyama, JCAP 05, 013 (2020), [Erratum: JCAP 05, E02 (2021)], eprint 2002.11932.
  * Feng et al. (2015) X.-H. Feng, H.-S. Liu, H. Lü, and C. N. Pope, JHEP 11, 176 (2015), eprint 1509.07142.
  * Kuang and Papantonopoulos (2016) X.-M. Kuang and E. Papantonopoulos, JHEP 08, 161 (2016), eprint 1607.04928.
  * Jiang et al. (2017) W.-J. Jiang, H.-S. Liu, H. Lu, and C. N. Pope, JHEP 07, 084 (2017), eprint 1703.00922.
  * Baggioli and Li (2017) M. Baggioli and W.-J. Li, JHEP 07, 055 (2017), eprint 1705.01766.
  * Feng and Liu (2019) X.-H. Feng and H.-S. Liu, Eur. Phys. J. C 79, 40 (2019), eprint 1811.03303.
  * Wang et al. (2019) X.-J. Wang, H.-S. Liu, and W.-J. Li, Eur. Phys. J. C 79, 932 (2019), eprint 1909.00224.
  * Zhang et al. (2023) D. Zhang, G. Fu, X.-J. Wang, Q. Pan, and J.-P. Wu, Eur. Phys. J. C 83, 316 (2023), eprint 2211.07074.
  * Bravo-Gaete and Santos (2022) M. Bravo-Gaete and F. F. Santos, Eur. Phys. J. C 82, 101 (2022), eprint 2010.10942.
  * Bravo-Gaete et al. (2022) M. Bravo-Gaete, F. F. Santos, and H. Boschi-Filho, Phys. Rev. D 106, 066010 (2022), eprint 2201.07961.
  * Hui and Nicolis (2013) L. Hui and A. Nicolis, Phys. Rev. Lett. 110, 241104 (2013), eprint 1202.1296.
  * Sotiriou and Zhou (2014a) T. P. Sotiriou and S.-Y. Zhou, Phys. Rev. Lett. 112, 251102 (2014a), eprint 1312.3622.
  * Sotiriou and Zhou (2014b) T. P. Sotiriou and S.-Y. Zhou, Phys. Rev. D 90, 124063 (2014b), eprint 1408.1698.
  * Babichev et al. (2017) E. Babichev, C. Charmousis, and A. Lehébel, JCAP 04, 027 (2017), eprint 1702.01938.
  * Rinaldi (2012) M. Rinaldi, Phys. Rev. D 86, 084048 (2012), eprint 1208.0103.
  * Cisterna and Erices (2014) A. Cisterna and C. Erices, Phys. Rev. D 89, 084038 (2014), eprint 1401.4479.
  * Babichev and Charmousis (2014) E. Babichev and C. Charmousis, JHEP 08, 106 (2014), eprint 1312.3204.
  * Babichev et al. (2018) E. Babichev, C. Charmousis, G. Esposito-Farèse, and A. Lehébel, Phys. Rev. Lett. 120, 241101 (2018), eprint 1712.04398.
  * Ben Achour and Liu (2019) J. Ben Achour and H. Liu, Phys. Rev. D 99, 064042 (2019), eprint 1811.05369.
  * Khoury et al. (2020) J. Khoury, M. Trodden, and S. S. C. Wong, JCAP 11, 044 (2020), eprint 2007.01320.
  * Bergliaffa et al. (2021) S. E. P. Bergliaffa, R. Maier, and N. d. O. Silvano (2021), eprint 2107.07839.
  * Yang et al. (2023) Z.-H. Yang, Y.-H. Lei, X.-M. Kuang, and J.-P. Wu (2023), eprint 2309.03565.
  * Kumar et al. (2022) J. Kumar, S. U. Islam, and S. G. Ghosh, Eur. Phys. J. C 82, 443 (2022), eprint 2109.04450.
  * Wang et al. (2023) X.-J. Wang, X.-M. Kuang, Y. Meng, B. Wang, and J.-P. Wu, Phys. Rev. D 107, 124052 (2023), eprint 2304.10015.
  * Azreg-Aïnou (2014a) M. Azreg-Aïnou, Phys. Rev. D 90, 064041 (2014a), eprint 1405.2569.
  * Azreg-Aïnou (2014b) M. Azreg-Aïnou, Eur. Phys. J. C 74, 2865 (2014b), eprint 1401.4292.
  * Walia et al. (2022) R. K. Walia, S. D. Maharaj, and S. G. Ghosh, Eur. Phys. J. C 82, 547 (2022), eprint 2109.08055.
  * Afrin and Ghosh (2022) M. Afrin and S. G. Ghosh, Astrophys. J. 932, 51 (2022), eprint 2110.05258.
  * Jha et al. (2022) S. K. Jha, M. Khodadi, A. Rahaman, and A. Sheykhi (2022), eprint 2212.13051.
  * Lin and Qian (2017) K. Lin and W.-L. Qian, Class. Quant. Grav. 34, 095004 (2017), eprint 1610.08135.
  * Lin et al. (2017) K. Lin, W.-L. Qian, A. B. Pavan, and E. Abdalla, Mod. Phys. Lett. A 32, 1750134 (2017), eprint 1703.06439.
  * Lin and Qian (2019) K. Lin and W.-L. Qian, Chin. Phys. C 43, 035105 (2019), eprint 1902.08352.
  * Siqueira and Richartz (2022) P. H. C. Siqueira and M. Richartz, Phys. Rev. D 106, 024046 (2022), eprint 2205.00556.
  * Hod (2012) S. Hod, Phys. Rev. D 86, 104026 (2012), [Erratum: Phys.Rev.D 86, 129902 (2012)], eprint 1211.3202.
  * Herdeiro and Radu (2014) C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. 112, 221101 (2014), eprint 1403.2757.
  * Benone et al. (2014) C. L. Benone, L. C. B. Crispino, C. Herdeiro, and E. Radu, Phys. Rev. D 90, 104024 (2014), eprint 1409.1593.