# Compact Binary Formation in Open Star Clusters III:
Probability of Binary Black Holes Hidden Inside of Gaia Black Holes

Ataru Tanikawa1 0000-0002-8461-5517 Long Wang2,3 0000-0001-8713-0366 Michiko
S. Fujii4 0000-0002-6465-2978 1Center for Information Science, Fukui
Prefectural University, 4-1-1 Matsuoka Kenjojima, Eiheiji-cho, Fukui 910-1195,
Japan  2School of Physics and Astronomy, Sun Yat-sen University, Daxue Road,
Zhuhai, 519082, China  3CSST Science Center for the Guangdong-Hong Kong-Macau
Greater Bay Area, Zhuhai, 519082, China  4Department of Astronomy, Graduate
School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo
113-0033, Japan  Corresponding author<EMAIL_ADDRESS>
###### Abstract

Gaia mission and its follow-up observations have discovered binaries
containing single BHs and visible stars without mass transfer, so-called Gaia
BHs. One important question is if Gaia BHs have binary BHs (BBHs), hereafter
Gaia BBHs, instead of single BHs. We have investigated how efficiently Gaia
BBHs are formed in open star clusters, one of the promising formation sites of
Gaia BHs, by means of gravitational $N$-body simulations. Limiting Gaia BHs’
periods to $10^{2}$-$10^{4}$ days, we have found that there are no Gaia BBHs
in the solar-metallicity environments, while the formation efficiency of Gaia
BBHs is not small ($\sim 10^{-6}M_{\odot}^{-1}$ or $\sim 10$ % of Gaia BHs) in
subsolar-metallicity environments. However, the probability of Gaia BBHs
hidden in Gaia BHs is only $\sim 1$ %. This is because most of the BBHs merge
within $1$ Gyr through gravitational wave radiation. Note that the ages of
discovered Gaia BHs are more than $1$ Gyr. If we extend Gaia BHs’ periods to
$10^{4}$-$10^{5}$ days, the probability becomes higher to $\sim 10$ %. In this
case, a large fraction of BBHs can have enough wide orbits not to merge within
the Hubble time. The probability would not be high for Gaia BHs already
discovered and to be discovered in the near future. Nevertheless, we have
shown the BH/BBH mass, visible star mass, and eccentricity distributions of
Gaia BHs and Gaia BBHs, which will be helpful for follow-up observations to
discover Gaia BBHs. A Gaia BH would be more likely to be a Gaia BBH if it has
younger age, longer period, lower-mass companion, more circular orbit, lower
metallicity, and more massive BH. Our results have implied that Gaia BH3 is
unlikely to be a Gaia BBH.

##  1 Introduction

Massive stars leave behind stellar-mass black holes (BHs) with $\sim
10$-$100M_{\odot}$ at the end of their lives. BHs are a laboratory for testing
strong gravity, a site of high-energy phenomena, and a probe to construct the
theory of massive single and binary star evolution. This motivates BH hunters
to search for various types of BHs in various ways: BH binaries with single
BHs and single visible stars by X-ray observations (see Casares et al., 2017,
for review), spectroscopic observations (Giesers et al., 2018, 2019; Shenar et
al., 2022), and astrometric observations (see El-Badry, 2024b, for review),
single BHs by gravitational microlensing (Lam et al., 2022; Sahu et al., 2022;
Howil et al., 2024), and binary BHs (BBHs) with two BHs by gravitational wave
observations (Abbott et al., 2023a).

Gaia mission (Gaia Collaboration et al., 2023b, a) offers chances to discover
BH binaries without mass transfer, and its follow-up observations have found
three BH binaries: Gaia BH1 (El-Badry et al., 2023b; Chakrabarti et al.,
2023), Gaia BH2 (Tanikawa et al., 2023; El-Badry et al., 2023a), and Gaia BH3
(Gaia Collaboration et al., 2024). Hereafter, we call these BH binaries “Gaia
BHs”. The formation process of Gaia BHs is an open question. Gaia BHs have
long periods ($100$-$4000$ days) and substantial eccentricities ($0.5$-$0.8$).
They cannot be formed through isolated binary evolution without 10 times more
efficient common envelope ejection than expected theoretically (Chawla et al.
2022; El-Badry et al. 2023b, a; Shikauchi et al. 2023; but see Kotko et al.
2024)111Actually, this is also true for neutron stars and white dwarf binaries
discovered from Gaia database (Ganguly et al., 2023; Geier et al., 2023; Zheng
et al., 2023; Yamaguchi et al., 2024; El-Badry et al., 2024a, b).. However,
Hirai & Mandel (2022) have claimed that common envelope ejection cannot be so
efficient for Gaia BH progenitors. Another possibility is dynamical formation
of Gaia BHs in dense star clusters. Since Gaia BH1 and BH2 are the Galactic
disk components (El-Badry et al., 2023b, a), they can be formed in open
clusters (Shikauchi et al., 2020; Rastello et al., 2023; Tanikawa et al.,
2024a; Di Carlo et al., 2024). Gaia BH3 is located in the ED-2 stream
(Balbinot et al., 2024), and thus it can be formed in a small globular cluster
(Kremer et al., 2018; Marín Pina et al., 2024). Note that Gaia BH3 can be
formed through isolated binary evolution because of its low metallicity (El-
Badry, 2024a; Iorio et al., 2024). Another noticeable scenario is Gaia BH
formation in triple star systems (El-Badry et al., 2023b, a; Generozov &
Perets, 2024). Gaia BHs have received substantial attention, because they are
good samples for studying BH natal kicks (Vigna-Gómez et al., 2024), and
interactions in wide binaries with compact objects (Cappelluti et al., 2024).

One important question is if any Gaia BH actually contains a BBH instead of a
single BH (Hayashi et al., 2022, 2023). Hereafter, we call it “Gaia BBH” (see
Figure 1). The presence of a Gaia BBH is testable by observing modulated
motion of a visible star (Hayashi et al., 2020; Hayashi & Suto, 2020; Liu et
al., 2022). Nagarajan et al. (2024) have already shown no convincing evidence
of a BBH in Gaia BH1. Since Gaia BH3 contains a massive BH ($\sim
33M_{\odot}$) (Gaia Collaboration et al., 2024), it may be expected to be a
Gaia BBH, however there is no constraint. No Gaia BBH has been discovered.
Nevertheless, if a large number of Gaia BHs will be found in the future, they
should contain a handful of Gaia BBHs.

Figure 1: Cartoon illustration of Gaia BH and Gaia BBH.

The discovery of only one Gaia BBH will have great impacts. If it contains a
BBH merging within the Hubble time, it will constrain the formation mechanism
of BBHs observed by gravitational wave observatories (Abbott et al., 2023a),
whose origins have been under debate (Abbott et al., 2023b). Gravitational
wave observation can reveal BH masses and spins of a BBH. On the other hand, a
Gaia BBH will provide its position, velocity, and metallicity via its visible
companion star. Its position and velocity will be helpful to make clear what
components (the Galactic disk, halo, bulge, and center) the Gaia BBH belong
to. Note that the sky localization and distance obtained from gravitational
wave observation are not helpful for elucidating the origin of the BBH, since
there are a huge number of galaxies in the possible position of the BBH. The
metallicity of the Gaia BBH will indicate what metallicity is important for
BBHs observed.

In this paper, we study the probability that Gaia BHs contain BBHs, or are
actually Gaia BBHs, and the difference of properties between Gaia BHs and Gaia
BBHs, based on a working hypothesis that Gaia BHs are dominantly formed in
open clusters. Since Gaia BH1 and BH2 are the Galactic disk components (El-
Badry et al., 2023b, a), they can be considered to be born in open clusters.
Although Gaia BH3 is a Galactic halo star (Gaia Collaboration et al., 2024),
it may be a member of the ED-2 stream, a disrupted star cluster with $2\times
10^{3}$–$4.2\times 10^{4}M_{\odot}$ (Balbinot et al., 2024), which is in the
mass range of open clusters. Tanikawa et al. (2024b) have claimed that neutron
star (NS) binaries discovered from the Gaia Database (so-called Gaia NSs)
cannot be formed through dynamical capture in open clusters, since open
clusters form more Gaia BHs than Gaia NSs, not in accordance with $\gtrsim 20$
Gaia NSs discovered (Geier et al., 2023; Zheng et al., 2023; El-Badry et al.,
2024a, b). Nevertheless, it cannot rule out the possibility that Gaia BHs are
dominantly formed in open clusters.

The structure of this paper is as follows. In section 2, we present our
numerical method to form Gaia BHs and Gaia BBHs in open clusters. In section
3, we show our results. In section 4, we make a conclusion and discussion.

##  2 Method

We follow the dynamical evolution of open clusters by means of gravitational
$N$-body simulations in the same way as Tanikawa et al. (2024b). Here, we
briefly describe our methods, especially related to the purpose of this paper.
We employ an $N$-body code PETAR (Wang et al., 2020a), highly optimized by
FDPS (Iwasawa et al., 2016, 2020). Binary and close-encounter orbits are
solved with high accuracy thanks to SDAR (Wang et al., 2020b). The model of
the Galaxy is constructed by GALPY (Bovy, 2015). The PETAR code is equipped
with the BSE code (Hurley et al., 2000, 2002; Banerjee et al., 2020) to treat
single and binary star evolutions. For our supernova model, we adopt the rapid
model (Fryer et al., 2012) with the correction of the moderate pair
instability and pulsational pair-instability supernovae (Belczynski et al.,
2020a). NSs and BHs receive natal kicks at their births, and their kick
velocities are distributed by a single Maxwellian with $265$ ${\rm
km~{}s}^{-1}$ (Hobbs et al., 2005). The kick velocities are decreased by
$1-f_{\rm fb}$ where $f_{\rm fb}$ is the fraction of the fallback mass (Fryer
et al., 2012). The BSE code handles common envelope evolution with the
$\alpha$ formalism (Webbink, 1984). In our simulations, we choose $\alpha_{\rm
CE}=1$ and $\lambda_{\rm CE}$ of Claeys et al. (2014). Our choice of
$\alpha_{\rm CE}=1$ and $\lambda_{\rm CE}$ prohibits Gaia BH formation from
primordial binaries.

Table 1: Summary of cluster models. Name | $Z$ | $M$ | $\rho$ | $f_{\rm b,<5M_{\odot}}$
---|---|---|---|---
|  | [$M_{\odot}$] | [$M_{\odot}~{}{\rm pc}^{-3}$] |
$Z=0.02$ | $0.02$ | $1000$ | $20$ |  $20$ %
$Z=0.01$ | $0.01$ | $1000$ | $20$ |  $20$ %
$Z=0.005$ | $0.005$ | $1000$ | $20$ |  $20$ %
$Z=0.002$ | $0.002$ | $1000$ | $20$ |  $20$ %
$Z=0.0002$ | $0.0002$ | $1000$ | $20$ |  $20$ %

The initial conditions of our cluster models are as follows (see also Table
1). Each cluster mass ($M$) is $10^{3}M_{\odot}$, and each cluster mass
density within its half-mass radius ($\rho$) is $20M_{\odot}~{}{\rm pc}^{-3}$
at the initial time. Each cluster circularly orbits at the distance of $8$ kpc
from the Galactic center. Initial binary fractions vary depending on primary
star masses ($m_{1}$). The initial binary fractions for $m_{1}\geq 5M_{\odot}$
($f_{\rm b,\geq 5M_{\odot}}$) and $m_{1}<5M_{\odot}$ ($f_{\rm b,<5M_{\odot}}$)
are $100$ and $20$ %, respectively. The initial mass function (IMF) of single
stars and primary stars in binaries is Kroupa’s IMF (Kroupa, 2001). Binary
parameters, such as secondary star masses, binary periods, and eccentricities,
follow initial conditions of Sana et al. (2012) for $m_{1}\geq 5M_{\odot}$,
and initial conditions of Kroupa (1995a, b) with modifications of Belloni et
al. (2017) for $m_{1}<5M_{\odot}$222Note that we adopt $f_{\rm
b,<5M_{\odot}}=20$ % despite that Kroupa (1995a, b) have adopted $f_{\rm
b,<5M_{\odot}}=100$ %. If $f_{\rm b,<5M_{\odot}}>20$ % in reality, we
underestimate the formation efficiency of Gaia BHs (Tanikawa et al., 2024b).
This may be the case for the formation efficiency of Gaia BBHs. However, we do
not expect that $f_{\rm b,<5M_{\odot}}$ significantly affects the probability
of Gaia BBHs hidden in Gaia BHs, since the formation efficiencies of Gaia BHs
and Gaia BBHs increase with $f_{\rm b,<5M_{\odot}}$ increasing.. We choose
various metallicities: $Z=0.02,0.01,0.005,0.002$, and $0.0002$. For each
metallicity, we generate $10^{4}$ cluster models. In other words, we
investigate clusters with the total mass of $5\times 10^{7}M_{\odot}$. We
construct these cluster models with MCLUSTER (Küpper et al., 2011). We finish
cluster simulations at $1$ Gyr.

Tanikawa et al. (2024b) have treated other cluster models with different
cluster masses ($M=200$, $500$, and $2000$ $M_{\odot}$), cluster mass
densities ($\rho=2$ and $200M_{\odot}~{}{\rm pc}^{-3}$), and initial binary
fractions for $m_{1}<5M_{\odot}$ ($f_{\rm b,<5M_{\odot}}=0$ and $50$ %) from
the above cluster models with $M=1000M_{\odot}$, $\rho=20M_{\odot}~{}{\rm
pc}^{-3}$, and $f_{\rm b,<5M_{\odot}}=20$ %. However, these other cluster
models have high metallicities ($Z=0.02$), and form no Gaia BBHs because of
their high metallicity. Even if we include these other cluster models in this
paper, it does not affect our conclusions. Thus, in this paper, we decide that
it is sufficient to see the $Z=0.02$ cluster model with $M=10^{3}M_{\odot}$,
$\rho=20M_{\odot}~{}{\rm pc}^{-3}$, and $f_{\rm b,<5M_{\odot}}=20$ % for
considering high-metallicity clusters.

We pick up Gaia BHs and Gaia BBHs from our cluster models according to the
following conditions.

  1. 1.

$10^{2}\leq P/{\rm day}\leq 10^{4}$, where $P$ is an orbital period of a Gaia
BH or Gaia BBH. Note that the orbital period of a Gaia BBH can be defined as
an outer binary period, since astrometric observation can determine only an
orbital period between BBH and a visible star.

  2. 2.

$m_{2}\leq 1.1M_{\odot}$, where $m_{2}$ is a visible star mass of a Gaia BH or
Gaia BBH. This constraint is motivated by the visible star masses of Gaia BH1
($0.93M_{\odot}$), BH2 ($1.07M_{\odot}$), and BH3 ($0.76M_{\odot}$).

  3. 3.

Outside of clusters. None of Gaia BHs are located in clusters.

In order to understand the properties of Gaia BHs and Gaia BBHs, we also
investigate BH binaries and triple systems having inner BBHs both of which
satisfy only conditions 2 and 3. In other words, such BH binaries have any
binary periods, and such triple systems have any outer binary periods. We call
the former “BH binaries” and the latter “triples with BBHs”. Note that BH
binaries and triples with BBHs include Gaia BHs and Gaia BBHs, respectively.

We summarize the terms used in this paper. Since Gaia BBHs are triple-star
systems, they consist of inner and outer binaries. When we mention orbits,
periods, and eccentricities of Gaia BBHs, they mean outer binary orbits,
periods, and eccentricities, respectively, if unstated. This is because they
correspond to orbits, periods, and eccentricities of Gaia BHs. When we focus
on BBH orbits, we call them “inner BBH orbits”. When we refer to BHs in Gaia
BHs and BBHs in Gaia BBHs collectively, we call them “unseen objects” or
“BH/BBH”.

##  3 Results

We make clear the overall features of Gaia BBHs formed in our cluster models.
We find $59$ Gaia BBHs at $1$ Gyr in our cluster models. The numbers of Gaia
BBHs are $0$, $1$, $11$, $24$, and $23$ for the $Z=0.02$, $Z=0.01$, $Z=0.005$,
$0.002$, and $0.0002$ models, respectively. All these BBHs are formed from
primordial binaries, not through dynamical capture, and capture visible stars
to form Gaia BBHs. We do not obtain any quadruple systems observed similarly
to Gaia BBHs: an outer binary consisting of a BBH and another binary with at
least one visible star.

Figure 2: Formation efficiencies of Gaia BBHs at $1$ Gyr as a function of
cluster metallicity. Filled points with the dotted curve indicate the
formation efficiency of all Gaia BBHs, while filled points with dashed and
solid curves show the formation efficiencies of Gaia BBHs with gravitational
wave merger timescales with $>1$ and $>10$ Gyr, respectively. For reference,
the formation efficiency of Gaia BHs is also indicated by open points with the
solid curve. For visibility, points with the dotted and solid curves are
shifted by $0.05$ dex rightward and leftward, respectively. We obtain error
bars, assuming that Gaia BHs and Gaia BBHs are formed according to a Poisson
process. Points with down arrows indicate no Gaia BBHs.

Figure 2 shows the formation efficiencies of Gaia BHs and Gaia BBHs as a
function of cluster metallicity. Gaia BHs are formed at an efficiency of $\sim
10^{-6}$–$10^{-5}M_{\odot}^{-1}$ as shown in Tanikawa et al. (2024b). The
formation efficiencies of all Gaia BBHs strongly depend on cluster
metallicity. They are $\sim 10^{-6}M_{\odot}^{-1}$ in low-metallicity
environments ($Z\leq 0.005$), while they are $10^{-7}M_{\odot}$ for $Z=0.01$
and $0$ for $Z=0.02$. This metallicity dependence results from the formation
efficiencies of merging BBHs in isolated binaries; a substantial amount of
BBHs are formed in low-metallicity environments, while little BBHs are formed
in the solar-metallicity environments (e.g. Santoliquido et al., 2020). Note
that BBHs in Gaia BBHs are formed from primordial binaries (or effectively
isolated binaries even in open clusters) as described above. In the solar-
metallicity environments, primordial binaries in close orbits merge when their
primary stars evolve to Hertzsprung gap stars (e.g. Tanikawa et al., 2022).
Thus, they cannot form merging BBHs.

A large fraction of these Gaia BBHs contain BBHs merging due to gravitational
wave radiation within the Hubble time. The merger timescale distribution is
logarithmically flat in a similar way to merger timescales of double compact
binaries (e.g. Totani et al., 2008). Thus, not all of them can be observed as
Gaia BBHs. After these BBHs merge, gravitational wave recoil kicks can
separate visible stars from BBH merger remnants, and thus the systems are
disrupted. Even if the systems survive, they can be observed as Gaia BHs. Most
of BBHs in Gaia BBHs merge within $1$ Gyr. The formation efficiencies of Gaia
BBHs with merger timescale of $1$ Gyr are only $3\times 10^{-7}M_{\odot}^{-1}$
for $Z\leq 0.005$. For $Z>0.005$, the formation efficiencies are much smaller
than that.

The formation efficiencies of Gaia BBHs with merger timescale of $10$ Gyr are
just around $10^{-7}M_{\odot}^{-1}$ for $Z\leq 0.01$. Eventually, the
probability of Gaia BBHs hidden in Gaia BHs is only $\sim 1$ %, since the ages
of Gaia BH1, BH2, and BH3 are thought to be a few Gyr (El-Badry et al., 2023b,
a; Gaia Collaboration et al., 2024, respectively). This probability may be
reduced further. Although we consider only inner BBH orbits to estimate the
merger timescales of BBHs in Gaia BBHs with a gravitational wave radiation
formula (Peters & Mathews, 1963; Peters, 1964), Kozai-Lidov mechanism (Kozai,
1962; Lidov, 1962) may accelerate the merger process (Seto, 2013; Antonini &
Rasio, 2016; Antonini et al., 2017; Silsbee & Tremaine, 2017; Rodriguez &
Antonini, 2018; Toonen et al., 2018; Martinez et al., 2020; Vigna-Gómez et
al., 2021; Arca Sedda et al., 2021; Martinez et al., 2022; Trani et al., 2022;
Liu et al., 2023).

Figure 3: Formation efficiencies of BH binaries and triples with BBHs for each
dex of orbital periods in the $Z=0.0002$ model. Open and filled circles
indicate BH binaries and triples with BBHs, respectively. Solid, dashed, and
dotted curves show BBHs with merger timescales of $>10$ Gyr, $>1$ Gyr, and
all. The shaded region corresponds to orbital periods of Gaia BHs and Gaia
BBHs. The meaning of error bars and down arrows are the same as in Figure 2.

A large fraction of Gaia BBHs have merger timescale smaller than the Hubble
time as shown in Figure 2. This is because Gaia BBHs with merger timescale
larger than the Hubble time tend to be dynamically and secularly unstable.
Even once such a Gaia BBH is formed, it is disrupted soon. This argument is
supported by Figure 3. The formation efficiencies of BH binaries increases
with their periods increasing, while the formation efficiencies of triples
with BBHs with merger timescale more than the Hubble time increases more
rapidly with their periods increasing. Thus, the probability of BBHs hidden in
BH binaries increases to about $10$ % when BH binaries’ periods are more than
$10^{4}$ days. Their orbits are stable even if their inner BBHs have enough
wide orbits not to merge within the Hubble time. However, we infer that Gaia
mission cannot find such triples with BBHs as well as BH binaries. Gaia
mission duration is about $10$ years or less than $4000$ days. It is difficult
to determine orbital solutions of binaries with $>2$ times larger periods than
an operation duration (e.g. Lucy, 2014; O’Neil et al., 2019)

Unfortunately, the probability of BBHs hidden in Gaia BHs should be small.
Nevertheless, when $100$ Gaia BHs would be discovered in the future, they may
contain about $1$ Gaia BBH. In preparation for that time, it would not be a
waste to make comparison between Gaia BHs and Gaia BBHs formed in our cluster
models. Hereafter, we show the properties of Gaia BHs and Gaia BBHs in our
simulations. In particular, we focus on the $Z=0.0002$ model. Gaia
Collaboration et al. (2024) suspect that Gaia BH3, which has a low-metallicity
companion ($[{\rm Fe/H}]=-2.76$), might be actually a Gaia BBH. We show
properties of Gaia BBHs formed in other $Z$ cluster models in Appendix A.

Figure 4: Cumulative distributions of BH and BBH masses in Gaia BHs and Gaia
BBHs, respectively, for the $Z=0.0002$ model. Grey and black curves indicate
Gaia BHs and Gaia BBHs, respectively. Solid, dashed and dotted curves show
Gaia BBHs with merger timescales of $>10$ Gyr, $>1$ Gyr, and all,
respectively. The dotted grey curve includes all Gaia BHs, while the solid
grey curve excludes Gaia BHs whose BHs are BBH merger remnants. The shaded
region with hatch pattern indicates the parameter range of Gaia BH3. Figure 5:
Cumulative distributions of visible star masses in Gaia BHs and Gaia BBHs for
the $Z=0.0002$ model. Curve types are the same as in Figure 4. Figure 6:
Cumulative distributions of eccentricities of Gaia BHs and Gaia BBHs for the
$Z=0.0002$ model. Curve types are the same as in Figure 4.

Figure 4 shows the cumulative distributions of BH/BBH masses of Gaia BHs and
Gaia BBHs. Before discussing the difference between Gaia BHs and Gaia BBHs, we
should make caveats on Gaia BHs with $>50M_{\odot}$ BHs. Although we obtain
them in our cluster simulations, they are unlikely to be formed in actual open
clusters for the following reason. These BHs are formed through BBH mergers.
In reality, such BBH mergers generate gravitational wave recoil kicks and
eject BBH merger remnants themselves from their host open clusters. Thus, they
should not capture visible stars or form Gaia BHs. However, the PETAR code has
not yet implemented gravitational wave recoil kicks. Thus, BBH merger remnants
are not ejected from the host open clusters and form Gaia BHs in our cluster
simulations. In the $Z=0.0002$ model, $14$ of $89$ Gaia BHs have BBH merger
remnants. This does not much affect the formation efficiency of Gaia BHs. On
the other hand, this changes the shape of the cumulative distribution of BH
masses (see the solid and dotted grey curves in Figure 4). This is because all
Gaia BHs with BH masses of $>50M_{\odot}$ have BBH merger remnants. Although a
part of BBH merger remnants might not be ejected from actual open clusters, we
regard the solid grey curves as correct in Figures 4, 5, and 6. This affects
Gaia BHs in the $Z=0.0002$ model most largely among our cluster models. The
fractions of such Gaia BHs decrease with metallicity increasing: $13$, $1.1$,
$0$, and $0$ % for $Z=0.002$, $0.005$, $0.01$, and $0.02$ models,
respectively. Moreover, none of Gaia BBHs contains BBH merger remnants.

Table 2: P-values of the Kolmogorov-Smirnov (K-S) test for BBH mass, visible star mass, and eccentricity distributions of Gaia BBHs with all merger timescales and $>10$ Gyr merger timescales (All-$10$), and with $>1$ and $>10$ Gyr merger timescales ($1$-$10$). | All-$10$ |  $1$-$10$
---|---|---
K-S P-values of BBH mass | $0.45$ | $1.0$
K-S P-values of visible star mass | $0.73$ | $1.0$
K-S P-values of eccentricity | $0.27$ | $1.0$

We first investigate differences among Gaia BBHs with different merger
timescales. For this purpose, we perform the Kolmogorov-Smirnov (K-S) test
between Gaia BBHs with all merger timescales and $>10$ Gyr merger timescales,
and between Gaia BBHs with $>1$ and $>10$ Gyr merger timescales as seen in
Table 2. We do not find any significant difference in BBH masses, visible star
masses, and eccentricities, although this reason may be why the number of Gaia
BBHs with $>10$ Gyr merger timescales is too small. Note that the numbers of
all Gaia BBHs with all merger timescales and $>1$ Gyr merger timescales, and
$>10$ Gyr merger timescales are $21$, $3$, and $1$, respectively.
Nevertheless, for now, we assume that Gaia BBHs with $>10$ Gyr merger
timescales have similar BBH masses, visible star masses, and eccentricities to
Gaia BBHs with all merger timescales.

Figure 4 shows that BBH masses in Gaia BBHs are 2 times larger than BH masses
in Gaia BHs. This is simply because each Gaia BBH contains two BHs, while each
Gaia BH contain only one BH. We have to remark that the maximum mass of single
BHs ($\sim 50M_{\odot}$) depends on the pair instability and pulsational pair
instability supernovae model we choose. The discovery of a gravitational wave
event GW190521 (Abbott et al., 2020a, b) reveals that there is a large
uncertainty in the lower edge of the pair instability mass gap (e.g.
Takahashi, 2018; Farmer et al., 2020; Kinugawa et al., 2021; Vink et al.,
2021; Tanikawa et al., 2021; Costa et al., 2022; Kawashimo et al., 2023). If a
Gaia BH with BH mass of $>50M_{\odot}$ will be discovered, it may be possibly
an actual Gaia BBH as shown in Figure 4. However, it strongly depends on pair
instability and pulsational pair instability supernovae models. Thus, in order
to confirm that it is a Gaia BBH, we need its radial velocity modulation as
Nagarajan et al. (2024) search for.

We draw the cumulative distributions of visible star masses in Gaia BHs and
Gaia BBHs in Figure 5. Their distributions are quite different. For example,
the median value of visible star masses in Gaia BBHs ($\sim 0.2M_{\odot}$) is
much smaller than the median value of visible star masses in Gaia BHs ($\sim
0.5M_{\odot}$). This can be also observed in the other metallicities (see
Figure 9). This is possibly due to Kozai-Lidov mechanism (Kozai, 1962; Lidov,
1962). Because Kozai-Lidov mechanism timescale is inversely proportional to
tertiary star masses (here, visible star masses), BBHs in Gaia BBHs with more
massive visible stars may merge at an early time. However, note that the
Kozai-Lidov mechanism may become stronger in our cluster simulations than in
reality, since the PETAR code does not take into account the relativistic
precession of BBHs (e.g. Liu et al., 2023).

Figure 6 shows the cumulative distributions of eccentricities of Gaia BHs and
Gaia BBHs. Gaia BBHs have more circular orbits than Gaia BHs. This is because
Gaia BBHs become more dynamically unstable and have shorter Kozai-Lidov
timescales if they have higher eccentricities.

We compare the BH mass, visible star mass, and eccentricity of Gaia BH3 with
those of Gaia BHs and Gaia BBHs obtained from our cluster simulations. For
this comparison, we adopt the $Z=0.0002$ model, the lowest-metallicity model
among our cluster models, since the metallicity of Gaia BH3 is $[{\rm
Fe/H}]=-2.76$. Overall, Gaia BH3 is likely to be Gaia BHs in the $Z=0.0002$
model. Gaia BBHs should be more massive in BH mass, less massive in visible
star mass, and less eccentric than Gaia BH3.

##  4 Conclusion and Discussion

We investigated the probability of Gaia BBHs, or BBHs hidden inside of Gaia
BHs. The probability is low, at most $1$ % for $Z\leq 0.01$. The formation
efficiency of Gaia BBHs is not so small, $10$ % of the formation efficiency of
Gaia BHs, for $Z\leq 0.005$. However, a large fraction of them experience BBH
mergers within the Hubble time. Such Gaia BBHs cannot be discovered, since
Gaia BHs discovered so far have the ages of a few Gyr. The reason why Gaia
BBHs have short merger timescales is the presence of visible stars, or the
tertiary stars of Gaia BBH systems. Although Gaia BBHs with longer merger
timescales can be formed, they should be soon disrupted, because their
configurations become more unstable. Thus, a small number of Gaia BBHs have
longer merger timescale than the Hubble time.

The fraction of Gaia BBHs to Gaia BHs is small, however not zero. Thus, we
compared the properties of Gaia BBHs with those of Gaia BHs. Since each Gaia
BBH contains two BHs, their expected unseen object masses are simply two times
larger than the unseen object masses of Gaia BHs. Gaia BBHs have lower-mass
visible stars and more circular orbits than Gaia BHs. In our simulations, both
Gaia BHs and Gaia BBHs are formed through the dynamical capture of visible
stars. Thus, they should have had similar visible star masses and orbits when
they were formed. The difference between the properties of Gaia BHs and Gaia
BBHs should come from the dependence of the lifetimes of Gaia BBHs on their
visible star masses and eccentricities. Gaia BBHs should be more robust
against dynamical and secular instability if they have lower-mass companions
and more circular orbits. They tend to have lower-mass visible stars and more
circular orbits.

We assess if Gaia BH3 is actually a Gaia BBH. Based on our simulation results,
Gaia BH3 is unlikely to be a Gaia BBH. Its BH mass is too small (see Figure
4), its visible star is too massive (see Figure 5), and its orbit is too
eccentric (see Figure 6) to be a Gaia BBH.

We show a guideline to discover a Gaia BBH from a large amount of Gaia BHs. A
Gaia BH satisfying the following properties could be a Gaia BBH.

  * •

Young Gaia BHs (see Figure 2). They can contain BBHs with short merger
timescales.

  * •

Gaia BHs with long periods, low-mass visible stars, and low eccentricities
(see Figures 3, 5, and 6, respectively). Their configurations are robust
against dynamical and secular instability.

  * •

Low-metallicity Gaia BHs (see Figure 2). There is no Gaia BBHs with the solar
metallicity. We make a caveat on this later.

  * •

Gaia BHs with more massive BHs than expected from theoretical models (see
Figure 4). This is simply because Gaia BBHs have two BHs. We discuss this
later.

As shown in Figure 2, we found no Gaia BBHs in the $Z=0.02$ model. Thus, we
claim that low-metallicity Gaia BHs are more likely to be Gaia BBHs than
solar-metallicity Gaia BHs. However, this may strongly depend on the binary
evolution model we adopt. The binary evolution model assumes that, when a star
in the Hertzsprung gap phase fills its Roche lobe and the mass transfer is
unstable, the star merges with its companion star, because such a star does
not have a steep density gradient between the helium core and hydrogen
envelope (Ivanova & Taam, 2004). This assumption largely reduces the formation
efficiency of BBHs in the solar-metallicity environments (Dominik et al.,
2012; Giacobbo et al., 2018; Tanikawa et al., 2022). If this assumption is not
appropriate, a large number of Gaia BBHs should be formed in the solar-
metallicity environments. Conversely, the discovery of a solar-metallicity
Gaia BBH would put strong constraints on BBH formation scenarios.

In the above guideline, we suggest that a Gaia BH can be a Gaia BBH if the BH
mass is unexpectedly large. However, this suggestion depends on the single-
star evolution and supernova models that we choose. Whenever unexpectedly
massive BHs have been discovered in the past, such as LB-1 (Liu et al.,
2019)333Several objections have been raised regarding the nature of the unseen
object of LB-1(Abdul-Masih et al., 2020; El-Badry & Quataert, 2020, 2021). and
GW190521 (Abbott et al., 2020a, b), modifications of single star evolution and
supernova models have been attempted (Takahashi, 2018; Belczynski et al.,
2020b; Farmer et al., 2020; Kinugawa et al., 2021; Vink et al., 2021; Tanikawa
et al., 2021; Costa et al., 2022; Kawashimo et al., 2023). Thus, it is
important to put observational constraints on the presence of Gaia BBHs (e.g.
Nagarajan et al., 2024). If a Gaia BBH is observationally discovered, it
should have significant impacts on single-star evolution and supernova models.

Finally, we estimate the total number of Gaia BBHs formed in open clusters in
the Milky Way galaxy. We assume that Gaia BBHs are located only in the thick
disk component of the Milky Way galaxy, because they have low metallicity.
Then, the number of such Gaia BBHs can be given by

$\displaystyle N_{\rm GaiaBBH}$ $\displaystyle\sim 60\left(\frac{\eta_{\rm
GaiaBBH}}{10^{-7}M_{\odot}}\right)\left(\frac{f_{\rm
cluster}}{0.1}\right)\left(\frac{M_{\rm thick}}{6\times
10^{9}M_{\odot}}\right),$

where $\eta_{\rm GaiaBBH}$ is the formation efficiency of Gaia BBHs in open
clusters (see Figure 2), $f_{\rm cluster}$ is the mass fraction of stars born
in open clusters (Misiriotis et al., 2006; Piskunov et al., 2007), and $M_{\rm
thick}$ are the mass of the thick disk component (Bland-Hawthorn & Gerhard,
2016). Currently, Gaia BH2 is the most distant among Gaia BHs, and its
distance is $\sim 1$ kpc (Tanikawa et al., 2023; El-Badry et al., 2023a).
Thus, we suppose that the Gaia’s detection limit for Gaia BHs is $1$ kpc. If
we assume that the thick disk size is $\sim 10$ kpc, the number of Gaia BBHs
is less than 1 within the Gaia’s detection limit. Gaia DR4 may extend the
detection limit, and the number of Gaia BBHs within the detection limit may
exceed 1. However, in order to discover Gaia BBHs, the completeness of the
detection of Gaia BBHs has to be nearly $100$ %. It may be hard to discover
Gaia BBHs if Gaia BHs are dominantly formed in open clusters.

## acknowledgments

This research could not be accomplished without the support by Grants-in-Aid
for Scientific Research (24K07040) from the Japan Society for the Promotion of
Science. M.F. is supported by The University of Tokyo Excellent Young
Researcher Program. L.W. thanks the support from the one-hundred-talent
project of Sun Yat-sen University, the Fundamental Research Funds for the
Central Universities, Sun Yat-sen University (22hytd09) and the National
Natural Science Foundation of China through grant 21BAA00619 and 12233013\.
Numerical simulations are carried out on Small Parallel Computers at Center
for Computational Astrophysics, National Astronomical Observatory of Japan,
the Yukawa Institute Computer Facility, and Cygnus/Pegasus at the CCS,
University of Tsukuba.

##  Appendix A Properties of Gaia BBHs

Figure 7 shows the formation efficiencies of BH binaries and triples with BBHs
as a function of orbital periods in the $Z=0.002$, $Z=0.005$ and $Z=0.01$
models. Figures 8, 9, and 10 draws the cumulative distributions of BH/BBH
masses, visible star masses, and eccentricities, respectively, for Gaia BHs
and Gaia BBHs in the $Z=0.002$, $Z=0.005$ and $Z=0.01$ models. We do not
include the $Z=0.02$ model in these Figures, since no Gaia BBHs are formed.
The results of these models are similar to those of the $Z=0.0002$ model as
seen in Figures 3, 4, 5, and 6, except for BH/BBH masses. BH/BBH masses
becomes smaller with metallicity increasing.

Figure 7: The same as Figure 3 except for the $Z=0.002$, $Z=0.005$ and
$Z=0.01$ models. Figure 8: The same as Figure 4 except for the $Z=0.002$,
$Z=0.005$ and $Z=0.01$ models. Figure 9: The same as Figure 5 except for the
$Z=0.002$, $Z=0.005$ and $Z=0.01$ models. Figure 10: The same as Figure 6
except for the $Z=0.002$, $Z=0.005$ and $Z=0.01$ models.

## References

  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) —. 2020b, ApJ, 900, L13
  * Abbott et al. (2023a) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023a, Physical Review X, 13, 041039
  * Abbott et al. (2023b) —. 2023b, Physical Review X, 13, 011048
  * Abdul-Masih et al. (2020) Abdul-Masih, M., Banyard, G., Bodensteiner, J., et al. 2020, Nature, 580, E11
  * Antonini & Rasio (2016) Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77
  * Arca Sedda et al. (2021) Arca Sedda, M., Li, G., & Kocsis, B. 2021, A&A, 650, A189
  * Balbinot et al. (2024) Balbinot, E., Dodd, E., Matsuno, T., et al. 2024, arXiv e-prints, arXiv:2404.11604
  * Banerjee et al. (2020) Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41
  * Belczynski et al. (2020a) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020a, A&A, 636, A104
  * Belczynski et al. (2020b) Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020b, ApJ, 890, 113
  * Belloni et al. (2017) Belloni, D., Giersz, M., Rocha-Pinto, H. J., Leigh, N. W. C., & Askar, A. 2017, MNRAS, 464, 4077
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29
  * Cappelluti et al. (2024) Cappelluti, N., Pacucci, F., & Hasinger, G. 2024, arXiv e-prints, arXiv:2406.07602
  * Casares et al. (2017) Casares, J., Jonker, P. G., & Israelian, G. 2017, X-Ray Binaries, ed. A. W. Alsabti & P. Murdin, 1499
  * Chakrabarti et al. (2023) Chakrabarti, S., Simon, J. D., Craig, P. A., et al. 2023, AJ, 166, 6
  * Chawla et al. (2022) Chawla, C., Chatterjee, S., Breivik, K., et al. 2022, ApJ, 931, 107
  * Claeys et al. (2014) Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83
  * Costa et al. (2022) Costa, G., Ballone, A., Mapelli, M., & Bressan, A. 2022, MNRAS, 516, 1072
  * Di Carlo et al. (2024) Di Carlo, U. N., Agrawal, P., Rodriguez, C. L., & Breivik, K. 2024, ApJ, 965, 22
  * Dominik et al. (2012) Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52
  * El-Badry (2024a) El-Badry, K. 2024a, arXiv e-prints, arXiv:2404.13047
  * El-Badry (2024b) —. 2024b, New A Rev., 98, 101694
  * El-Badry & Quataert (2020) El-Badry, K., & Quataert, E. 2020, MNRAS, 493, L22
  * El-Badry & Quataert (2021) —. 2021, MNRAS, 502, 3436
  * El-Badry et al. (2023a) El-Badry, K., Rix, H.-W., Cendes, Y., et al. 2023a, MNRAS, 521, 4323
  * El-Badry et al. (2023b) El-Badry, K., Rix, H.-W., Quataert, E., et al. 2023b, MNRAS, 518, 1057
  * El-Badry et al. (2024a) El-Badry, K., Simon, J. D., Reggiani, H., et al. 2024a, The Open Journal of Astrophysics, 7, 27
  * El-Badry et al. (2024b) El-Badry, K., Rix, H.-W., Latham, D. W., et al. 2024b, arXiv e-prints, arXiv:2405.00089
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91
  * Gaia Collaboration et al. (2023a) Gaia Collaboration, Arenou, F., Babusiaux, C., et al. 2023a, A&A, 674, A34
  * Gaia Collaboration et al. (2023b) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023b, A&A, 674, A1
  * Gaia Collaboration et al. (2024) Gaia Collaboration, Panuzzo, P., Mazeh, T., et al. 2024, arXiv e-prints, arXiv:2404.10486
  * Ganguly et al. (2023) Ganguly, A., Nayak, P. K., & Chatterjee, S. 2023, ApJ, 954, 4
  * Geier et al. (2023) Geier, S., Dorsch, M., Dawson, H., et al. 2023, A&A, 677, A11
  * Generozov & Perets (2024) Generozov, A., & Perets, H. B. 2024, ApJ, 964, 83
  * Giacobbo et al. (2018) Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959
  * Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15
  * Giesers et al. (2019) Giesers, B., Kamann, S., Dreizler, S., et al. 2019, A&A, 632, A3
  * Hayashi & Suto (2020) Hayashi, T., & Suto, Y. 2020, ApJ, 897, 29
  * Hayashi et al. (2023) Hayashi, T., Suto, Y., & Trani, A. A. 2023, ApJ, 958, 26
  * Hayashi et al. (2022) Hayashi, T., Trani, A. A., & Suto, Y. 2022, ApJ, 939, 81
  * Hayashi et al. (2020) Hayashi, T., Wang, S., & Suto, Y. 2020, ApJ, 890, 112
  * Hirai & Mandel (2022) Hirai, R., & Mandel, I. 2022, ApJ, 937, L42
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
  * Howil et al. (2024) Howil, K., Wyrzykowski, Ł., Kruszyńska, K., et al. 2024, arXiv e-prints, arXiv:2403.09006
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * Iorio et al. (2024) Iorio, G., Torniamenti, S., Mapelli, M., et al. 2024, arXiv e-prints, arXiv:2404.17568
  * Ivanova & Taam (2004) Ivanova, N., & Taam, R. E. 2004, ApJ, 601, 1058
  * Iwasawa et al. (2020) Iwasawa, M., Namekata, D., Nitadori, K., et al. 2020, PASJ, 72, 13
  * Iwasawa et al. (2016) Iwasawa, M., Tanikawa, A., Hosono, N., et al. 2016, PASJ, 68, 54
  * Kawashimo et al. (2023) Kawashimo, H., Sawada, R., Suwa, Y., et al. 2023, arXiv e-prints, arXiv:2306.01682
  * Kinugawa et al. (2021) Kinugawa, T., Nakamura, T., & Nakano, H. 2021, MNRAS, 501, L49
  * Kotko et al. (2024) Kotko, I., Banerjee, S., & Belczynski, K. 2024, arXiv e-prints, arXiv:2403.13579
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Kremer et al. (2018) Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 855, L15
  * Kroupa (1995a) Kroupa, P. 1995a, MNRAS, 277, 1491
  * Kroupa (1995b) —. 1995b, MNRAS, 277, 1507
  * Kroupa (2001) —. 2001, MNRAS, 322, 231
  * Küpper et al. (2011) Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300
  * Lam et al. (2022) Lam, C. Y., Lu, J. R., Udalski, A., et al. 2022, ApJ, 933, L23
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719
  * Liu et al. (2022) Liu, B., D’Orazio, D. J., Vigna-Gómez, A., & Samsing, J. 2022, Phys. Rev. D, 106, 123010
  * Liu et al. (2019) Liu, J., Zhang, H., Howard, A. W., et al. 2019, Nature, 575, 618
  * Liu et al. (2023) Liu, S., Wang, L., Hu, Y.-M., Tanikawa, A., & Trani, A. A. 2023, arXiv e-prints, arXiv:2311.05393
  * Lucy (2014) Lucy, L. B. 2014, A&A, 563, A126
  * Marín Pina et al. (2024) Marín Pina, D., Rastello, S., Gieles, M., et al. 2024, arXiv e-prints, arXiv:2404.13036
  * Martinez et al. (2022) Martinez, M. A. S., Rodriguez, C. L., & Fragione, G. 2022, ApJ, 937, 78
  * Martinez et al. (2020) Martinez, M. A. S., Fragione, G., Kremer, K., et al. 2020, ApJ, 903, 67
  * Misiriotis et al. (2006) Misiriotis, A., Xilouris, E. M., Papamastorakis, J., Boumis, P., & Goudis, C. D. 2006, A&A, 459, 113
  * Nagarajan et al. (2024) Nagarajan, P., El-Badry, K., Triaud, A. H. M. J., et al. 2024, PASP, 136, 014202
  * O’Neil et al. (2019) O’Neil, K. K., Martinez, G. D., Hees, A., et al. 2019, AJ, 158, 4
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Peters & Mathews (1963) Peters, P. C., & Mathews, J. 1963, Physical Review, 131, 435
  * Piskunov et al. (2007) Piskunov, A. E., Schilbach, E., Kharchenko, N. V., Röser, S., & Scholz, R. D. 2007, A&A, 468, 151
  * Rastello et al. (2023) Rastello, S., Iorio, G., Mapelli, M., et al. 2023, MNRAS, 526, 740
  * Rodriguez & Antonini (2018) Rodriguez, C. L., & Antonini, F. 2018, ApJ, 863, 7
  * Sahu et al. (2022) Sahu, K. C., Anderson, J., Casertano, S., et al. 2022, ApJ, 933, 83
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Santoliquido et al. (2020) Santoliquido, F., Mapelli, M., Bouffanais, Y., et al. 2020, ApJ, 898, 152
  * Seto (2013) Seto, N. 2013, Phys. Rev. Lett., 111, 061106
  * Shenar et al. (2022) Shenar, T., Sana, H., Mahy, L., et al. 2022, Nature Astronomy, 6, 1085
  * Shikauchi et al. (2020) Shikauchi, M., Kumamoto, J., Tanikawa, A., & Fujii, M. S. 2020, PASJ, 72, 45
  * Shikauchi et al. (2023) Shikauchi, M., Tsuna, D., Tanikawa, A., & Kawanaka, N. 2023, ApJ, 953, 52
  * Silsbee & Tremaine (2017) Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39
  * Takahashi (2018) Takahashi, K. 2018, ApJ, 863, 153
  * Tanikawa et al. (2024a) Tanikawa, A., Cary, S., Shikauchi, M., Wang, L., & Fujii, M. S. 2024a, MNRAS, 527, 4031
  * Tanikawa et al. (2023) Tanikawa, A., Hattori, K., Kawanaka, N., et al. 2023, ApJ, 946, 79
  * Tanikawa et al. (2021) Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2021, MNRAS, 505, 2170
  * Tanikawa et al. (2024b) Tanikawa, A., Wang, L., & Fujii, M. S. 2024b, The Open Journal of Astrophysics, 7, 39
  * Tanikawa et al. (2022) Tanikawa, A., Yoshida, T., Kinugawa, T., et al. 2022, ApJ, 926, 83
  * Toonen et al. (2018) Toonen, S., Perets, H. B., & Hamers, A. S. 2018, A&A, 610, A22
  * Totani et al. (2008) Totani, T., Morokuma, T., Oda, T., Doi, M., & Yasuda, N. 2008, PASJ, 60, 1327
  * Trani et al. (2022) Trani, A. A., Rastello, S., Di Carlo, U. N., et al. 2022, MNRAS, 511, 1362
  * Vigna-Gómez et al. (2021) Vigna-Gómez, A., Toonen, S., Ramirez-Ruiz, E., et al. 2021, ApJ, 907, L19
  * Vigna-Gómez et al. (2024) Vigna-Gómez, A., Willcox, R., Tamborra, I., et al. 2024, Phys. Rev. Lett., 132, 191403
  * Vink et al. (2021) Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146
  * Wang et al. (2020a) Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020a, MNRAS, 497, 536
  * Wang et al. (2020b) Wang, L., Nitadori, K., & Makino, J. 2020b, MNRAS, 493, 3398
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Yamaguchi et al. (2024) Yamaguchi, N., El-Badry, K., Fuller, J., et al. 2024, MNRAS, 527, 11719
  * Zheng et al. (2023) Zheng, L.-L., Sun, M., Gu, W.-M., et al. 2023, Science China Physics, Mechanics, and Astronomy, 66, 129512


# Anomalous circular phonon dichroism in transition metal dichalcogenides

Wen-Yu Shan  Department of Physics, School of Physics and Materials Science,
Guangzhou University, Guangzhou 510006, China

###### Abstract

Magnetic field can generally induce circular phonon dichroism based on the
formation of Landau levels of electrons. Here we study the magnetization-
induced circular phonon dichroism in transition metal dichalcogenides, without
forming the Landau levels. We find that, instead of the conventional
deformation potential coupling, the pseudogauge-type electron-phonon coupling
plays an essential role in the emergence of the phenomenon. As a concrete
example, a large dichroism signal is obtained in monolayer MoTe2 on a EuO
substrate, even without considering the Rashba spin-orbit coupling. Due to the
two-dimensional spin-valley-coupled band structure, MoTe2 shows a reciprocal
and nonreciprocal absorption of circularly polarized acoustic phonons upon
reversing the direction of phonon propagation and magnetization, respectively.
By varying the gate voltage, a tunable circular phonon dichroism can be
realized, which paves a way toward new physics and applications of two-
dimensional acoustoelectronics.

Introduction.—Recent years have seen a surge of interest in investigating
topological properties in the nonelectronic systems, e.g., photonic, magnonic
and phononic materials. For phonons, the concepts of band topology and
geometry have brought into new ingredients: chiral phonons Zhang and Niu
(2015); Zhu  _et al._ (2018); Chen  _et al._ (2019); Maity  _et al._ , angular
momentum Zhang and Niu (2014); Hamada  _et al._ (2018); Zhang and Murakami ;
Komiyama and Murakami (2021), orbital magnetic moments of phonons Juraschek
_et al._ (2017); Juraschek and Spaldin (2019); Cheng  _et al._ (2020);
Juraschek  _et al._ (2020), phonon angular momentum Hall effect Park and Yang
(2020), phonon rotoelectric effect Hamada and Murakami (2020) and so on. In
metals, the interplay between phonons and electrons with nontrivial band
topology or geometry may further induce distinctive features, such as phonon
helicity Hu  _et al._ (2021) and phonon magnetochiral effect Nomura  _et al._
(2019); Sengupta  _et al._ (2020).

Circular dichroism, the differential absorption between left- and right-handed
circularly polarized light, has been widely used in examining topological
phases of matter Souza and Vanderbilt (2008); Wang  _et al._ (2011); Tran  _et
al._ (2017); Liu  _et al._ (2018); Pozo  _et al._ (2019); Repellin and Goldman
(2019); Schüler  _et al._ (2020). A phononic analog, namely, circular phonon
dichroism (CPD), is later proposed in three-dimensional Weyl semimetals Liu
and Shi (2017). However, a direct analogy between phonons and photons is not
that obvious. The reasons are twofold. First, the photon wave vector is
usually much smaller than the Fermi wave vector of electrons, thus only
inducing the interband transition of electrons; whereas the phonon wave vector
may be comparable to that of electrons, giving rise to either interband or
intraband transition (see Fig. 1 (b)). Second, light waves consist of only
transverse modes, whereas acoustic waves in solids have both longitudinal and
transverse modes. Particularly, when dealing with two-dimensional (2D)
materials, one has to mix longitudinal and transverse in-plane modes to create
circular phonons Sonntag  _et al._ (2021), in marked contrast to the case of
light. This indicates that 2D circular phonon dichroism is intrinsically
different from the circular dichroism of light, where the former has received
far less attention.

Experimentally, several works have unveiled the effect of Landau levels of
electrons on the phonon dispersion or circular dichroism in graphene Kossacki
_et al._ (2012); Kumaravadivel  _et al._ (2019); Sonntag  _et al._ (2021),
such as the magnetophonon resonance. Nevertheless, the treatment of Landau
levels inevitably induces topology, even into an originally trivial system. In
this sense, the CPD can not resolve the real band topology or geometry of the
underlying system. Another way of breaking time-reversal symmetry is to
introduce the magnetic exchange interaction, which does not require the
formation of Landau levels and could retain the basic topology or geometry of
the band structure. Up to now, the intrinsic magnetization-induced CPD in 2D
materials like monolayer transition metal dichalcogenides, remains unknown.
This generalization of magnetization bears similarities to the case of
anomalous Hall effect, hence the name $anomalous$ $circular$ $phonon$
$dichroism$. The distinct spin-valley-coupled band structure of transition
metal dichalcogenides may further contribute to the anomalous behaviors of CPD
and their nonreciprocal relations. Therefore studying this new type of CPD
would be desirable for a better understanding and manipulation of band
geometry or topology in 2D materials.

In this paper, we explore the magnetization-induced CPD in monolayer
transition metal dichalcogenides. To allow this effect, the pseudogauge-type
electron-phonon coupling is necessary instead of the conventional deformation
potential coupling. We obtain a large dichroism signal in monolayer MoTe2 on a
EuO substrate, even in the absence of Rashba spin-orbit coupling. Due to the
unique spin-valley coupling, we find that MoTe2 shows a reciprocal
(nonreciprocal) absorption of circularly polarized acoustic phonons upon
reversing the direction of phonon propagation (magnetization). Our study
refreshes our knowledge on the effect of electron-phonon coupling on phonon
dynamics, and paves the way toward acoustoelectronics for 2D materials.

Figure 1:  Schematics of (a) the setup and (b) electronic band structure of
monolayer MoTe2. In (a), 2H-phase monolayer MoTe2 is deposited on the EuO
substrate. In (b), yellow (green) region corresponds to spin-up (-down) bands.
Transition process of electrons due to acoustic phonons (photons) is indicated
by the blue solid (red dashed) line.

Model Hamiltonian.—We take the pristine 2H-phase transition metal
dichalcogenides MoTe2 on a EuO substrate as a prototype (see Fig. 1 (a)). The
effective electronic Hamiltonian is given by
$\mathcal{H}_{e}=\sum_{\bm{k}}\psi^{+}(\bm{k})[H_{0}+H_{soc}+H_{ex}+H_{R}]\psi(\bm{k})$,
where Qi  _et al._ (2015); Habe and Koshino (2017)

$\begin{split}H_{0}&=\hbar
v(\tau\sigma_{x}k_{x}+\sigma_{y}k_{y})+\frac{\Delta}{2}\sigma_{z},\\\
H_{soc}&=\tau s_{z}(\lambda_{c}\sigma_{+}+\lambda_{v}\sigma_{-}),\\\
H_{ex}&=-\bm{s}\cdot\bm{n}(B_{c}\sigma_{+}+B_{v}\sigma_{-}),\\\
H_{R}&=\lambda_{R}(\tau s_{y}\sigma_{x}-s_{x}\sigma_{y}).\end{split}$ (1)

$\psi^{+}(\bm{k})$ and $\psi(\bm{k})$ are the creation and annihilation
operator of electrons. $H_{soc}$, $H_{ex}$ and $H_{R}$ correspond to the
Ising-type spin-orbit coupling, proximity-induced exchange and Rashba
interaction, respectively. $\bm{s}$ and $\bm{\sigma}$ are Pauli matrices
acting on spin $\\{\uparrow,\downarrow\\}$ and orbit subspace
$\\{|d_{z^{2}}\rangle,\frac{1}{\sqrt{2}}(|d_{x^{2}-y^{2}}\rangle+i\tau|d_{xy}\rangle)\\}$,
and $\sigma_{\pm}=\frac{1}{2}(\sigma_{0}\pm\sigma_{z})$. $\tau=\pm 1$ labels
valley $K_{\pm}$. $\lambda_{c/v}$ describes the spin splitting of the
conduction and valence bands, respectively. $B_{c/v}$ is the effective Zeeman
field experienced by the conduction and valence bands, arising from the
exchange coupling with the magnetic substrate. The out-of-plane $z$-direction
magnetization $\bm{n}=\bm{e}_{z}$ is considered (see Fig. 1 (a)). For the
moment, we set $\lambda_{R}=0$ in order to have analytical expressions and
intuitive physical picture. The role of $\lambda_{R}$ will be clarified later.
The electronic band structure upon magnetization is schematically shown in
Fig. 1 (b), where the signature of spin-valley coupling can be seen
explicitly. The Fermi level is pinned at the valence bands, where the effect
of spin-valley coupling is manifest. The magnetization $B_{c/v}$ shifts the
opposite-spin states from different valleys in opposite directions, and thus
breaks the time-reversal symmetry.

For the phononic part, we consider two branches of in-plane acoustic phonon
modes. Due to the low sound velocity $c_{l/t}$ ($l/t$ for
longitudinal/transverse phonon polarization), the acoustic phonon energy
$\omega_{l/t}=\hbar c_{l/t}|\bm{q}|$ is much smaller than the valence band
splitting of electrons, i.e., $2(B_{v}\pm\lambda_{v})$, where $\bm{q}$ is the
phonon wave vector. As a result, only intraband transitions of electrons are
triggered by acoustic phonons (see Fig. 1 (b)). By contrast, optical phonon
modes with larger energy, may enable either intraband or interband
transitions. Nonetheless, the basic physical picture should be similar. For
simplicity, we further study the long-wavelength limit of phonon modes, which
allows us to neglect the intervalley scattering process of electrons.

Based on the theory of elasticity lan ; Suzuura and Ando (2002); Cazalilla
_et al._ (2014), the electron-acoustic-phonon coupling in MoTe2 contains two
terms: $\mathcal{H}_{e-ph}=\mathcal{H}_{e-ph}^{d}+\mathcal{H}_{e-ph}^{p}$,
where $\mathcal{H}_{e-ph}^{d}$ ($\mathcal{H}_{e-ph}^{p}$) refers to the
deformation (pseudogauge) potential coupling. $\mathcal{H}_{e-ph}^{d/p}$ has a
general form Shan (2020); sup

$\begin{split}\mathcal{H}_{e-ph}^{d/p}&=\sum_{\bm{k},\bm{q}}\psi^{+}(\bm{k}+\bm{q})[\bm{u}(\bm{q})\cdot\hat{\bm{T}}_{d/p}(\bm{q})]\psi(\bm{k}),\end{split}$
(2)

where $\bm{u}(\bm{q})$ is the Fourier transform of the in-plane collective
displacement $\bm{u}(\bm{r})$ for acoustic modes Suzuura and Ando (2002) and
$\hat{\bm{T}}(\bm{q})$ is the Fourier transform of the “effective” force
operator $\hat{\bm{T}}(\bm{r})$ acting on atoms by electrons. For the
deformation potential $\mathcal{H}_{e-ph}^{d}$,
$\hat{\bm{T}}_{d}(\bm{q})=ig_{d}\bm{q}$, which is independent of the valley
index $\tau$. For the pseudogauge potential $\mathcal{H}_{e-ph}^{p}$, the
force operator becomes valley-dependent, that is,
$\hat{\bm{T}}^{\tau=-1}_{p}(\bm{q})=ig_{p}[\bm{q}\cdot\bm{\sigma},(\bm{q}\times\bm{\sigma})_{z}]$
and
$\hat{\bm{T}}^{\tau=1}_{p}(\bm{q})=\mathcal{K}[\hat{\bm{T}}^{\tau=-1}_{p}(-\bm{q})]$.
$\mathcal{K}$ is the complex conjugation operator. The relation between
$\hat{\bm{T}}_{p}^{\tau=1}(\bm{q})$ and $\hat{\bm{T}}_{p}^{\tau=-1}(-\bm{q})$
preserves the time-reversal symmetry of electron-phonon coupling in the
absence of magnetization.

Phonon equation of motion.—For the phonon dynamics, we consider the phonon
equation of motion in the frequency-momentum ($\omega$, $\bm{q}$) domain Liu
and Shi (2017)

$\omega^{2}u_{\alpha}(\bm{q})=\sum_{\beta}[\Phi_{\alpha\beta}(\bm{q})+\hbar\chi_{\alpha\beta}(\bm{q},\omega)]u_{\beta}(\bm{q}),\;$
(3)

where $\alpha,\beta=x,y$ and $\Phi(\bm{q})$ is the dynamical matrix.
$\chi_{\alpha\beta}(\bm{q},\omega)$ is a retarded response function arising
from the electron-phonon coupling and follows at each valley sup

$\begin{split}&\chi_{\alpha\beta}^{\tau}(\bm{q},\omega+i\delta)=\sum_{n,m}\int\frac{\hbar
d^{2}\bm{k}}{\rho(2\pi)^{2}}\frac{f_{\tau,m,\bm{k}}-f_{\tau,n,\bm{k}-\bm{q}}}{\omega+i\delta+E_{\tau,m,\bm{k}}-E_{\tau,n,\bm{k}-\bm{q}}}\\\
&\times\langle\tau,m,\bm{k}|\hat{T}_{\alpha}^{\tau}(\bm{q})|\tau,n,\bm{k}-\bm{q}\rangle\langle\tau,n,\bm{k}-\bm{q}|\hat{T}_{\beta}^{\tau}(-\bm{q})|\tau,m,\bm{k}\rangle.\end{split}$
(4)

$E_{\tau,m,\bm{k}}$ and $|\tau,m,\bm{k}\rangle$ are the dispersion and
electronic wave function of Hamiltonian (1), respectively. $f_{\tau,m,\bm{k}}$
($f_{\tau,n,\bm{k}-\bm{q}}$) is the Fermi distribution function, $\rho$ is the
2D mass density, $\delta$ is a positive infinitesimal. Since only intraband
transitions (band indices $m=n$) of electrons are allowed by acoustic modes in
the low-temperature limit, $m,n$ reduce to the ones intersected by the Fermi
level, i.e., spin-split valence bands at valley $K_{\pm}$ (see Fig. 1 (b)).

Circular phonon dichroism.—Our main interest lies in the anti-Hermitian part
of $\chi(\bm{q},\omega)$, that is, $-2i\omega\gamma(\bm{q},\omega)$, where
$\gamma(\bm{q},\omega)$ is a Hermitian matrix satisfying
$\gamma^{+}(\bm{q},\omega)=\gamma(\bm{q},\omega)$. This matrix corresponds to
the non-Hermitian part of the phonon self-energy, which physically originates
from the phonon absorption by electrons. In the basis of
$\\{\hat{x},\hat{y}\\}^{T}$, $\gamma$ matrix has the form

$\begin{split}\gamma(\bm{q},\omega)=\left[\begin{array}[]{cc}D(\bm{q},\omega)+\bar{D}(\bm{q},\omega)&\bar{A}(\bm{q},\omega)+iA(\bm{q},\omega)\\\
\\\
\bar{A}(\bm{q},\omega)-iA(\bm{q},\omega)&D(\bm{q},\omega)-\bar{D}(\bm{q},\omega)\\\
\end{array}\right].\end{split}$ (5)

Different from the Weyl semimetals Liu and Shi (2017), new terms
$\bar{D}(\bm{q},\omega)$ and $\bar{A}(\bm{q},\omega)$ occur in monolayer MoTe2
as a result of $D_{3h}$ point-group symmetry. For the left- and right-handed
circularly polarized phonons, $|u_{L/R}\rangle=$ $\frac{1}{\sqrt{2}}[1$ $\pm
i]^{T}$, the damping (absorption) coefficients read $\gamma^{L/R}=$
$D(\bm{q},\omega)\mp A(\bm{q},\omega)$. The relative difference between
$\gamma^{L}$ and $\gamma^{R}$ defines the $circular$ $phonon$ $dichroism$
(CPD). One can see that the behavior of CPD is totally determined by
$A(\bm{q},\omega)/D(\bm{q},\omega)$.

For longitudinal or transverse phonons, the polarization is linear as
$|u_{l}\rangle=$ $[\cos\phi_{\bm{q}}$ $\sin\phi_{\bm{q}}]^{T}$ and
$|u_{t}\rangle=$ $[-\sin\phi_{\bm{q}}$ $\cos\phi_{\bm{q}}]^{T}$, where the
angular variable $\phi_{\bm{q}}=\tan^{-1}(q_{y}/q_{x})$. The damping
coefficients are given by $\gamma^{l/t}=$ $D(\bm{q},\omega)$ $\pm$ $\cos
2\phi_{\bm{q}}\bar{D}(\bm{q},\omega)$ $\pm$ $\sin
2\phi_{\bm{q}}\bar{A}(\bm{q},\omega)$, which explicitly depends on the phonon
propagation direction $\bm{q}$. Here, different from the circular phonons, the
damping coefficients $\gamma^{l/t}$ for the linear phonons depend on the
parameters $\bar{D}(\bm{q},\omega)$ and $\bar{A}(\bm{q},\omega)$.

Specifically for the deformation potential $\mathcal{H}_{e-ph}^{d}$, $\gamma$
matrix is proportional to sup

$\begin{split}\gamma(\bm{q},\omega)\propto\left[\begin{array}[]{cc}q_{x}^{2}&q_{x}q_{y}\\\
\\\ q_{x}q_{y}&q_{y}^{2}\\\ \end{array}\right].\end{split}$ (6)

This immediately leads to $A(\bm{q},\omega)=0$, meaning that the CPD vanishes
when only the deformation potential coupling is taken in account. Meanwhile,
$\gamma^{t}=0$, suggesting that there is no absorption for the transverse
phonon modes. This agrees with the fact that the deformation potential only
couples electrons to the longitudinal phonon modes Shan (2020).

Figure 2:  (a) Angular $\phi_{\bm{q}}$-dependence of the damping coefficients
$\gamma^{l}$ ($\gamma^{t}$) for the longitudinal (transverse) acoustic phonon
modes. (b)-(d) Relations of the circular phonon dichroism $A/D$ versus the
phonon wave vector $q/k_{F}^{\tau=-1}$, Fermi energy $E_{F}$ and gap function
$\Delta$, respectively. In (b), the Fermi energy is fixed: $E_{F}=-0.48$ eV.
The inset shows the details of the cyan region, and the two peaks are due to
the longitudinal and transverse phonon modes, respectively. The peaks in the
cyan region are given by a summation of valleys $K_{\pm}$, whereas the peaks
near $q/k_{F}^{\tau=-1}=2$ are only determined by valley $K_{-}$.
$k_{F}^{\tau=-1}$ is the Fermi wave vector at valley $K_{-}$. The red dot
refers to the case of (a). In (c), different values of $q$ are adopted. The
locations of the Fermi level at the peaks are shown in the inset: both valley
$K_{\pm}$ are intersected at smaller $E_{F}$; only valley $K_{-}$ is
intersected at larger $E_{F}$. In (d), the black dot indicates the value of
$\Delta$ in (a)-(c): $\Delta=1.05$ eV. Parameters: $\lambda_{v}=0.11$ eV,
$\lambda_{c}=0.029$ eV, $\hbar v=2.33$ eV$\cdot\AA$, $B_{c}=0.206$ eV,
$B_{v}=0.17$ eV Qi  _et al._ (2015), longitudinal and transverse sound
velocity $c_{l}=3.64\times 10^{3}$ m/s and $c_{t}=2.21\times 10^{3}$ m/s Rano
_et al._ (2020), mass density $\rho=9.40\times 10^{-6}$ kg/m2 Keum  _et al._
(2015) and the electron-phonon coupling constant $g_{p}=0.32$ eV Chen  _et
al._ (2016).

For the pseudogauge potential $\mathcal{H}_{e-ph}^{p}$, the situation is more
complex. When only focusing on the acoustic modes, analytical expressions for
all elements of $\gamma(\bm{q},\omega)$ matrix can be obtained sup . For
example, when a single valence band at valley $K_{\tau}$ is intersected by the
Fermi level, $A(\bm{q},\omega)/D(\bm{q},\omega)$ reduces to

$\begin{split}\frac{A^{\tau}(\bm{q},\omega)}{D^{\tau}(\bm{q},\omega)}&=\tau\omega\frac{\Delta-\tau\lambda_{c}+\tau\lambda_{v}+B_{c}-B_{v}}{2[\omega
x_{F}^{\tau}-(\hbar
vk_{F}^{\tau})^{2}]}\Theta(k_{F}^{\tau}-\frac{q}{2}-k_{0}^{\tau}),\end{split}$
(7)

where the Heaviside step function $\Theta(\cdots)$ constrains the magnitude of
phonon wave vector $q=|\bm{q}|$. $k_{F}^{\tau}$ is a valley-dependent Fermi
wave vector of electrons. $k_{0}^{\tau}=$ $\frac{\omega}{2\hbar
v}\sqrt{1+\frac{(\Delta-\tau\lambda_{c}+\tau\lambda_{v}+B_{c}-B_{v})^{2}}{(\hbar
vq)^{2}-\omega^{2}}}$ and $x_{F}^{\tau}=$
$\sqrt{(\frac{\Delta-\tau\lambda_{c}+\tau\lambda_{v}+B_{c}-B_{v}}{2})^{2}+(\hbar
vk_{F}^{\tau})^{2}}$. However, such a simple relation fails when both valleys
are intersected by the Fermi level, given that
$A(\bm{q},\omega)=\sum_{\tau}A^{\tau}(\bm{q},\omega)$ and
$D(\bm{q},\omega)=\sum_{\tau}D^{\tau}(\bm{q},\omega)$. On the other hand, both
$\bar{A}(\bm{q},\omega)$ and $\bar{D}(\bm{q},\omega)$ become
$\phi_{\bm{q}}$-dependent sup : $\bar{A}(\bm{q},\omega)=-F(q,\omega)\sin
4\phi_{\bm{q}}$ and $\bar{D}(\bm{q},\omega)=F(q,\omega)\cos 4\phi_{\bm{q}}$,
with a $\phi_{\bm{q}}$-independent factor $F(q,\omega)$. By substituting these
into $\gamma^{l/t}$, we find for linearly polarized phonons,

$\begin{split}\gamma^{l/t}=D(q,\omega)\pm F(q,\omega)\cos
6\phi_{\bm{q}}.\end{split}$ (8)

One can see that $\gamma^{l/t}$ has a six-fold ($C_{6}$) rotational symmetry
on $\phi_{\bm{q}}$ (see Fig. 2 (a)), which is different from the three-fold
($C_{3}$) rotational symmetry of the underlying crystals. The reason for the
symmetry mismatch is due to the reciprocal behaviors of $\gamma^{l/t}$ upon
reversing the direction of phonon propagation, i.e.,
$\bm{q}\rightarrow-\bm{q}$, as shown in Table 1. For phonons, $c_{l/t}\ll v$,
giving rise to sup  $D(q,\omega)\approx-F(q,\omega)$. As a result,
$\gamma^{l}\approx 2D(q,\omega_{l})\sin^{2}3\phi_{\bm{q}}$ and
$\gamma^{t}\approx 2D(q,\omega_{t})\sin^{2}3(\phi_{\bm{q}}-\frac{\pi}{6})$.
This means that there is an angular shift $\frac{\pi}{6}$ in $\phi_{\bm{q}}$
between $\gamma^{l}$ and $\gamma^{t}$, as shown in Fig. 2 (a).

Table 1:  Transformation properties of parameters $D$, $\bar{D}$, $A$ and $\bar{A}$. Transformation | $D(\bm{q},\omega)$ | $\bar{D}(\bm{q},\omega)$ | $A(\bm{q},\omega)$ | $\bar{A}(\bm{q},\omega)$
---|---|---|---|---
$\bm{q}\rightarrow-\bm{q}$ | $+$ | $+$ | $+$ | $+$
$\bm{n}\rightarrow-\bm{n}$ | $+$ | $+$ | $-$ | $+$

For circularly polarized phonons, numerical results of
$A(\bm{q},\omega)/D(\bm{q},\omega)$ as functions of the rescaled phonon wave
vector $q/k_{F}^{\tau=-1}$, Fermi energy $E_{F}$ and $\Delta$ are shown in
Fig. 2 (b)-(d), respectively. The Fermi wave vector $k_{F}^{\tau=-1}$ rather
than $k_{F}^{\tau=1}$ is selected since the valence band edge of valley
$K_{-}$ is higher than $K_{+}$, as shown in Fig. 1. In Fig. 2 (b), a non-
monotonic behavior of $A/D$ as $q$ increases can be seen explicitly. The jumps
at $q/k_{F}^{\tau=-1}\approx 0.47$ and $1.98$ originate from the sudden
vanishing of valley $K_{+}$ and $K_{-}$, respectively, as required by the
factor $\Theta(k_{F}^{\tau}-\frac{q}{2}-k_{0}^{\tau})$ in Eq. (7). Such a
factor can be understood as a result of the energy and momentum conservation
for the electron-phonon scattering process. For acoustic phonons, the electron
scattering approximately occurs on the Fermi surface. In this sense, the
phonon wave vector $q$ must be smaller than the maximum value of momentum
transfer of electrons, that is, $q<2k_{F}^{\tau}$. $k_{0}^{\tau}$ is a small
offset wave vector arising from the acoustic phonon dispersion $\omega$. As
seen in the inset of Fig. 2 (b), there are actually two adjacent peaks (jumps)
in the highlighted region corresponding to the $l$ and $t$ mode, respectively,
since $k_{0}^{\tau}$ is different for $\omega=\omega_{l/t}$. As the sound
velocity $c_{l}>c_{t}$, $k_{0}^{\tau}$ is larger for the longitudinal mode,
leading to a smaller transition value of $q$. In Fig. 2 (c), the locations of
the Fermi level for the peaks are indicated in the inset. The peaks at the
lower (higher) Fermi level are dominated by valley $K_{+}$ ($K_{-}$), which
exhibit opposite signs of $A/D$. For each valley, the magnitude $|A/D|$
increases when the Fermi level is tuned toward the band edge. Different values
of $q$ are also compared. We find that by adopting a smaller $q$, the peaks
are shifted to a higher Fermi level, as $k_{F}^{\tau}$ becomes smaller. The
peaks also show a larger magnitude and become sharper, particularly for the
second peaks. Therefore this provides a means of tuning the sign and magnitude
of the CPD. In Fig. 2 (d), the value of $\Delta$ adopted in Fig. 2 (a)-(c) is
indicated. We can see that the magnitude $|A/D|$ is basically enhanced when
$\Delta$ increases, expect for the discontinuous points. That is the reason
why we propose monolayer transition metal dichalcogenides as candidate
materials, which have large band gap and thus large CPD signals. For an order-
of-magnitude estimate, we consider the parameters corresponding to the red dot
in Fig. 2 (b), which also refer to the case of Fig. 2 (a). We find
$D=1.90\times 10^{7}$/s and $A=-7.81\times 10^{5}$/s. This yields a difference
of the attenuation between the left- and right-handed circularly polarized
waves, that is, $(\gamma^{L}-\gamma^{R})/\bar{c}\sim 534$/m, where
$\bar{c}=(c_{l}+c_{t})/2$ is the average sound velocity. Such difference is
much larger than that of the Weyl semimetals Liu and Shi (2017), and should be
observable in ultrasonic experiments.

Nonreciprocal absorption.—Given that both the space-inversion and time-
reversal symmetry are broken in our system, the absorption of circularly
polarized phonons is expected to be nonreciprocal. To see this, we consider in
Table 1 the transformation properties of parameters $D$, $\bar{D}$, $A$ and
$\bar{A}$ upon reversing the direction of phonon propagation $\bm{q}$ or
magnetization $\bm{n}$. We find that $D$, $\bar{D}$ and $\bar{A}$ are even
functions of $\bm{q}$ and $\bm{n}$, whereas $A$ is an even (odd) function of
$\bm{q}$ ($\bm{n}$). Accordingly, the absorption coefficients of circular
phonons $\gamma^{L/R}$ remain unchanged under the transformation
$\bm{q}\rightarrow-\bm{q}$, whereas $\gamma^{L/R}$ interchange with each other
under the transformation $\bm{n}\rightarrow-\bm{n}$. This represents a
reciprocal and nonreciprocal CPD upon reversing the direction of phonon
propagation and magnetization, respectively. Such result is similar to that of
the Faraday rotation of light polarization Freiser (1968), where the rotation
angle only depends on the magnetic field direction. However, the origin is
different. The absorption of circularly polarized photons is actually
nonreciprocal when $\bm{q}\rightarrow-\bm{q}$, but the chirality of circular
photons also depends on the light propagation direction $\bm{q}$. As a result,
the reciprocal $\bm{q}$-dependence of the rotation angle is recovered. On the
other hand, due to the 2D nature, the chirality of in-plane circular phonons
is independent of the phonon propagation direction $\bm{q}$, giving rise to
the reciprocal absorption. This also indicates that there is no directional
dichroism Fuchs (1965) or phonon magnetochiral effect Nomura  _et al._ (2019);
Sengupta  _et al._ (2020) in our system.

Figure 3:  (a) Electronic band structure of monolayer MoTe2 and (b) circular
phonon dichroism $A/D$ versus the Fermi energy $E_{F}$ for different strength
of Rashba spin-orbit coupling.

Roles of Rashba spin-orbit coupling.—Now we take into account the Rashba term
$H_{R}$ and treat it as a perturbation. Analytical expressions are calculated
sup  and numerical results of electronic band structure and CPD are shown in
Fig. 3. In Fig. 3 (a), we find that the Rashba term shifts the conduction
(valence) band edge to higher (lower) energy at valley $K_{-}$, whereas it
hardly changes the band structure at valley $K_{+}$. This explains the
phenomenon of peak shift in Fig. 3 (b) since the peaks are always close to the
band edge. By increasing the strength of Rashba spin-orbit coupling, the
magnitude of CPD can be enhanced, which provides us a knob to tune the CPD.
For a realistic strength $\lambda_{R}=0.072$ eV Qi  _et al._ (2015), the
behavior of $A/D$ is similar to the case without Rashba spin-orbit coupling,
thus validating our above treatment. Particularly, we find that the
introduction of $H_{R}$ does not change the reciprocal behaviors of absorption
coefficients $\gamma^{L/R}$ under $\bm{q}\rightarrow-\bm{q}$. Therefore, to
obtain the nonreciprocity, additional ingredients should be taken into
account, such as the cyclotron motion of electrons Heil  _et al._ (1982);
Kossacki  _et al._ (2012); Kumaravadivel  _et al._ (2019); Liu  _et al._
(2020); Sonntag  _et al._ (2021) or phonon-magnon coupling Shan  _et al._
(2020).

Discussion and conclusion.—We have studied the circular phonon dichroism in
magnetic two-dimensional materials, i.e., monolayer MoTe2 in proximity to the
EuO substrate. Large dichroism signal is obtained for the pseudogauge-type
electron-phonon coupling, even without introducing the Landau levels or Rashba
spin-orbit coupling. Such a signal is reciprocal (nonreciprocal) upon
reversing the direction of phonon propagation (magnetization). By varying the
gate voltage, the CPD signal can be tuned through the role of Fermi level and
Rashba spin-orbit coupling.

The proposed CPD effect can also be applied to other transition-metal
dichalcogenides with spin-valley-coupled band structure, or their van der
waals heterostructures. The effect can be detected by the pulse-echo technique
tru ; lut  based on the different absorption coefficients between left- and
right-handed circularly polarized phonons. An alternative detection is the
Raman spectroscopy analysis Kossacki  _et al._ (2012); Sonntag  _et al._
(2021) of phonon polarization by injecting a linearly polarized acoustic
waves.

Acknowledgments. This work is supported by the National Natural Science
Foundation of China (NSFC, Grant No. 11904062). We also acknowledge the
support of a startup grant from Guangzhou University.

## References

  * Zhang and Niu (2015) L. Zhang and Q. Niu, Phys. Rev. Lett. 115, 115502 (2015).
  * Zhu  _et al._ (2018) H. Zhu, J. Yi, M.-Y. Li, J. Xiao, L. Zhang, C.-W. Yang, R. A. Kaindl, L.-J. Li, Y. Wang, and X. Zhang, Science 359, 579 (2018).
  * Chen  _et al._ (2019) X. Chen, X. Lu, S. Dubey, Q. Yao, S. Liu, X. Wang, Q. Xiong, L. Zhang, and A. Srivastava, Nature Physics 15, 221 (2019).
  * (4) I. Maity, A. A. Mostofi, and J. Lischner, arXiv:2108.03965 .
  * Zhang and Niu (2014) L. Zhang and Q. Niu, Phys. Rev. Lett. 112, 085503 (2014).
  * Hamada  _et al._ (2018) M. Hamada, E. Minamitani, M. Hirayama, and S. Murakami, Phys. Rev. Lett. 121, 175301 (2018).
  * (7) T. Zhang and S. Murakami, arXiv:2107.04243 .
  * Komiyama and Murakami (2021) H. Komiyama and S. Murakami, Phys. Rev. B 103, 214302 (2021).
  * Juraschek  _et al._ (2017) D. M. Juraschek, M. Fechner, A. V. Balatsky, and N. A. Spaldin, Phys. Rev. Materials 1, 014401 (2017).
  * Juraschek and Spaldin (2019) D. M. Juraschek and N. A. Spaldin, Phys. Rev. Materials 3, 064405 (2019).
  * Cheng  _et al._ (2020) B. Cheng, T. Schumann, Y. C. Wang, X. S. Zhang, D. Barbalas, S. Stemmer, and N. P. Armitage, Nano Lett. 20, 5991 (2020).
  * Juraschek  _et al._ (2020) D. M. Juraschek, P. Narang, and N. A. Spaldin, Phys. Rev. Research 2, 043035 (2020).
  * Park and Yang (2020) S. Park and B.-J. Yang, Nano Lett. 20, 7694 (2020).
  * Hamada and Murakami (2020) M. Hamada and S. Murakami, Phys. Rev. B 101, 144306 (2020).
  * Hu  _et al._ (2021) L.-H. Hu, J. Yu, I. Garate, and C.-X. Liu, Phys. Rev. Lett. 127, 125901 (2021).
  * Nomura  _et al._ (2019) T. Nomura, X.-X. Zhang, S. Zherlitsyn, J. Wosnitza, Y. Tokura, N. Nagaosa, and S. Seki, Phys. Rev. Lett. 122, 145901 (2019).
  * Sengupta  _et al._ (2020) S. Sengupta, M. N. Y. Lhachemi, and I. Garate, Phys. Rev. Lett. 125, 146402 (2020).
  * Souza and Vanderbilt (2008) I. Souza and D. Vanderbilt, Phys. Rev. B 77, 054438 (2008).
  * Wang  _et al._ (2011) Y. H. Wang, D. Hsieh, D. Pilon, L. Fu, D. R. Gardner, Y. S. Lee, and N. Gedik, Phys. Rev. Lett. 107, 207602 (2011).
  * Tran  _et al._ (2017) D. T. Tran, A. Dauphin, A. G. Grushin, P. Zoller, and N. Goldman, Sci. Adv. 3, e1701207 (2017).
  * Liu  _et al._ (2018) Y. Liu, S. A. Yang, and F. Zhang, Phys. Rev. B 97, 035153 (2018).
  * Pozo  _et al._ (2019) O. Pozo, C. Repellin, and A. G. Grushin, Phys. Rev. Lett. 123, 247401 (2019).
  * Repellin and Goldman (2019) C. Repellin and N. Goldman, Phys. Rev. Lett. 122, 166801 (2019).
  * Schüler  _et al._ (2020) M. Schüler, U. D. Giovannini, H. Hübener, A. Rubio, M. A. Sentef, and P. Werner, Sci. Adv. 6, eaay2730 (2020).
  * Liu and Shi (2017) D. Liu and J. Shi, Phys. Rev. Lett. 119, 075301 (2017).
  * Sonntag  _et al._ (2021) J. Sonntag, S. Reichardt, B. Beschoten, and C. Stampfer, Nano Lett. 21, 2898 (2021).
  * Kossacki  _et al._ (2012) P. Kossacki, C. Faugeras, M. Kühne, M. Orlita, A. Mahmood, E. Dujardin, R. R. Nair, A. K. Geim, and M. Potemski, Phys. Rev. B 86, 205431 (2012).
  * Kumaravadivel  _et al._ (2019) P. Kumaravadivel, M. T. Greenaway, D. Perello, A. Berdyugin, B. J., J. Wengraf, S. Liu, J. H. Edgar, A. K. Geim, L. Eaves, and R. K. Kumar, Nat. Commun. 10, 3334 (2019).
  * Qi  _et al._ (2015) J. Qi, X. Li, Q. Niu, and J. Feng, Phys. Rev. B 92, 121403 (2015).
  * Habe and Koshino (2017) T. Habe and M. Koshino, Phys. Rev. B 96, 085411 (2017).
  * (31) L. D. Landau and E. M. Lifschitz, _Theory of Elasticity_ , Pergamon Press, Oxford, 1959.
  * Suzuura and Ando (2002) H. Suzuura and T. Ando, Phys. Rev. B 65, 235412 (2002).
  * Cazalilla  _et al._ (2014) M. A. Cazalilla, H. Ochoa, and F. Guinea, Phys. Rev. Lett. 113, 077201 (2014).
  * Shan (2020) W.-Y. Shan, Phys. Rev. B 102, 241301 (2020).
  * (35) See Supplementary Material for calculation details.
  * Rano  _et al._ (2020) B. R. Rano, I. M. Syed, and S. H. Naqib, Results in Physics 19, 103639 (2020).
  * Keum  _et al._ (2015) D. H. Keum, S. Cho, J. H. Kim, D.-H. Choe, H.-J. Sung, M. Kan, H. Kang, J.-Y. Hwang, S. W. Kim, H. Yang, K. J. Chang, and Y. H. Lee, Nature Physics 11, 482 (2015).
  * Chen  _et al._ (2016) F. C. Chen, X. Luo, R. C. Xiao, W. J. Lu, B. Zhang, H. X. Yang, J. Q. Li, Q. L. Pei, D. F. Shao, R. R. Zhang, L. S. Ling, C. Y. Xi, W. H. Song, and Y. P. Sun, Appl. Phys. Lett. 108, 162601 (2016).
  * Freiser (1968) M. Freiser, IEEE Trans. Magn. 4, 152 (1968).
  * Fuchs (1965) R. Fuchs, Phil. Mag. 11, 647 (1965).
  * Heil  _et al._ (1982) J. Heil, B. Lüthi, and P. Thalmeier, Phys. Rev. B 25, 6515 (1982).
  * Liu  _et al._ (2020) Z. Liu, K. Guo, G. Hu, Z. Shi, Y. Li, L. Zhang, H. Chen, L. Zhang, P. Zhou, H. Lu, M.-L. Lin, S. Liu, Y. Cheng, X. L. Liu, J. Xie, L. Bi, P.-H. Tan, L. Deng, C.-W. Qiu, and B. Peng, Sci. Adv. 6, eabc7628 (2020).
  * Shan  _et al._ (2020) P. J. Shan, D. A. Bas, I. Lisenkov, A. Matyushov, N. X. Sun, and M. R. Page, Sci. Adv. 6, eabc5648 (2020).
  * (44) R. Truell, C. Elbaum, and B. B. Chick, _Ultrasonic Methods in Solid State Physics_ (Academic Press, New York, 1969).
  * (45) B. Lüthi, _Physical Acoustics in the Solid State_ (Springer-Verlag, Berlin, 2004).