# Thermodynamic determination of the equilibrium first-order phase-transition
line hidden by hysteresis in a phase diagram

Keisuke Matsuura<EMAIL_ADDRESS>RIKEN Center for Emergent Matter
Science, Wako 351-0198, Japan  Yo Nishizawa  Department of Applied Physics and
Quantum-Phase Electronics Center (QPEC), University of Tokyo, Tokyo 113-8656,
Japan  Markus Kriener  RIKEN Center for Emergent Matter Science, Wako
351-0198, Japan  Takashi Kurumaji  Department of Advanced Materials Science,
University of Tokyo, Kashiwa 277-8561, Japan  Hiroshi Oike  RIKEN Center for
Emergent Matter Science, Wako 351-0198, Japan  Department of Applied Physics
and Quantum-Phase Electronics Center (QPEC), University of Tokyo, Tokyo
113-8656, Japan  PRESTO, Japan Science and Technology Agency (JST), Kawaguchi
332-0012, Japan  Yoshinori Tokura  RIKEN Center for Emergent Matter Science,
Wako 351-0198, Japan  Department of Applied Physics and Quantum-Phase
Electronics Center (QPEC), University of Tokyo, Tokyo 113-8656, Japan  Tokyo
College, University of Tokyo, Tokyo 113-8656, Japan  Fumitaka Kagawa
<EMAIL_ADDRESS>RIKEN Center for Emergent Matter Science, Wako
351-0198, Japan  Department of Applied Physics and Quantum-Phase Electronics
Center (QPEC), University of Tokyo, Tokyo 113-8656, Japan  Department of
Physics, Tokyo Institute of Technology, Tokyo 152-8551, Japan

###### Abstract

Phase diagrams form the basis for the study of material science, and the
profiles of phase-transition lines separating different thermodynamic phases
include comprehensive information about thermodynamic quantities, such as
latent heat. However, in some materials exhibiting field-induced first-order
transitions (FOTs), the equilibrium phase-transition line is hidden by the
hysteresis region associated with the FOT; thus, it cannot be directly
determined from measurements of resistivity, magnetization, etc. Here, we
demonstrate a thermodynamics-based method for determining the hidden
equilibrium FOT line. This method is verified for the FOT between
antiferromagnetic and ferrimagnetic states in magneto-electric compounds
(Fe0.95Zn0.05)2Mo3O8. The equilibrium FOT line determined based on the
Clausius–Clapeyron equation exhibits a reasonable profile in terms of the
third law of thermodynamics, and it shows marked differences from the
midpoints of the hysteresis region. Our findings highlight that care should be
taken for referring to the hysteresis midpoint line when discussing field-
induced latent heat or magnetocaloric effects.

Condensed matter often changes its structural/electric/magnetic states in
response to changes in temperature, $T$, and external fields such as magnetic
fields, $H$, electric fields, $E$, and pressures, $P$. These changes are
concisely summarized in phase diagrams. In addition to separating different
phases, the profile of a phase-transition line includes comprehensive
information on thermodynamic quantities associated with the phase transition.
For example, in the $P$–$T$ phase diagram of 3He, the slope of the first-order
transition (FOT), $\frac{dT^{*}}{dP^{*}}$ (the asterisk represents a value on
the phase boundary), between the solid and liquid phases becomes negative
below 0.32 K, indicating that the solid 3He has larger entropy than the liquid
3He, which is a behaviour known as the Pomeranchuk effect 1. Additionally, in
the $P$–$T$ phase diagrams in certain organic conductors, the slope of the FOT
line between antiferromagnetic (AFM) and superconducting (SC) phases is
negative, whereas the slope between quantum spin liquid (QSL) and SC phases is
almost perpendicular to the pressure axis. This observation indicates that the
QSL phase has larger entropy than the AFM phase 2; 3. Furthermore, for the
case of an FOT, the Clausius–Clapeyron equation enables a quantitative
estimation of the entropy change (or latent heat) accompanying the FOT,
$\Delta S^{*}$ (or $T^{*}$$\Delta S^{*}$), if the change in the other
extensive quantity is given (the magnetization change, $\Delta M^{*}$, for the
case of an $H$–$T$ phase diagram). For instance, accurate evaluation of
magnetic-field-induced latent heat has been a central issue in the study of
magnetocaloric effects 4, and thus, comparison with theoretical values derived
from the correct phase diagram provides important insights. Thus, the accurate
determination of a phase-boundary profile can be an important issue in many
contexts.

However, in a real material, an FOT does not necessarily occur at the
equilibrium FOT line but in a parameter region away from it. This behaviour
originates from the fact that an FOT generally involves nucleation and growth,
which are nonequilibrium kinetic processes. As a result, two hysteresis lines,
which we define as the lower-field and higher-field boundaries of the
hysteresis region, are often drawn instead of one equilibrium FOT line,
particularly in an experimentally determined phase diagram. In fact, the
precise determination of the equilibrium FOT line is often not straightforward
in experiments.

We consider two prototypical cases of a field-temperature phase diagram that
exhibits an FOT (Figs. 1a and b). The first case represents an FOT between
“up” and “down” phases connected by a symmetry operation, such as time-
reversal or space-inversion operations (Fig. 1a). Typical ferromagnetic (FM)
and ferroelectric materials are classified into this case if uniform magnetic
and electric fields are chosen as the external fields, respectively. The free
energies of the two states obviously degenerate in the absence of a symmetry-
breaking field 5, and the equilibrium FOT line lies on the zero-field line by
definition. The hysteresis lines in positive- and negative-field regions
should be symmetric regarding the zero-field line, and thus, one may determine
the equilibrium FOT points by taking the midpoints of the two hysteresis
lines. In the second case, by contrast, the competing two phases do not
degenerate at zero field (Fig. 1b). A typical example includes a magnetic-
field-induced FOT from AFM to FM phases. In such a case, there is no symmetry
reason for the midpoints of the two hysteresis lines to agree with the
equilibrium FOT points. To the best of the authors’ knowledge, the validity of
considering the midpoints as the FOT points has never been carefully discussed
for a field-induced FOT.

The uncertainty of a field-induced equilibrium FOT line may be crucial for a
material exhibiting a wide hysteresis; for instance, various materials, such
as Gd5Ge4 6, LaFe12B6 7, doped CeF2 alloys 8; 9, doped Mn2Sb 10; 11, doped
manganese oxides 12; 13; 14, and martensitic materials 15; 16 are known to
exhibit distinct hysteresis broadening at low temperatures. Similar hysteresis
broadening is also observed in a pressure-induced liquid-liquid transition of
an aqueous solution 17. In the present study, by focusing on (Fe1-yZny)2Mo3O8
(Fig. 1c), we exemplify how the equilibrium FOT line is determined for a
material exhibiting a broad hysteresis and show that the midpoints of the two
hysteresis lines appreciably deviate from the equilibrium FOT line. This
finding indicates that considering the midpoints of hysteresis as the
equilibrium FOT line potentially leads to an erroneous conclusion on
thermodynamic quantities, such as field-induced latent heat accompanying the
FOT, even at a qualitative level.

### Introduction to (Fe1-yZny)2Mo3O8

Our target material, (Fe1-yZny)2Mo3O8, is a polar crystal (space group
$P6_{3}mc$) with a linear magnetoelectric effect 18; 19. The magnetic
properties are dominated by two kinds of Fe2+ sites surrounded by oxygen
tetrahedra (A-sites) and octahedra (B-sites) (Fig. 1c); local spins of three
Mo4+ ions form a nonmagnetic spin-trimer singlet 20 and thus have no
contribution to the magnetism. According to a previous neutron study on the
mother compound Fe2Mo3O8 21, the magnetic moments at octahedral sites are
slightly larger than those at tetrahedral sites; hence, the spin configuration
shown in Fig. 1d has no net magnetization (antiferromagnetic; AFM), whereas
that shown in Fig. 1e creates an appreciable macroscopic magnetization
(ferrimagnetic; FRI). The ground state at zero magnetic field is the AFM
phase, which is replaced with the FRI phase when a sufficiently high magnetic
field is applied. For undoped Fe2Mo3O8, the magnetic-field-induced FOT between
the AFM to the FRI phases occurs far above 15 T below 30 K 18.

The molecular-field analysis 21 argues that the intralayer A–B, interlayer A–B
and interlayer A–A magnetic interactions are all AFM, whereas the interlayer
B–B magnetic interaction is weakly FM. In addition, the magnitudes of the
interlayer A–B and interlayer A–A magnetic interactions are comparable,
resulting in a delicate energy balance between the AFM and FRI phases in this
system. Doped nonmagnetic Zn ions selectively occupy tetrahedral sites 22 and
weaken the effective interlayer A–A magnetic interactions. Thus, the Zn-doping
stabilizes the FRI phase, and accordingly the transition field to the FRI
phase decreases upon Zn doping. In this study, we chose 5% Zn-doping
($y=0.05$) so that both the AFM and FRI phases are accessible in the feasible
magnetic field range of 14 T 18.

### Notes on thermodynamic analysis

#### Clausius–Clapeyron equation

On an equilibrium FOT line separating two phases, $\alpha$ and $\beta$, the
thermodynamic potentials of the two phases are equal to each other by
definition. This condition leads to the Clausius–Clapeyron equation, which is
given in the following form for the case of an equilibrium FOT line in an
$H$–$T$ phase diagram:

$\frac{dT^{*}}{dH^{*}}=-\frac{M^{*}_{\rm\alpha}-M^{*}_{\rm\beta}}{S^{*}_{\rm\alpha}-S^{*}_{\rm\beta}}\equiv-\frac{\Delta
M^{*}_{\rm\alpha\beta}}{\Delta S^{*}_{\rm\alpha\beta}},$ (1)

where ($H^{*}$, $T^{*}$) represents an arbitrary point on the equilibrium FOT
line and $M^{*}_{\rm\alpha}$ (or $M^{*}_{\rm\beta}$) and $S^{*}_{\rm\alpha}$
(or $S^{*}_{\rm\beta}$) are the magnetization and entropy values of the
$\alpha$ (or $\beta$) phase on the equilibrium FOT line, respectively. Note
that any equilibrium FOT line in a $H$–$T$ phase diagram must satisfy this
equation. This requirement means that by following the Clausius–Clapeyron
equation, an equilibrium FOT line can be drawn sequentially by integrating
$\frac{dT^{*}}{dH^{*}}$, as detailed in Results section.

To determine the slope of the FOT, $\frac{dT^{*}}{dH^{*}}$, the entropy and
magnetization changes, $\Delta S^{*}_{\rm\alpha\beta}$ and $\Delta
M^{*}_{\rm\alpha\beta}$, respectively, accompanying the magnetic-field-induced
FOT should be given. To this end, the $S$–$H$ and $M$–$H$ curves should be
derived for the $\alpha$ and $\beta$ phases, separately. Although the $M$–$H$
curves are readily measurable, the $S$–$H$ curves are not because the
magnetic-field dependence of the specific heat, $c_{p}$, does not provide the
$S$–$H$ curve. To obtain the $S$–$H$ curves and thus $\Delta
S^{*}_{\rm\alpha\beta}$, we measure the $c_{p}(T)$ at predetermined magnetic
fields: $H$ = 0, 0.1, 1, 2, 3, 4, 5, 6 and 7 T (we use the notation
$H\equiv\mu_{0}H$ as an external magnetic field, where $\mu_{0}$ is magnetic
permeability in vacuum). Thus, the $S_{H}(T)$ are derived by using
$S_{H}(T)=\int_{0}^{T}\frac{c_{p}(H,T^{\prime})}{T^{\prime}}dT^{\prime}$, and
discretized $S_{T}(H)$ data points at a given $T$ are obtained.

#### Application of the Maxwell relation

In principle, if the $S$–$T$ curves are given for many magnetic field points,
isothermal $S$–$H$ curves at different temperatures can be derived with
sufficient accuracy. However, from an experimental point of view, the
measurement of specific heat requires much more time than that of
magnetization; therefore, collecting $c_{p}(T)$ at many magnetic-field points
is not practical. To overcome this technical difficulty, we refer to the
Maxwell relation, $\left(\frac{\partial S}{\partial
H}\right)_{T}=\left(\frac{\partial M}{\partial T}\right)_{H}$, which enables
compensation for the limited number of $S_{T}(H)$ data points by referring to
easily measurable magnetization data. Thus, isothermal $S$–$H$ curves are
obtainable with improved accuracy.

For the purpose of deriving $\Delta S^{*}$ accompanying an equilibrium FOT,
however, the Maxwell relation should be carefully used. To begin, the Maxwell
relation cannot be used at an ideal FOT in principle because its derivation
postulates that $S$ and $M$ are a differentiable continuous function of $H$
and $T$. This prerequisite is obviously not satisfied at an ideal FOT, which
is accompanied by discontinuity in $S$ and $M$. One may assume that an actual
FOT often exhibits relatively continuous change in $S$ and $M$ as a result of
the evolution of the two-phase mixture (see Supplementary Fig. 1), and thus,
the application of the Maxwell relation to the transition region appears fine.
Nevertheless, the integration of $\left(\frac{\partial S}{\partial
H}\right)_{T}=\left(\frac{\partial M}{\partial T}\right)_{H}$ with respect to
$H$ across the transition region still tends to overestimate the intrinsic
value of $\Delta S^{*}$; that is, the integration of the Maxwell relation
across the FOT, frequently used in the analysis of the magnetocaloric effect4;
23, cannot be applied to the precise determination of $\Delta S^{*}$ (for more
details, see Supplementary Note 1). Therefore, we restrict the use of the
Maxwell relation to when the magnetic state is considered a single phase of
either the AFM or FRI phases. $\Delta S^{*}$ cannot be obtained from the
Clausius-Clapeyron equation either 23, because the determination of the exact
profile of the equilibrium FOT line itself is the purpose of this study.

Fig. 1 : Field-temperature phase diagram
of$\mathbf{(Fe_{0.95}Zn_{0.05})_{2}Mo_{3}O_{8}.}$ a,b Phase diagrams with
temperature and external field axes in the system with field-induced first-
order phase transitions: (a) ferroic-order case in which two states are
connected with symmetry operations, such as time-reversal and space-inversion
symmetries, and (b) two-phase-competing case in which two thermodynamic phases
indicated by $\alpha$ and $\beta$ are not degenerate in free energy at zero
external field. The solid lines in (a) and (b) represent magnetic fields at
which the phase transition is observed during magnetic field sweeps at a given
rate. c The crystal structures of (Fe1-yZny)2Mo3O8. Fe2+ ions at the A-sites
(brown) and B-sites (blue) are surrounded by oxygen tetrahedra and octahedra,
respectively. Mo4+ ions (gray) form the nonmagnetic spin-trimer. d,e The
magnetic structure of the (d) antiferromagnetic (AFM) and (e) ferrimagnetic
(FRI) phases. f The isothermal magnetization curves of (Fe0.95Zn0.05)2Mo3O8 at
selected temperatures. Each curve is shifted by 2 $\mu_{\rm B}/{\rm f.u.}$ for
clarity. The transition from the AFM to the FRI phases and from the FRI to the
AFM phases are indicated by red and blue triangles, respectively. g The
temperature-magnetic-field phase diagrams of (Fe0.95Zn0.05)2Mo3O8. The
transitions from the AFM to the FRI phases and from the FRI to the AFM phases
are indicated by red and blue circles, respectively. Gray closed and open
circles represent the transition between PM and AFM and the crossover between
PM to FRI, respectively, and they were determined from $M$–$T$ curves. Green
triangles represent the midpoints between two transitions from the AFM (FRI)
to the FRI (AFM) phases. In g,f, the gray-hatched areas represents the
hysteresis region.

## Results

### Field-temperature phase diagram

Figure 1f shows the isothermal magnetization curves of the target compound,
(Fe0.95Zn0.05)2Mo3O8. The red (blue) triangles indicate the transition fields
from the AFM to the FRI phases (the FRI to the AFM phases) upon increasing
(decreasing) magnetic field. The two hysteresis lines accompanying the
magnetic-field-induced FOT in a positive magnetic field are drawn by tracing
these two transition-fields. The experimental phase diagram is thus obtained,
as shown in Fig. 1g.

This material exhibits distinct hysteresis broadening at low temperatures. The
hysteresis width, which is the difference between two transition fields, is as
small as 0.4 T at 40 K; however, it significantly increases to $\approx$ 20 T
at 6.5 K. In Fig. 1g, the midpoints of the two hysteresis lines are plotted
for reference. Note that the midpoint line starts to bend below 18 K, and its
slope is appreciably positive at 6.5 K. If one considers the midpoint line as
an approximate equilibrium FOT line, the midpoint line should bend again at
low temperatures to be perpendicular to the magnetic-field axis to satisfy the
third law of thermodynamics (Nernst–Planck hypothesis 24). Thus, the
Clausius–Clapeyron equation concludes that $\Delta S^{*}$ ($\equiv S^{*}_{\rm
FRI}-S^{*}_{\rm AFM}$) and/or $\Delta M^{*}$ ($\equiv M^{*}_{\rm
FRI}-M^{*}_{\rm AFM}$) should exhibit complicated behaviour below 18 K. This
finding, however, appears unusual, and its validity should be carefully
considered by testing whether the midpoint line accurately represents the
equilibrium FOT line.

###  $\bm{M}$–$\bm{T}$ curves

We aimed to determine the equilibrium FOT line by following the equilibrium
thermodynamics for a single phase. It was therefore important to avoid the
data analysis in the $(H,T)$ region where the magnetic state is considered
AFM–FRI phase mixture. In particular, extensive attention has to be paid to
the magnetic state during cooling because we collected the $c_{p}(T)$ upon
decreasing temperature. To identify the $(H,T)$ region of the two-phase
mixture during cooling, we measured the $M$–$T$ curves at various magnetic
fields during field-cooling (FC) and field-warming (FW) processes; in FC
process, the magnetic fields were applied at 100 K ($>T_{\rm c}\approx$ 56 K)
and the magnetization was measured from 100 to 2 K while retaining the
magnetic fields; then, in FW process, the magnetization was measured from 2 to
100 K at the same magnetic fields. The data at selected magnetic fields are
shown in Fig. 2a. Overall, a transition from the PM to the AFM phase was
observed below 1.5 T, whereas that to the FRI phase was observed above 3 T. At
an intermediate field ranging from 1.6 to 3 T, a sharp, continuous change in
magnetization and the associated thermal hysteresis were observed, signifying
the thermally induced FOT between the AFM and FRI phases (see also Fig. 1g).

### Identification of the $\bm{(H,T)}$ region of single magnetic phase during
field cooling

Fig. 2 : Identification of single-phase AFM and FRI phases from
$\bm{M}$-$\bm{T}$ curves. a Temperature dependencies of the magnetization upon
field-cooling (FC; solid lines) and field-warming (FW; dashed lines)
processes. b Schematic phase diagram highlighting the single thermodynamic AFM
and FRI phases. The black dots represent data points used for the present
analysis while carefully avoiding the “two-phase mixture” (gray) region. Fig.
3 : Temperature dependences of the specific heats and entropies under magnetic
fields. a Temperature dependencies of the specific heats under various
magnetic fields. The specific heat measurement was performed during the field-
cooling process. Inset: Enlarged view around the transitions. The $c_{p}(T)$
at 0 and 0.1 T are almost identical and no difference is distinguishable on
this scale. b–j The temperature dependencies of the entropies under various
magnetic fields.

From the $M$–$T$ curves, we determined the $(H,T)$ regions where the two-phase
coexistence could occur during FC, as summarized in Fig. 2b, according to the
following two criteria. First, in the temperature range where the transition
progressed, we determined that the data were affected by the two-phase
coexistence. Second, especially below 30 K, even if the transition did not
appear to progress, we determined that the two-phase mixture was present when
its magnetization value exhibited an intermediate value between the AFM and
FRI single phases. The lowest-temperature state below 2 T consisted
exclusively of the AFM single phase (see Supplementary Note 2 and
Supplementary Fig. 2, in which the $M$–$T$ curve under FW at 2 T after zero FC
was compared with that under FC at 2 T). Similarly, above 3 T, the
magnetization at the lowest temperature was insensitive to the magnetic field
(see also Fig. 1f), indicating the FRI single phase. In contrast, at 2.25 T,
although thermal hysteresis was not observed below 25 K, the magnetization
value at the lowest temperature was intermediate between the AFM and FRI
single phase values. Thus, the magnetic state at 2.25 T below 25 K was
considered the two-phase mixture.

The $(H,T)$ region that fell under “two-phase mixture” criteria is shown in
grey in Fig. 2b. The data from this region were not used in the following
single-phase thermodynamic analysis to avoid a possible influence of the two-
phase mixture.

Fig. 4 : Derivation of isothermal $\bm{\left(\frac{\partial S}{\partial
H}\right)_{T}}$ curves via the Maxwell relation. a $\frac{dM}{dT}$–$T$ curves
measured upon decreasing temperature. The data indicated by the dotted lines
were not used for the Maxwell relation. b,c Isothermal $\frac{dM}{dT}$–$H$
curves of the (b) AFM and (c) FRI phases. The solid lines represent the best-
fit results (see the main text). The legend shown in (b) also applies to (c).
Fig. 5 : Determination of $\bm{\Delta S^{*}}$ and $\bm{\Delta M^{*}}$. a
Isothermal entropy curves at selected temperatures. Closed circles are the
data obtained from the specific heat measurements. Blue and red solid lines
indicate the fittings for the AFM and the FRI phases, respectively, based on
the $M$–$T$ curve and the Maxwell relation. The error bars are estimated from
the least-squares fitting. b Isothermal magnetization curves at selected
temperatures. Each magnetization curve is shifted by 10 J/(mol$\cdot$T) for
visibility. c Temperature dependence of the entropy and magnetization changes
on the determined phase boundary. d Temperature dependence of the ratio of
$\Delta S^{*}$ and $\Delta M^{*}$. The error bars in (c) and (d) are estimated
from the standard errors of the parameters when performing the least-squares
fitting of the cubic function to the $\frac{dM}{dT}$ data.

### Derivation of $\bm{S}$–$\bm{H}$ curves

Figure 3a shows the temperature dependences of the specific heat, $c_{p}$, at
various magnetic fields measured upon field cooling. A clear lambda-like
anomaly is observed in the $c_{p}(T)$ at approximately 55 K at 0 T, signifying
a second-order transition from the paramagnetic (PM) to the AFM phases. The
application of a magnetic field larger than 2 T suppresses the $\lambda$ peak,
reflecting that a change from the PM to the FRI phases is crossover under a
finite magnetic field (see Fig. 1g). Similar behaviours have been reported in
previous studies 19; 25; 26. The $S$–$T$ profiles are obtained by integrating
$c_{p}/T$ with respect to $T$. As described in Supplementary Information
(Supplementary Note 3), the value of entropy at 2 K is obtained at each
magnetic field by extrapolating the $c_{p}(T)$ to zero temperature with
$T^{3}$ behaviour. The resulting $S(T)$ data are shown in Figs. 3b–j. Note
that a discontinuous entropy change associated with an FOT cannot be
accurately evaluated by integrating $c_{p}/T$ because
$\frac{dS}{dT}=\frac{c_{p}}{T}$ is valid only when $dS/dT$ is well defined
(i.e. $S$ is a differential continuous function of $T$). If the latent heat
associated with the FOT can be accurately captured by specific heat
measurements, the entropy above the transition temperature can be determined;
however, this is not the case. Therefore, the $S$–$T$ curve (i.e., an entropy
counting from zero temperature) at 2 T is only valid up to 27 K (Fig. 3e).
From these $S_{H}(T)$ data set, the discretized $S_{T}(H)$ points are derived
(for instance, see Fig. 5a). The comparison of $S_{H}(T)$ between 0 and 7 T
tells that the FRI phase has a higher entropy than the AFM phase
(Supplementary Fig. 4).

Then, we aim to derive the $S$–$H$ curves with accuracy from the limited
number of available $S_{T}(H)$ data points by using the Maxwell relation. We
derive the $\left(\frac{\partial M}{\partial T}\right)_{H}-T$ curves from the
$M$–$T$ curves arising from a single-phase region (either the AFM or FRI
phases). The results at selected magnetic fields are shown in Fig. 4a with
solid lines (for the data above 50 K, see Supplementary Note 4). Then, from
the continuous $\left(\frac{\partial M}{\partial T}\right)_{H}-T$ curves, we
obtain the discretized $\left(\frac{\partial M}{\partial T}\right)_{H}$ data
points as functions of $H$ at selected temperatures for the AFM and FRI single
phases, as shown in Figs. 4b and 4c, respectively.

To obtain the isothermal $S$–$H$ curves, $\left(\frac{\partial S}{\partial
H}\right)_{T}$ should be integrated with respect to $H$. To this end, we first
fit the magnetic field dependence of $\left(\frac{\partial M}{\partial
T}\right)_{H}=\left(\frac{\partial S}{\partial H}\right)_{T}$ for the AFM and
FRI phases separately by the following polynomial functions:

$\left(\frac{\partial M_{\rm AFM}(H,T)}{\partial
T}\right)_{H}=\left(\frac{\partial S_{\rm AFM}(H,T)}{\partial
H}\right)_{T}=a_{3}(T)H^{2}+a_{2}(T)H+a_{1}(T)$ (2) $\left(\frac{\partial
M_{\rm FRI}(H,T)}{\partial T}\right)_{H}=\left(\frac{\partial S_{\rm
FRI}(H,T)}{\partial H}\right)_{T}=b_{3}(T)(H-7)^{2}+b_{2}(T)(H-7)+b_{1}(T).$
(3)

The fittings are successful as shown in Figs. 4b and 4c, from which the
parameters $a_{i}$ and $b_{i}$ ($i=1,2,3$) are determined. The $S$–$H$ curves
are thus given as follows:

$S_{\rm AFM}(H,T)=S_{\rm AFM}(0,T)+\int_{0}^{H}\left(\frac{\partial S_{\rm
AFM}}{\partial
H^{\prime}}\right)_{T}dH^{\prime}=\frac{a_{3}(T)}{3}H^{3}+\frac{a_{2}(T)}{2}H^{2}+a_{1}(T)\mu_{0}H+a_{0}(T),$
(4) $S_{\rm FRI}(H,T)=S_{\rm FRI}(7,T)+\int_{7}^{H}\left(\frac{\partial S_{\rm
FRI}}{\partial
H^{\prime}}\right)_{T}dH^{\prime}=\frac{b_{3}(T)}{3}(H-7)^{3}+\frac{b_{2}(T)}{2}(H-7)^{2}+b_{1}(T)(H-7)+b_{0}(T).$
(5)

Note that $a_{i}$ and $b_{i}$($i$=1,2,3) are given from the fitting of the
$\left(\frac{\partial M}{\partial T}\right)_{H}$ data points. Finally, to
determine the remaining parameters $a_{0}$ (for the AFM phase) and $b_{0}$
(for the FRI phase), we refer to the discretized $S_{T}(H)$ points, which are
independently obtained from the $c_{p}$–$T$ curves. Because there is only one
unknown parameter for each phase, the problem of the limited $S_{T}(H)$ points
is minimized. As shown in Fig. 5a, the fitting $S$–$H$ curves based on the
magnetization data, equations (4) and (5), well reproduce the $S_{T}(H)$ data
points, which are obtained from the specific heat measurement. This good
agreement corroborates that the present method based on the single-phase
thermodynamics successfully avoids the experimental artifacts in the derived
$S$–$H$ curve27; 28; 29 that are often caused by a hysteretic behaviour and a
kinetic effect associated with a FOT (see Supplementary Note 1). In this way,
the continuous $S$–$H$ curves are derived at various temperatures for the AFM
and FRI single phases, separately. The resulting two curves are not connected
to each other, representing the entropy discontinuity accompanying the FOT.
Overall, the entropy difference decreases as the temperature decreases.
Nevertheless, to determine $\Delta S^{*}$ more precisely, the equilibrium
transition field, $H^{*}$, should be determined at each temperature. This
issue is discussed below.

###  $\bm{M}$–$\bm{H}$ curves

To apply the Clausius–Clapeyron equation, $\Delta M^{*}$ should also be
derived from the isothermal $M$–$H$ curves. In contrast to the $S$–$H$ curves,
the $M$–$H$ curves are straightforwardly measured, as shown in Fig. 5b. To
obtain the magnetization of a single phase, we first prepare the AFM and FRI
single phases by ZFC and 7 T-FC from 100 K ($>T_{\rm c}\approx 56$ K),
respectively; then, the magnetization is measured while increasing and
decreasing the magnetic fields, respectively; the cooling process from 100 K
is performed after each $M$–$H$ measurement. The two branches are separated
from each other, and the transition to the other branch occurs when the
magnetic field reaches a value of the hysteresis line (Figs. 1f,g). From these
$M$–$H$ curves, one can determine $\Delta M^{*}$, if the equilibrium
transition field $H^{*}$ is given at each temperature.

### Derivation of the equilibrium first-order transition line

The determination of $\Delta S^{*}$ and $\Delta M^{*}$ and that of the
equilibrium transition field $H^{*}$ are an intertwined issue. If $H^{*}$ is
given at a certain temperature, $\Delta S^{*}$ and $\Delta M^{*}$ at the same
temperature are determined; if $\Delta S^{*}$ and $\Delta M^{*}$ are given at
a certain temperature, $H^{*}$ at a nearby temperature is determined. To
perform this sequential determination, the starting point for drawing the
equilibrium FOT line must be determined first.

The starting point is chosen at ($H^{*}_{0}$, $T^{*}_{0}$) = (1.81 T, 45 K)
because the hysteresis width is as small as 0.1 T and it is reasonable to
approximate the midpoint of the two hysteresis lines as the equilibrium FOT
point. First, the values of $\Delta S^{*}$ and $\Delta M^{*}$ at (1.81 T, 45
K) are determined by referring to the isothermal $S$–$H$ and $M$–$H$ curves,
respectively; the slope of the FOT line at ($H^{*}_{0}$, $T^{*}_{0}$),
$\frac{dH^{*}}{dT^{*}}|_{T^{*}_{0}}$, is calculated from $-\frac{\Delta
S^{*}}{\Delta M^{*}}$. The next transition field $H^{*}_{1}$ at the nearby
temperature $T^{*}_{1}$ ($<T^{*}_{0}$) is calculated as follows:

$\displaystyle H^{*}_{1}$
$\displaystyle=H^{*}_{0}+\left(\frac{dH^{*}}{dT^{*}}\right)_{T^{*}_{0}}\times(T^{*}_{1}-T^{*}_{0})$
$\displaystyle=H^{*}_{0}-\left(\frac{\Delta S^{*}}{\Delta
M^{*}}\right)_{T^{*}_{0},H^{*}_{0}}\times(T^{*}_{1}-T^{*}_{0})$ (6)

We repeat this procedure sequentially, and the transition field $H^{*}_{n}$ at
$T^{*}_{n}$ ($n\geq 1$) is obtained by the following equations:

$\displaystyle\Delta S^{*}(H^{*}_{n},T^{*}_{n})$ $\displaystyle=S_{\rm
FRI}(H^{*}_{n},T^{*}_{n})-S_{\rm AFM}(H^{*}_{n},T^{*}_{n})$ (7)
$\displaystyle\Delta M^{*}(H^{*}_{n},T^{*}_{n})$ $\displaystyle=M_{\rm
FRI}(H^{*}_{n},T^{*}_{n})-M_{\rm AFM}(H^{*}_{n},T^{*}_{n})$ (8) $\displaystyle
H^{*}_{n+1}$
$\displaystyle=H^{*}_{n}+\left(\frac{dH^{*}}{dT^{*}}\right)_{T^{*}_{n}}\times(T^{*}_{n+1}-T^{*}_{n})$
$\displaystyle=H^{*}_{n}-\left(\frac{\Delta S^{*}}{\Delta
M^{*}}\right)_{T^{*}_{n},H^{*}_{n}}\times(T^{*}_{n+1}-T^{*}_{n})$ (9)

The decrement of the temperature step, $T^{*}_{n+1}-T^{*}_{n}$, is determined
by the temperatures at which the isothermal $M$–$H$ curves are measured, and
it is $-$2.5 K in the present study. Thus, $\Delta S^{*}$ and $\Delta M^{*}$ ,
$-\frac{\Delta S^{*}}{\Delta M^{*}}$, and the equilibrium FOT line are
determined sequentially for temperatures below 45 K, as displayed in Figs.
5c,d and Fig. 6, respectively. $\Delta S^{*}$ is sensitive to temperature and
monotonically decreases to zero as the temperature decreases to zero
temperature, whereas $\Delta M^{*}$ depends only weakly on temperature (Fig.
5c). Thus, the ratio, $-\frac{\Delta S^{*}}{\Delta M^{*}}$, monotonically
approaches to zero at zero temperature (Fig. 5d).

The obtained equilibrium FOT line (Fig. 6) becomes appreciably distinct from
the midpoints of the two hysteresis lines, especially below 18 K. Note that
the hysteresis distinctly broadens below 18 K. Toward zero temperature, the
obtained equilibrium FOT line monotonically becomes perpendicular to the
magnetic-field axis to satisfy the third law of thermodynamics. The time
evolution measurements of the net magnetization, $M(t)$, indicates that the
equilibrium phase at 2.3 T and 10 K is the AFM phase, rather than the FRI
phase (Supplementary Note 5 and Supplementary Fig. 6), further confirming that
the profile of the obtained equilibrium FOT line is consistent with the
observation. The asymptotic behaviour of the FOT line indicates that nothing
peculiar is involved below 18 K, as demonstrated in Figs. 5c and 5d, in
contrast to the implications drawn from the profile of the midpoint line.

Fig. 6 : Determination of the equilibrium first-order transition line.
Determined phase boundary indicated by orange circles with error bars, which
are superimposed on the low-temperature enlarged view of Fig. 1g.

## Discussion

Finally, we show that if the midpoint line is regarded as the equilibrium FOT
line, it can lead to even qualitatively incorrect conclusion regarding
thermodynamic quantities such as latent heat. By using equation (1), the
latent heat $\Delta q^{*}$ accompanying the FOT from the AFM to the FRI can be
calculated as $\Delta q^{*}=T^{*}\Delta S^{*}=-T^{*}\Delta
M^{*}\left(\frac{dH^{*}}{dT^{*}}\right)$ for either the midpoint or
equilibrium FOT lines. The calculated values of the latent heat for each case
are compared in Fig. 7. At high temperatures above 25 K, $\Delta q^{*}$ is
positive, and its values derived from the midpoint and equilibrium FOT lines
agree well with each other. However, the two curves show distinct behaviour
below 25 K: $\Delta q^{*}$ derived from the equilibrium FOT line approaches
zero asymptotically toward zero temperature, whereas that derived from the
midpoint line shows nonmonotonous temperature dependence and even a large
negative value, $\approx-$14 J/mol at 10 K. Thus, the adoption of the midpoint
line as the equilibrium FOT line causes the erroneous latent heat even at a
qualitative level, especially at low temperatures, demonstrating that caution
should be exercised when referring to the midpoint lines.

The experimental observation of FOTs is known to be susceptible to kinetic
effects, which are beyond the framework of the equilibrium thermodynamics. We
have demonstrated that the equilibrium FOT line can nevertheless be determined
using the single-phase thermodynamics with avoiding analysis of the phase-
mixed state. The understanding of the temperature-dependent
agreement/disagreement between the midpoint and equilibrium FOT lines remains
an open question in this study. The large deviation of the midpoint line from
the equilibrium line at low temperatures appears to be related to distinct
hysteresis broadening. This observation may suggest that the details of the
phase evolution in the FOT may vary with temperature.

Fig. 7 : The comparison of calculated latent heats between midpoint and
equilibrium first-order-transition lines. Temperature dependences of the
calculated latent heats based on the midpoint line (red triangles) and the
equilibrium FOT line (blue circles).

## Data availability

The datasets used and/or analysed during the current study available from the
corresponding author on reasonable request.

## Methods

### Sample preparation

A single crystal of (Fe0.95Zn0.05)2Mo3O8 was grown by the chemical-vaper
transport method. The details about the growth of this sample were reported in
the literature 30; 31; 18. All magnetization and specific heat measurements
were performed using the same sample with a relatively small mass ($\approx$ 3
mg).

### Magnetization measurement

The magnetization along the $c$–axis was measured in Quantum Design MPMS-XL,
MPMS-3, and PPMS 14 T.

### Specific heat measurement

Specific heat measurements were performed in a field-cooling procedure, and
the standard relaxation method available in the heat capacity option of
Quantum Design PPMS 14 T was used. A small amount of grease (Apiezon N) was
used to ensure good thermal contact between the sample stage (platform) and
the sample.

## References

  * (1) Richardson, R. C.  The Pomeranchuk effect.  Rev. Mod. Phys. 69, 683–690 (1997).
  * (2) Pustogow, A., Bories, M., Löhle, A., Rösslhuber, R., Zhukova, E., Gorshunov, B., Tomić, S., Schlueter, J. A., Hübner, R., Hiramatsu, T., Yoshida, Y., Saito, G., Kato, R., Lee, T.-H., Dobrosavljević, V., Fratini, S. & Dressel, M.  Quantum spin liquids unveil the genuine Mott state.  Nat. Mater. 17, 773–777 (2018).
  * (3) Furukawa, T., Kobashi, K., Kurosaki, Y., Miyagawa, K. & Kanoda, K.  Quasi-continuous transition from a Fermi liquid to a spin liquid in $\kappa$-(ET)2Cu2(CN)3.  Nat. Commun. 9, 1–7 (2018).
  * (4) Franco, V., Blázquez, J., Ingale, B. & Conde, A.  The magnetocaloric effect and magnetic refrigeration near room temperature: materials and models.  Annu. Rev. Mater. Res. 42, 305–342 (2012).
  * (5) Chaikin, P. M., Lubensky, T. C. & Witten, T. A.  Principles of condensed matter physics (Cambridge university press, Cambridge, 1995).
  * (6) Roy, S., Chattopadhyay, M., Banerjee, A., Chaddah, P., Moore, J., Perkins, G., Cohen, L., Gschneidner Jr, K. & Pecharsky, V.  Devitrification of the low temperature magnetic-glass state in Gd5Ge4.  Phys. Rev. B 75, 184410 (2007).
  * (7) Fujieda, S., Fukamichi, K. & Suzuki, S.  Itinerant-electron metamagnetic transition in LaFe12B6.  J. Magn. Magn. Mater. 421, 403–408 (2017).
  * (8) Manekar, M., Chaudhary, S., Chattopadhyay, M., Singh, K., Roy, S. & Chaddah, P.  First-order transition from antiferromagnetism to ferromagnetism in Ce(Fe0.96Al0.04)2.  Phys. Rev. B 64, 104416 (2001).
  * (9) Kumar, K., Pramanik, A., Banerjee, A., Chaddah, P., Roy, S., Park, S., Zhang, C. & Cheong, S.-W.  Relating supercooling and glass-like arrest of kinetics for phase separated systems: Doped CeFe2 and (La, Pr, Ca)MnO3.  Phys. Rev. B 73, 184435 (2006).
  * (10) Kushwaha, P., Rawat, R. & Chaddah, P.  Metastability in the ferrimagnetic–antiferromagnetic phase transition in Co substituted Mn2Sb.  J. Phys. Condens. Matter. 20, 022204 (2008).
  * (11) Singh, V., Karmakar, S., Rawat, R. & Kushwaha, P.  Giant negative magnetoresistance and kinetic arrest of first-order ferrimagnetic-antiferromagnetic transition in Ge doped Mn2Sb.  J. Appl. Phys. 125, 233906 (2019).
  * (12) Tokura, Y.  Critical features of colossal magnetoresistive manganites.  Rep. Prog. Phys. 69, 797 (2006).
  * (13) Rawat, R., Mukherjee, K., Kumar, K., Banerjee, A. & Chaddah, P.  Anomalous first-order transition in Nd0.5Sr0.5MnO3: an interplay between kinetic arrest and thermodynamic transitions.  J. Phys. Condens. Matter. 19, 256211 (2007).
  * (14) Matsuura, K., Oike, H., Kocsis, V., Sato, T., Tomioka, Y., Kaneko, Y., Nakamura, M., Taguchi, Y., Kawasaki, M., Tokura, Y. & Kagawa, F.  Kinetic pathway facilitated by a phase competition to achieve a metastable electronic phase.  Phys. Rev. B 103, L041106 (2021).
  * (15) Ito, W., Ito, K., Umetsu, R. Y., Kainuma, R., Koyama, K., Watanabe, K., Fujita, A., Oikawa, K., Ishida, K. & Kanomata, T.  Kinetic arrest of martensitic transformation in the NiCoMnIn metamagnetic shape memory alloy.  APL 92, 021908 (2008).
  * (16) Niitsu, K., Date, H. & Kainuma, R.  Thermal activation of stress-induced martensitic transformation in Ni-rich Ti-Ni alloys.  Scr. Mater. 186, 263–267 (2020).
  * (17) Suzuki, Y.  Direct observation of reversible liquid–liquid transition in a trehalose aqueous solution.  Proc. Natl. Acad. Sci. USA 119, e2113411119 (2022).
  * (18) Kurumaji, T., Ishiwata, S. & Tokura, Y.  Doping-tunable ferrimagnetic phase with large linear magnetoelectric effect in a polar magnet Fe2Mo3O8.  Phys. Rev. X 5, 031034 (2015).
  * (19) Wang, Y., Pascut, G. L., Gao, B., Tyson, T. A., Haule, K., Kiryukhin, V. & Cheong, S.-W.  Unveiling hidden ferrimagnetism and giant magnetoelectricity in polar magnet Fe2Mo3O8.  Sci. Rep. 5, 1–8 (2015).
  * (20) Cotton, F. A.  Metal atom clusters in oxide systems.  Inorg. Chem. 3, 1217–1220 (1964).
  * (21) Bertrand, D. & Kerner-Czeskleba, H.  Structural and magnetic study of iron-group molybdates.  J. Phys.(Paris) 36, 379–390 (1975).
  * (22) Varret, F., Czeskleba, H., Hartmann-Boutron, F. & Imbert, P.  Étude par effet Mössbauer de l’ion Fe2+ en symétrie trigonale dans les composés du type (Fe, M)2Mo3O8 (M= Mg, Zn, Mn, Co, Ni) et propriétés magnétiques de (Fe, Zn)2Mo3O8.  J. Phys. 33, 549–564 (1972).
  * (23) Xu, K., Li, Z., Zhang, Y.-L. & Jing, C.  An indirect approach based on Clausius–Clapeyron equation to determine entropy change for the first-order magnetocaloric materials.  Phys. Lett. A 379, 3149–3154 (2015).
  * (24) Callen, H. B.  Thermodynamics and an introduction to thermostatistics; 2nd ed. (Wiley, New York, 1985).
  * (25) Ideue, T., Kurumaji, T., Ishiwata, S. & Tokura, Y.  Giant thermal Hall effect in multiferroics.  Nat. Mater. 16, 797–802 (2017).
  * (26) Csizi, B., Reschke, S., Strinić, A., Prodan, L., Tsurkan, V., Kezsmarki, I. & Deisenhofer, J.  Magnetic and vibronic terahertz excitations in Zn-doped Fe2Mo3O8.  Phys. Rev. B 102, 174407 (2020).
  * (27) Liu, G., Sun, J., Shen, J., Gao, B., Zhang, H., Hu, F. & Shen, B.  Determination of the entropy changes in the compounds with a first-order magnetic transition.  APL 90, 032507 (2007).
  * (28) Tocado, L., Palacios, E. & Burriel, R.  Entropy determinations and magnetocaloric parameters in systems with first-order transitions: Study of MnAs.  J. Appl. Phys. 105, 093918 (2009).
  * (29) Amaral, J. & Amaral, V.  On estimating the magnetocaloric effect from magnetization measurements.  J. Magn. Magn. Mater. 322, 1552–1557 (2010).
  * (30) Strobel, P., Le Page, Y. & McAlister, S.  Growth and physical properties of single crystals of Fe${}^{\rm II}_{2}$Mo${}^{\rm IV}_{3}$O8.  J. Solid State Chem. 42, 242–250 (1982).
  * (31) Strobel, P. & Le Page, Y.  Growth and morphology of single crystals of hexagonal molybdates (IV) M2Mo3O8 (M= Mn, Fe, Co, Ni).  J. Cryst. Growth 61, 329–338 (1983).
  * (32) Momma, K. & Izumi, F.  VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data.  J. Appl. Crystallogr. 44, 1272–1276 (2011).

## Acknowledgments

The authors acknowledge enlightening discussions with Y. Taguchi and S. Imajo.
This work was supported by JSPS KAKENHI (Grant No. 21K14398, No. 21H04442, and
No. 18H05225) and JST CREST (No. JPMJCR1874). K.M. was supported by the
Special Postdoctoral Researcher Program of RIKEN. The crystal structure was
visualized by VESTA 3 32.

## Author contributions

K.M. and Y.N. performed the magnetization measurement and analyzed the data.
K.M. and Y.N. performed the specific heat measurement and analyzed the data,
with the help of M.K. T.K. grew the single crystal used for this study. K.M.
and F.K. wrote the manuscript. All authors have discussed the results and
commented on the manuscript.

## Competing interests

The authors declare no competing interests.


# Young Star Clusters Dominate the Production of Detached Black Hole-Star
Binaries

Ugo Niccolò Di Carlo<EMAIL_ADDRESS>Scuola Internazionale Superiore di
Studi Avanzati (SISSA) Via Bonomea 265, I-34136 Trieste, Italy  Department of
Physics and Astronomy, University of North Carolina at Chapel Hill, 120 E.
Cameron Ave, Chapel Hill, NC, 27599, USA  McWilliams Center for Cosmology,
Department of Physics, Carnegie Mellon University, 5000 Forbes Avenue,
Pittsburgh, PA 15213  Poojan Agrawal Department of Physics and Astronomy,
University of North Carolina at Chapel Hill, 120 E. Cameron Ave, Chapel Hill,
NC, 27599, USA  McWilliams Center for Cosmology, Department of Physics,
Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213  Carl L.
Rodriguez Department of Physics and Astronomy, University of North Carolina at
Chapel Hill, 120 E. Cameron Ave, Chapel Hill, NC, 27599, USA  McWilliams
Center for Cosmology, Department of Physics, Carnegie Mellon University, 5000
Forbes Avenue, Pittsburgh, PA 15213  Katelyn Breivik Center for Computational
Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY, 10010, USA

###### Abstract

The recent discovery of two detached black hole-star (BH-star) binaries from
Gaia’s third data release has sparkled interest in understanding the formation
mechanisms of these systems. We investigate the formation of these systems by
dynamical processes in young open star clusters (SCs) and via isolated binary
(IB) evolution, using a combination of direct $N$-body models and population
synthesis simulations. By comparing dynamical and isolated systems created
using the same model of binary stellar evolution, we find that dynamical
formation in SCs is nearly 40 times as efficient per unit of star formation at
producing BH-star binaries compared to IB evolution. We expand this analysis
to the full Milky Way (MW) using a FIRE-2 hydrodynamical simulation of a MW-
mass galaxy. Even assuming that only $10\%$ of star formation goes into SCs,
we find that approximately 4 out of every 5 BH-star systems are formed
dynamically, and that the MW contains a total of $\sim 2\times 10^{5}$ BH-star
systems. Many of these dynamically-formed systems have larger orbital periods,
eccentricities, and black hole masses than their isolated counterparts. For
binaries older than 100 Myr, we show that any detectable system with $e\gtrsim
0.5$ or $M_{\rm BH}\gtrsim 10\,\mathrm{M_{\odot}}$ can _only_ be formed
through dynamical processes. Our MW model predicts between 61 and 210 such
detections from the complete DR4 Gaia catalog, with the majority of systems
being dynamically formed in massive and metal-rich SCs. Finally, we compare
our populations to the recently discovered Gaia BH1 and Gaia BH2, and conclude
that the dynamical scenario is the most favorable formation pathway for both
systems.

##  1 Introduction

Current estimates anticipate that the Milky Way (MW) harbors a population of
$\sim 10^{7}-10^{9}$ stellar black holes (BHs) and $\sim 10^{4}$ massive stars
which are likely BH progenitors (Garmany et al., 1982; Reed, 2003). While the
exact binary fraction of stellar BHs is unknown, the high binary fraction of
BH-progenitor stars (e.g., Sana et al., 2012) suggests that a substantial
quantity of their BH descendants likely exist in binaries as well. Most of the
current observational evidence for BH binary systems in the MW comes from
X-ray binary systems, consisting of 20 dynamically confirmed BHs in X-ray
binaries, with an additional $\sim 50$ X-ray sources identified as strong
candidates for containing BHs (McClintock et al., 2006; Remillard &
McClintock, 2006; Corral & Santana, 2016). However, population synthesis
models indicate that the majority of BHs in binary systems within the MW are
likely to be dormant BHs, with large orbital periods that preclude their
involvement in X-ray binaries (Portegies Zwart et al., 1997; Corral & Santana,
2016; Breivik et al., 2017; Chawla et al., 2022), making them harder to be
detected. Over the past few years, it has become possible to identify dormant
BHs through the motion of their luminous companions, with the first of these
binaries being identified through the radial velocity of their luminous
components in the globular cluster (GC) NGC 3201 (Giesers et al., 2018, 2019).
More recently, proper motion data from ESO Gaia’s Data Release 3 (DR3, Gaia
Collaboration et al., 2022) has enabled the discovery of two dormant BHs in
BH-star binary systems in the galactic field (Chakrabarti et al., 2022; El-
Badry et al., 2023a, b). Combined with follow-up radial velocity measurements
to fully characterize the orbit (El-Badry et al., 2023a, b), these systems
offer new insights into the formation, evolution, and characteristics of BHs
in our galaxy. The first system, Gaia BH1, consists of a Sun-like main
sequence (MS) star orbiting a black hole with mass
$M_{\mathrm{BH}}=9.62\,\mathrm{M_{\odot}}$, eccentricity $e=0.45$ and orbital
period $P=185.6\,$days (El-Badry et al., 2023b). Gaia BH2 is a $\sim
1\,\mathrm{M_{\odot}}$ red giant orbiting a a black hole with mass
$M_{\mathrm{BH}}=8.9\,\mathrm{M_{\odot}}$, eccentricity $e=0.52$ and orbital
period $P=1277\,$days (El-Badry et al., 2023a). According to El-Badry et al.
(2023b) and El-Badry et al. (2023a), the formation histories of both Gaia BH1
and Gaia BH2 are extremely difficult to explain by isolated binary (IB)
evolution, due to the incompatibility between their orbital properties and the
outcomes predicted for common envelope evolution. Star clusters (SCs), on the
other hand, are highly favored environments for the formation of such systems;
not only can dynamical interactions significantly alter the orbital properties
of primordial binaries allowing them to follow other evolutionary pathways,
but BH-star binaries and their progenitors may also dynamically assemble
through binary exchanges and/or multi-body encounters. Despite this, the thin-
disk-like Galactic orbits and near-solar metallicities of the Gaia BHs suggest
that they most likely did not form dynamically in an old, metal-poor GC.
Young/open SCs, on the other hand, are dynamical systems with Galactic orbits
and metallicities perfectly compatible with those of Gaia BHs. Additionally,
they are the birthplace of the vast majority of massive stars, which are the
progenitors of compact objects (e.g. Lada & Lada, 2003a; Portegies Zwart et
al., 2010). Hence, the majority of BHs in the MW have probably spent the first
part of their life in young/open SCs undergoing dynamical interactions.

In this paper, we investigate the dynamical formation of binary systems
consisting of a star and a BH in young/open SCs in the MW and via IB
evolution. We use a large set of $3\times 10^{3}$ $N$-body simulations of
young/open SCs (from Di Carlo et al., 2020b), along with population synthesis
simulations of IB evolution. We describe the population of BH-star binaries in
the MW, predicting the number and properties of the intrinsic population, the
expected detections in Gaia DR3, and in future Gaia data releases. Finally, we
compare our systems with Gaia BH1 and Gaia BH2 to understand their possible
formation pathway, and conclude that both systems (and especially BH2) are
more compatible with a dynamical formation scenario. Throughout this paper, we
use “star” to refer to every luminous stellar type except for white dwarfs and
neutron stars (further discriminating between MS stars and giants when
appropriate).

##  2 Methods

###  2.1 Population of Star Clusters

The star cluster simulations used in this paper were performed with the same
code and methodology described in Di Carlo et al. (2019) and Di Carlo et al.
(2020b). Dynamics is treated by the direct summation $N$-Body code nbody6++gpu
(Wang et al., 2015), coupled with the population synthesis code mobse (Mapelli
et al., 2017; Giacobbo et al., 2018), an upgraded version of BSE (Hurley et
al., 2002). mobse includes prescriptions for core-collapse supernovae (Fryer
et al., 2012), winds of massive stars (Giacobbo & Mapelli, 2018), electron-
capture supernovae (Giacobbo & Mapelli, 2019), natal kicks with fallback
(Giacobbo & Mapelli, 2020) and (pulsational) pair instability supernovae
(Mapelli et al., 2020). Orbital evolution due to gravitational-wave emission
is treated with the equations from Peters (1964).

In this work, we analyze $3\times 10^{3}$ simulations of young/open SCs, which
are extensively discussed in Di Carlo et al. (2020a, b, 2021).The initial
distribution of stars in space is modeled with fractal initial conditions
(Küpper et al., 2011; Ballone et al., 2020, 2021; Torniamenti et al., 2021),
to mimic the asymmetry and clumpiness of star forming regions (Goodwin &
Whitworth, 2004). The initial stellar mass $M_{\rm SC}$ of every SC ranges
from $10^{3}$ M⊙ to $3\times{}10^{4}$ M⊙. Each cluster mass is drawn from a
$dN/dM_{\rm SC}\propto M_{\rm SC}^{-2}$ distribution, i.e. the initial mass
function of YSCs in the MW described in Lada & Lada (2003b). We calculate the
initial half-mass radius $r_{\rm h}$ using the relationship from Marks &
Kroupa (2012).

The SC initial conditions are generated using mcluster (Küpper et al., 2011).
We adopt the initial mass function from Kroupa (2001), with a minimum stellar
mass of 0.1 M⊙ and a maximum stellar mass of 150 M⊙. The initial total binary
fraction is $f_{\mathrm{bin}}=0.4$. Mass ratios are drawn from a distribution
$\mathcal{P}(q)\propto{}q^{-0.1}$ (where $q=m_{2}/m_{1}\in{}[0.1,1]$, Sana et
al. 2012). All the stars with $\,m\,\geq{}5$ M⊙ are in binary systems, while
the stars with $m\,<\,5$ M⊙ are paired stochastically until the initial binary
fraction is reached. This is consistent with the multiplicity properties of
O/B-type stars (Sana et al. 2012; Moe & Di Stefano 2017). Eccentricities and
orbital periods are drawn from the distributions from Sana et al. (2012).

We set the efficiency of common envelope ejection to $\alpha=5$. We adopt the
rapid core-collapse supernovae model, described in Fryer et al. (2012). Natal
kicks are randomly drawn from a Maxwellian velocity distribution. A one-
dimensional root mean square velocity $\sigma=15\,\mathrm{km/s}$ is adopted
for both core-collapse supernovae and electron-capture supernovae. We have
simulated SCs with three different metallicities: $Z=0.02,$ 0.002 and 0.0002.
Each SC is simulated for 100 Myr in a static solar neighbourhood-like tidal
field (Wang et al., 2016). From each simulation, we extract all the binaries
composed of a BH and a luminous companion which are still in the SC at the end
of the simulations, as well as all the binaries of the same type which escape
the SC.

###  2.2 Synthetic MW of Binaries and Clusters

Following Chawla et al. (2022), we seed both IBs and SCs following the star-
formation history of a cosmological zoom-in simulation of a MW-mass galaxy
(specifically the m12i galaxy from the Latte simulations, Wetzel et al., 2016)
from the Feedback In Realistic Environments (FIRE-2) suite of galaxies
(Hopkins, 2015; Hopkins et al., 2018). We extract from the publicly-available
data every star particle formed in the simulation (which have masses at
formation of $7070M_{\odot}$), including the particle’s metallicity and
formation time. These are then binned into low ($[\rm{Fe}/\rm{H}]<-1.5$),
intermediate ($-1.5<[\rm{Fe}/\rm{H}]<-0.5$) and high ($-0.5<[\rm{Fe}/\rm{H}]$)
metallicity bins, designed to match the metallicities of the original
nbody6++gpu cluster simulations. Each metallicity bin is then further divided
into 100 equal intervals of cosmic time between 0 and 13.78 Gyr, providing us
with the metallicity-dependent star formation rate of our MW-analogue galaxy.

To convert this star formation per unit time into stars and clusters, we begin
by assuming that all star formation occurs in clusters (Lada & Lada, 2003b).
Since the lowest-mass clusters from Di Carlo et al. (2021) begin at
$1000M_{\odot}$, we assume that all star formation clumps between
$100-1000M_{\odot}$ rapidly dissolve into isolated stars and binaries, while
clusters above $1000M_{\odot}$ remain bound and undergo dynamical
encounters.111For an extensive analysis on BH-star binaries in SCs with masses
smaller than 1000 M⊙, we refer to Rastello et al. (2023). Along with our
cluster initial mass function ($dN/dM_{\rm SC}\propto M_{\rm SC}^{-2}$), this
suggests that $10\%$ of star formation should occur in open clusters.

We seed our clusters in the m12i galaxy in the following way: for each bin of
star formation, we draw randomly with replacement from the cluster sample,
subtracting that cluster’s mass from the total mass of star particles formed
in that bin, until we have completely turned that 10% of the star formation
into clusters. Each time we draw a cluster, we randomly pick an m12i star
particle from the same time/metallicity bin, and assign that particle’s exact
formation time and present-day position in the m12i galaxy to the BH-star
binaries produced by that cluster. Note that this procedure generates more
clusters that present in the original catalog, meaning that many of our
dynamical sources are sampled repeatedly.

The remaining 90% of star formation is turned into massive binaries in a
similar fashion: for each time/metallicity bin, the remaining star formation
is converted directly into star particles, assuming the primary masses (that
will become BHs) follow a $1/M^{2.3}$ initial mass function (IMF) starting at
$18M_{\odot}$ (Kroupa, 2001) and up to $150\,$M⊙, while the secondaries are
drawn from a uniform mass ratio distribution (Mazeh et al., 1992; Goldberg &
Mazeh, 1994) from $0.1M_{\odot}$ up to the primary masses. The binary orbital
periods are drawn randomly from a flat-in-log distribution (Abt, 1983), while
the eccentricities are drawn from a thermal distribution. Note that these
initial conditions (other than the IMF) are different from those assumed in
§2.1: they are instead chosen to be consistent with Chawla et al. (2022), with
a specific focus on forming BH-star binaries. These systems are then evolved
to the present day as described in the next section. For the m12i galaxy, this
procedure yields approximately $9.4\times 10^{7}$ binaries with at least one
BH progenitor.

###  2.3 Evolving Binaries to the Present Day

For consistency with El-Badry et al. (2023a, b), we evolve both IB and SC
binaries to the present day using the binary population synthesis code COSMIC
(Breivik et al., 2020). The binaries are evolved using the same stellar and
binary evolution parameters and recipes used in Di Carlo et al. (2020b),
described in §2.1, to ensure that binary evolution is treated identically in
both populations. IBs are evolved starting from the Zero-Age Main Sequence
(ZAMS) with absolute metallicities of 0.02, 0.002 and 0.0002 (to match the
original cluster simulations), from the time when they were formed ($T_{\rm
form}$) in the galaxy until the current age of the Universe
($13.78\,\rm{Gyr}$).

The evolution of SC binaries is restarted from the last recorded state of the
binary in the $N$-body simulation—that is, either the time of binary ejection
or the end of the cluster simulation at 100 Myr. Restarting the evolution of
cluster binaries yields additional complications, as the binary evolution in
COSMIC for systems having undergone mass loss or accretion depend on the
current mass of the star, its effective initial mass, and the change in its
effective age due to the mass loss/transfer (referred to as $M(t)$, $M_{0}$,
and the stellar “epoch” in Hurley et al., 2000, respectively). $M_{0}$ and
epoch are not part of the default nbody6++gpu output, but since we are
primarily interested in BHs and MS stars, this did not present a major
difficulty: in COSMIC, $M_{0}$ is equal to $M_{t}$ for both BHs and MS stars,
while the epoch is only relevant for the evolution calculations of MS stars.
We assume $\rm{epoch}=0\,\rm{Myr}$ for MS stars when restarting their
evolution. As a test, we ran the entire population with epoch values of 0 and
100 Myr (the minimum and maximum values it could have given the SC integration
time) and found the difference in the population-level statistics to be
negligible (only 1% systems show a relative difference greater than $\sim$1%
in their final masses and orbital periods between $\rm{epoch}=0\,\rm{Myr}$ and
$\rm{epoch}=100\,\rm{Myr}$). While restarting post-MS stars in COSMIC is
substantially more complicated, none of the luminous companion stars evolved
beyond the MS before the end of the SC simulations (as opposed to evolving
onto the giant branch later in the galactic field). In our analysis, we
distinguish between BH-MS binaries and BH-Giant binaries. We follow the
criterion from Drout et al. (2012) and classify post-MS stars with $T_{\rm
eff}\leq 4800\,K$ as giants. This includes stars with BSE stellar types
associated with red giants, core helium burning stars and asymptotic giant
branch stars (types 3, 4 and 5).

For our analysis, we only select systems with an age $t_{\mathrm{age}}\geq
100$ Myr. This choice is justified for three reasons: first, our SCs are
evolved for 100 Myr, and therefore we are unable to accurately track systems
within SCs which formed in the past 100 Myr. Second, the star-formation
history of the FIRE-2 simulations can be larger than that observed in MW-mass
galaxies at late times by a factor of a few (§4.5 of Hafen et al., 2022). This
is true in particular for the m12i galaxy we consider, whose $z=0$ star-
formation rate ($6M_{\odot}/\rm{yr}$, Wetzel et al., 2016) is $\gtrsim 3-4$
times the best observational estimate of star formation in the MW ($1.65\pm
0.19M_{\odot}/\rm{yr}$, Licquia & Newman, 2015). Finally, one of the main
focuses of this paper is to understand the formation channels of Gaia BH1 and
Gaia BH2, which are both older than $\sim 1$ Gyr (and likely significantly
older, especially in the case of Gaia BH2, El-Badry et al., 2023b, §3.8). We
briefly discuss the implications of this cutoff (and what our results would
look like including such systems) in §5.3.

###  2.4 Gaia Detectability

Following Chawla et al. (2022), we further determine whether Gaia can
astrometrically resolve each orbit based on the motion of the luminous
companion based on an optimistic and pessimistic criteria. In the optimistic
case, we assume that Gaia can resolve any orbit for which the star’s motion is
at least as large as Gaia’s single pointing precision, $\sigma_{G}$, when it
is projected on the 2d plane of the sky. In the pessimistic case, we require
the projected star’s orbit to be at least three times as large as
$\sigma_{G}$.

We also consider the cut in the parallax signal-to-noise ratio (SNR),
$\omega/\sigma_{\omega}$, which was applied in in Gaia DR3. This cut was
applied as a function of orbital period as shown in Figure 1. The results are
shown for the SC BH-star binaries in red and IB binaries in blue. Dynamically-
formed BH-star binaries are preferentially located above the SNR cut while IB
BH-star binaries which lie above the cut are mostly due to random-chance
placements that are unrealistically close to Earth.

To increase our statistical sample of “detected” systems, we consider
different positions for Earth in the m12i simulation and combine the results.
We use 100 different equally-spaced starting positions along a ring in the
disk with radius $8.5\,\rm{Kpc}$. For what follows, the results for our
synthetic Gaia population are divided by 100 when quoting the number of
detectable binaries and all other relevant quantities.

Figure 1: The orbital period vs parallax SNR, $\omega/\sigma_{\omega}$ are
shown for BH-star binaries formed in clusters (red) and in isolated binaries
(blue) from all of our $100$ MW-like simulation realizations. Dots show BH-MS
binaries and crosses show BH-Giant binaries. The dashed black like shows the
parallax SNR cut placed as a function of orbital period for Gaia DR3 such that
only binaries above the line are included in the data set.

##  3 Results

###  3.1 Intrinisic Population

We first explore the intrinsic population of BH-star binary systems in the
galaxy; that is, the full population before considering our Gaia detectability
criterion. In all our BH-star systems, the luminous member is either a MS or
giant star. The exact numbers of systems per formation channel and stellar
type is reported in Table 1. Based on the m12i simulation, we predict that the
MW hosts a total of $\sim 2\times 10^{5}$ BH-star systems from the IB and the
young/open SC channels, with $\sim 81$% of the systems formed dynamically and
$\sim 19$% formed in isolation. This means that according to our models, SCs
produce nearly 4 times more systems than IBs. However, we must emphasize that
we have assumed that only 10% of the star formation occurs in young SCs. If we
take this into account, it is obvious that SCs are _dramatically_ more
efficient than IBs at producing BH-star systems, with SCs producing $\sim 40$
times as many systems overall, per unit mass, as IB evolution. This is driven
largely by the production of BH-MS star systems; for evolved systems however,
the situation is reversed, with IBs producing 2.4 times more BH-Giant binaries
than SCs.

Figure 2 shows the distribution of the orbital parameters of the systems in
the intrinsic population. SCs are much more efficient at producing systems
with large orbital periods for both BH-MS and BH-Giant binaries. IB systems
have eccentricities up to $\sim 0.84$, while the most eccentric dynamical
system has an eccentricity of $\sim 0.99$. IBs produce fewer systems with more
massive stars: the maximum value of $M_{*}$ for IB is $\sim 5$ M⊙, while it is
$\sim 7$ M⊙ for SCs (though this does depend on the 100 Myr cutoff we have
employed, as described in §5.3). The maximum BH mass is $\sim 50$ M⊙ for IBs
and $\sim 67$ M⊙ for the SC channel. The SC channel is thus able to form
systems with eccentricities and masses not accessible from IB evolution.

Channel | All | BH-MS | BH-Giant
---|---|---|---
Formation Efficiency $\mathrm{[M_{\odot}^{-1}]}$
Isolated | $4.55\times 10^{-7}$ | $4.21\times 10^{-7}$ | $3.40\times 10^{-8}$
Dynamical | $1.81\times 10^{-5}$ | $1.80\times 10^{-5}$ | $1.28\times 10^{-7}$
Total | $2.22\times 10^{-6}$ | $2.18\times 10^{-6}$ | $4.34\times 10^{-8}$
Intrinsic Population
Isolated | 35229 | 32596 | 2633
Dynamical | 155733 | 154633 | 1100
Total | 190962 | 187229 | 3733
Gaia DR3 Detections
Isolated | 0.7 | 0.4 | 0.3
Dynamical | 6.6 | 5.2 | 1.4
Total | 7.3 | 5.6 | 1.7
Gaia Detections (including future data releases)
Isolated | 8-32 | 5-19 | 3-13
Dynamical | 53-178 | 40-129 | 13-49
Total | 61-210 | 45-148 | 16-62
Table 1: Formation efficiency (i.e. number of systems in the intrinsic
population produced per unit simulated stellar mass), number of total BH-star
binary systems (intrinsic population), number of expected Gaia DR3 detections
and number of pessimistic-optimistic expected Gaia detections including future
data releases. Column 1: formation channel of the binary; column 2: value for
all BH-star binaries; column 3: value for BH-main sequence binaries only;
column 4: value for BH-Giant binaries only. Decimal values come from our
average over 100 different observations (see §2.4). For Gaia DR3 detections,
our optimistic and pessimistic values are identical. Figure 2:  Orbital
parameters of systems in the intrinsic population. From top to bottom we show
the distributions of binary orbital periods, eccentricities, stellar masses
$M_{*}$, and BH masses $M_{\mathrm{BH}}$, with SC binaries shown in red lines
and IB shown in blue lines. Solid lines represent BH-MS binaries, while dashed
lines show BH-Giant binaries. The yellow and purple vertical lines represent
the values of Gaia BH1 and Gaia BH2 respectively.

###  3.2 Population Detectable by Gaia

The Gaia mission provides a unique opportunity to detect a subset of the
intrinsic population of BH-star binary systems in the Milky Way. Table 1
presents the number of systems in our data-set according to their
detectability by Gaia. The table displays two subsets: the first containing
the number of systems detectable by Gaia DR3 in our model, while the second
includes systems detectable by both Gaia DR3 and future Gaia data releases. It
is evident how SCs are significantly more efficient at producing detectable
binaries. In particular, we see that all the 7 expected Gaia DR3 detections
are expected to come from the dynamical formation channel. If we also consider
detections with future Gaia data releases, we predict that Gaia will detect
between 8 and 32 binaries from IBs, and between 53 and 178 binaries from the
dynamical channel, i.e. a factor of $\sim 6$ more dynamical systems. In
general, our models produce more detectable BH-MS binaries than BH-Giant
binaries.

Figure 3 shows the distributions and correlations of the orbital parameters
for the detectable population. It is evident that only SCs are able to produce
detectable systems with large eccentricities and with large BH masses. In
particular, the detection of a system with $e\gtrsim 0.5$ or with
$M_{\mathrm{BH}}\gtrsim 10$ M⊙ would be a smoking gun of dynamical formation,
as IBs are unable to form detectable systems with these characteristics in our
models. The SC channel is more efficient at producing detectable binaries with
large orbital periods as well. However, we reiterate that these results are
for systems older than 100 Myr; we explore the implications of this cutoff in
§5.3.

Figure 3:  Correlations and distributions of orbital period, eccentricity,
mass of the star and mass of the BH of detectable systems. Dynamical systems
are shown in red, while isolated systems are shown in blue. Gaia BH1 is shown
in yellow, while Gaia BH2 is shown in violet. BH-MS systems are represented by
filled circles in the scatter panels and by solid lines in the histogram
panels. BH-Giant systems are represented by crosses in the scatter panels and
by dashed lines in the histogram panels. Scatter plots show systems from all
the 100 observers around the Galaxy (see §2.4), while histograms show the
number of detections $N_{\mathrm{det}}=N/100$.

###  3.3 Dynamical Formation and Star Cluster Properties

The results from our analysis of the intrinsic and detectable populations
strongly suggest that the majority of these systems are formed dynamically in
SCs. Given that, it is informative to understand how many detectable systems
produced by SC evolution come from primordial binaries (i.e. binaries present
in the SC since the beginning of the simulation that may have been altered
through dynamical processes) or from binaries that were assembled later
dynamically. We find that approximately $94\%$ of detections stem from these
dynamically-assembled binaries, while the remaining $6\%$ arise from
primordial binaries; $\sim 20$% of the dynamically assembled binaries
initially assemble as a star-star system (that later evolves to a BH-star
system), while the remaining $\sim 80$% form dynamically from a BH and a star.

We can also ask which SCs produce the majority of the BH-star binaries. In
Figure 4, we show the number of BH-star binaries produced by SCs of a given
mass. Even though the cluster catalog is sampled following a $\propto M_{\rm
SC}^{-2}$ cluster mass function (containing many more low-mass clusters), the
intrinsic population exhibits a relatively uniform distribution across SC
masses, while the detectable population predominantly arises from SCs with
$M_{\mathrm{SC}}\gtrsim 3\times 10^{3}$ M⊙. For insights into the formation of
BH-MS binaries in SCs with lower masses ($M_{\mathrm{SC}}=300-1000$ M⊙),
Rastello et al. (2023) show a comprehensive analysis that compares results
from low- and high-mass SCs.

The overwhelming majority of observable binary systems originate from SCs with
high metallicity (Z=0.02), with $\sim 99.84\%$ of them forming in such
clusters. This is because of the metallicity-dependent star formation history
of the MW. Massive metal-rich SCs are thus very efficient at producing
detectable binaries, and we expect the majority of BH-star systems in the MW
and current/future Gaia detections to have formed in such environments. In a
recent study, Tanikawa et al. (2023) estimated that the MW hosts around
$1.6\times 10^{4}$ BH-star systems by simulating SCs with
$M_{\mathrm{SC}}=1000$ M⊙. In contrast, our analysis suggests a much larger
number of $\sim 1.9\times 10^{5}$ BH-star systems, highlighting the crucial
role of more massive SCs to comprehend the formation and evolution of these
systems.

We also highlight that $\sim 85\%$ of the BH-star binaries formed in SCs are
retained by their host SC at the end of the simulations. Despite the
relatively low initial mass of the SCs considered in this study, most of them
have not undergone complete disruption by the end of the simulations and would
likely survive longer (see e.g. Torniamenti et al., 2022). While it is
possible that some of the retained binaries may eventually disrupt due to
dynamical encounters, it is also highly likely that new BH-star binaries will
form as a result of ongoing dynamical interactions.

Figure 4:  Distributions of the initial SC masses $M_{\mathrm{SC}}$ which
produce BH-star binaries. The intrinsic population is shown by the solid red
line, while the detectable population is shown by the dotted red line. The
initial masses of the simulated SCs are drawn from a $dN/dM_{\rm SC}\propto
M_{\rm SC}^{-2}$ distribution (Lada & Lada, 2003b), meaning that we simulated
more low mass SCs and less high mass SCs.

##  4 Forming Gaia BH1 and Gaia BH2

In order to understand how the BH-star binary systems detected by Gaia might
have formed, we compare their properties to those of our simulated
populations. For both the IB and the SC channel, we find the most similar
system to Gaia BH1 and Gaia BH2 in the detectable population, and describe
their formation history. The identification of the most similar system entails
assessing the fractional differences in key parameters, namely the orbital
period, primary mass, and secondary mass, between the simulated binaries and
the Gaia BHs. The fractional difference is defined as
$|x-x_{\mathrm{GAIA}}|/x_{\mathrm{GAIA}}$, where $x$ is the parameter of the
object in the simulated population and $x_{\mathrm{GAIA}}$ is the
corresponding parameter of the Gaia BH. The system that exhibits the smallest
maximum fractional difference across these parameters (i.e. the one that
deviates the least across all the parameters considered) is regarded as the
most similar to the Gaia BH. Due to the relatively limited statistical data
available for the SC channel, we make the assumption that eccentricities are
subject to randomization through dynamical encounters in SCs, where the
eccentricity distribution is expected to follow a thermal profile. Thus, we
exclude the eccentricity in the evaluation of the most similar systems. The
formation histories of these systems is summarized in Figure 5.

Figure 5: Formation history of the most similar systems to Gaia BH1 (left
panels) and Gaia BH2 (right panels) in our simulations. Top panels show
systems from the isolated binary evolution channel; bottom panels show systems
from the star clusters channel. Main sequence stars (with label MS) are
represented as yellow stars; core helium burning stars (label cHeB) are
visualized as red stars; naked helium stars (label Naked He) are represented
as blue stars; black holes (label BH) are shown as black circles. The mass of
each object is shown next to them. The orbital period in days of the final
binaries is shown underneath. The time axis and the size of the objects and
orbits are not to scale. Even though Gaia BH-like systems can be formed by
both the formation channels, our results on rates and ages presented in §5.2
strongly suggest that the dynamical scenario is the primary formation pathway
for the Gaia BHs.

###  4.1 Gaia BH1

The parameters of the most similar IB system to Gaia BH1 fall within
approximately $25\%$ of those of Gaia BH1. The system has an orbital period
$P\sim 232$ days, $M_{\mathrm{BH}}\simeq 8.4$M⊙ and $M_{*}\simeq 1.05$M⊙. The
initial conditions of the system are $P\simeq 3486$ days,
$M_{\mathrm{1}}\simeq 55.1$M⊙, $M_{2}\simeq 1.04$M⊙ and a metallicity of
$Z=0.02$. The system undergoes a common envelope episode that significantly
shrinks its orbit, and then turns into a BH-MS binary after 4.8 Myr. The
binary then evolves unperturbed to the present day. The system has an age of
$\sim 1.5$ Gyr.

The closest system from the SC channel has the parameters within $\sim 32\%$
of Gaia BH1. The system has an orbital period $P\sim 130$ days,
$M_{\mathrm{BH}}\simeq 6.6$M⊙ and $M_{*}\simeq 0.6$M⊙. The system assembles
dynamically inside a SC with metallicity $Z=0.02$ and initial mass
$M_{\mathrm{SC}}=3204$ M⊙. The binary forms $\sim 5$ Myr after the beginning
of the simulation. The primary turns into a giant, triggering a common
envelope episode that shrinks the orbit. Afterwards, the envelope is ejected,
and the primary turns into a BH. $\sim 10$ Myr after its formation, the system
is dynamically ejected from the SC. The age of this system is $\sim 3$ Gyr.

Both the IB and SC models struggle to reproduce a system that very closely
resembles Gaia BH1. Rastello et al. (2023) suggests that the formation of Gaia
BH1 might be explained by less massive SCs, with masses ranging from 300M⊙ to
1000M⊙. These lower mass star clusters could provide a favorable environment
for the formation of a system with properties similar to Gaia BH1. Based on
our analysis of the age distribution and of our detection rates for BH-MS
systems in Gaia DR3, we draw the conclusion that the dynamical formation
channel is the most favorable for Gaia BH1, as we better explain in §5.2.

###  4.2 Gaia BH2

The closest IB system to Gaia BH2 is a BH-Giant binary with $P\sim 1044$ days,
$M_{\mathrm{BH}}\simeq 8.0$M⊙ and $M_{*}\simeq 1.3$M⊙. The values lie within
$\sim 19\%$ of the parameters of Gaia BH2. The initial conditions of the
system are $P\simeq 3967$ days, $M_{\mathrm{1}}\simeq 45.8$M⊙, $M_{2}\simeq
1.36$M⊙. The system has a metallicity of $Z=0.02$, and undergoes a common
envelope episode throughout its evolution. The system has an age of $\sim 4.1$
Gyr, which is not compatible with the age of Gaia BH2.

The best candidate from the SC channel has the parameters within $13\%$ of the
Gaia BH2 parameters. The system has an orbital period $P\sim 1197$ days,
$M_{\mathrm{BH}}\simeq 8.1$M⊙ and $M_{*}\simeq 1.2$M⊙. The binary is
dynamically assembled in a SC with metallicity $Z=0.02$, initial mass
$M_{\mathrm{SC}}=14284$ M⊙ and an age of $\sim 6$ Gyr. The system forms as a
BH-MS binary $\sim 39$ Myr after the beginning of the simulation, and it later
evolves to a BH-Giant. The system is retained by its host SC at the end of the
simulation.

Both channels produce systems that exhibit close parameter matches to Gaia
BH2. Among these, the dynamical scenario yields the closest resemblance. If we
also take into account that Gaia BH2 is estimated to be older than 5 Gyr, and
that our Gaia DR3 detections for BH-Giant systems come entirely from SCs, we
conclude that Gaia BH2 has likely formed via the dynamical formation channel.
We provide a more extensive analysis on this in §5.2.

##  5 Discussion

###  5.1 Dynamical versus Isolated

We see from Table 1 that the dynamical channel produces $\sim 4.3$ times more
BH-star systems than the isolated channel. Since we consider that only 10% of
the total star formation occurs into SCs, we find that the number of systems
produced per unit stellar mass by SCs is $\sim 40$ times larger than the
number of systems produced by isolated binary evolution. This discrepancy can
be attributed to several mechanisms. Stellar dynamics in SCs significantly
affects the formation and evolution of binaries. Multi-body encounters can
dynamically assemble BH-star or star-star binaries which can later evolve into
BH-star binaries, resulting in combinations of orbital parameters that are
inaccessible in IB evolution. Furthermore, three-body encounters between a
binary and a single object can modify the eccentricity and orbital period of
binaries, triggering or preventing mass transfer episodes (e.g., Roche-lobe
overflow and common envelope evolution) that may or may not occur if the
binary evolved in isolation. Common envelope episodes significantly reduce the
orbital period of binaries. From Figure 2 and Figure 3, we see how the
dynamical channel is extremely more efficient at producing binaries with
larger orbital periods. Indeed, we find that 99.4% of all the simulated IB
systems which become BH-star binaries undergo at least one common envelope
episode throughout their evolution. This can be easily avoided in SCs; as we
said in §3.3, the majority of detectable SC binaries dynamically assembles as
a BH-star binary already, which means that they do not undergo a common
envelope episode222The only way for such binaries to undergo a common envelope
episode and still be classified as BH-star binaries is if the following
sequence of events occurs: i) the star becomes a giant; ii) its expansion
triggers a common envelope; iii) the envelope is ejected by the BH companion,
leaving behind a BH and a naked He star. We do not find any BH-naked He star
binary in our populations.. The dynamical channel is thus not only more
efficient, but can also leave unique fingerprints on the properties of the
binaries.

###  5.2 Gaia BHs formation channels

Our models provide several hints on how the dynamical scenario is a more
favorable formation pathway for Gaia BH1 and Gaia BH2. According to our
results in Table 1, all the expected Gaia DR3 detections come from SCs,
suggesting that Gaia BHs have more likely formed dynamically. In order to form
systems like Gaia BH1 or Gaia BH2 through IB evolution, El-Badry et al.
(2023b) and El-Badry et al. (2023a) highlight the necessity of assuming an
exceedingly large and potentially unrealistic value for the common envelope
efficiency ($\alpha\approx 12.8$ for Gaia BH1). Lower values of $\alpha$ would
have led to the merger of the binary components, preventing the formation of
Gaia BHs through this channel. As discussed in §5.1, dynamics offers an
alternative pathway for the formation of BH-star systems that completely
bypasses the common envelope phase.

The age distribution of the detected systems is shown in Figure 6. The
efficiency of both channels is comparable at 0.1 Gyr, but the dynamical
channel is much more efficient at forming detectable systems with older ages.
Gaia BH1 has an age $\gtrsim$1 Gyr, while Gaia BH2 has an age $\gtrsim$5 Gyr,
indicating that both Gaia BHs, particularly Gaia BH2, have likely formed
dynamically, according to our models.

Rastello et al. (2023) show that both lower mass (300-1000M⊙) and higher mass
(1000-30000 M⊙) young SCs are efficient in forming Gaia BH1-like systems. Gaia
BHs may also have formed via other channels not taken into account in this
study, such as formation in isolated hierarchical triples or hierarchical
triples formed in SCs (see e.g. Trani et al., 2022). Another plausible channel
could be the dynamical formation in a globular cluster (GC); while the
Galactic orbit of both Gaia BHs is not aligned with any GC in the MW (El-Badry
et al., 2023b, a), these systems might have dynamically formed in GCs which
have already completely disrupted.

###  5.3 BH-star ages and the 100 Myr age cutoff

Figure 6:  Distributions of the ages of systems in our datasets. The intrinsic population is shown by solid lines, while Gaia detectable systems are shown by dotted lines. Dynamical systems are shown in red, while isolated systems are shown in blue. The two dashed vertical lines show the lower limits of the age of Gaia BH1 and Gaia BH2. The grey shaded area on the left represents the systems not taken into account due to the age cutoff described in §2.3. Channel | All | BH-MS | BH-Giant
---|---|---|---
Formation Efficiency $\mathrm{[M_{\odot}^{-1}]}$
Isolated | $5.13\times 10^{-7}$ | $4.79\times 10^{-7}$ | $3.47\times 10^{-8}$
Dynamical | $1.84\times 10^{-5}$ | $1.83\times 10^{-5}$ | $1.30\times 10^{-7}$
Total | $2.30\times 10^{-6}$ | $2.26\times 10^{-6}$ | $4.42\times 10^{-8}$
Intrinsic Population
Isolated | 39737 | 37048 | 2689
Dynamical | 158329 | 157214 | 1115
Total | 198066 | 194262 | 3804
Gaia DR3 Detections
Isolated | 87 | 82 | 5
Dynamical | 14 | 12 | 2
Total | 101 | 94 | 7
Gaia Detections (including future data releases)
Isolated | 286-455 | 272-430 | 14-25
Dynamical | 88-272 | 73-220 | 15-52
Total | 374-727 | 345-650 | 29-77
Table 2: Same as Table 1, but without the 100 Myr formation time cutoff.

Throughout this paper, we have argued that SCs dominate the production of old
BH-star systems in the MW. We show this explicitly in Figure 6, which displays
the distribution of ages of the systems, and reveals that the dynamical
formation channel appears to be more efficient at producing older systems,
particularly those older than 1 Gyr.

Of course, it is also obvious in Figure 6 that our 100 Myr cutoff has removed
a large number of young systems, particularly from the IB channel, that
contribute significantly to the detectable population. For comparison, we
present the results obtained without applying the 100 Myr age cutoff (which
was justified in §2.3), starting with Table 2 which reports the number of
systems in the intrinsic and detectable populations. Although there is an
increase in the number of systems in the isolated intrinsic population, the
dynamical channel remains dominant, constituting approximately $80\%$ of the
total intrinsic population.

However, we observe a significant increase in the number of isolated
detectable systems. The expected number of IB detections by Gaia DR3 increases
from 0 to 87, while the total number of IB Gaia detections increases from 8-32
to 286-455. Conversely, the number of detectable systems from the dynamical
channel barely increases. This disparity can be attributed to the fact that
systems in the isolated population are predominantly very young (see Figure 6)
and preferentially host bright massive stars that are significantly easier to
detect. This trend is evident from Figure 7 and Figure 8, which illustrate the
significant impact of the cutoff on the $M_{*}$ distribution. The difference
comes almost entirely from a huge population of young massive O-stars. The IB
detectable systems younger than 20 Myr have an average BH mass of 8 M⊙ and an
average star mass of 41 M⊙, and constitute $\sim 86\%$ of the total number of
detectable systems in the whole IB catalogue without the 100 Myr age cutoff.
We emphasize that our predictions for the formation channels of both Gaia BH1
and Gaia BH2 should remain unaffected by this cutoff, as they both are older
than 1 Gyr.

Figure 7:  Distributions of orbital parameters of systems in the intrinsic
population, including systems with age $<100$ Myr. From top to bottom: orbital
period, orbital eccentricity $e$, star mass $M_{*}$, BH mass
$M_{\mathrm{BH}}$. SC binaries are shown by red lines, IB are shown by blue
lines. Solid lines represent BH-MS binaries, while dashed lines show BH-Giant
binaries. The yellow and purple vertical lines represent the values of Gaia
BH1 and Gaia BH2 respectively. Figure 8:  Same as Figure 7, but for systems
detectable in future Gaia data releases.

##  6 Conclusions

The MW is very likely populated by a large number of stellar BHs, many of
which probably reside in binary systems. The recent discovery of two BH-star
binaries in Gaia DR3 has sparkled new interest in understanding the formation
channels and population characteristics of such systems. In this paper, we
have investigated the formation of BH-star binaries in young SCs and via IB
evolution.

According to our simulations, the MW harbors a total of $\sim 2\times 10^{5}$
BH-star systems of which $\sim 81$% formed dynamically and $\sim 19$% formed
in isolation. We find that dynamical formation in SCs is nearly 40 times more
efficient per unit stellar mass at producing BH-star binaries compared to
isolated binary evolution. Dynamical systems tend to have larger orbital
periods and eccentricities than the isolated ones. We expect that a total of 7
BH-star systems are present in the Gaia DR3 data, all of which come from the
dynamical channel. We also predict between 61 and 210 detections from the
whole Gaia mission, $\sim 83$% of which come from SCs. Overall, dynamics
enhances dramatically the number of BH-star binaries, both in the intrinsic
and the detectable populations.

We compare our detectable populations with Gaia BHs, and we conclude that our
models support the dynamical scenario as the primary formation pathway for
Gaia BH1 and Gaia BH2. Our results suggest that identifying systems with an an
eccentricity greater than $\sim 0.5$, or with black hole mass
$M_{\mathrm{BH}}$ exceeding 10 M⊙ would provide compelling evidence for their
dynamical formation. With new detections from future Gaia data releases
expected over the next few years, these results will help disentangle the
formation channels of BH-star systems and the characteristics of the BH
population in the MW.

Finally, our analysis also revealed the presence of some BH-Giant systems with
3 M⊙ $<M_{\mathrm{BH}}<$ 5 M⊙, i.e. below the minimum BH mass allowed by the
rapid supernova mechanism (Fryer et al., 2012). These peculiar systems form
through accretion of matter by neutron stars from their binary companions,
leading to the formation of BHs through accretion-induced collapse (AIC). None
of these systems are detectable by Gaia according to our analysis, but we have
not yet compared these systems (identified only in the IB channel) to their
dynamical counterparts. Efforts to better characterize these systems are
currently underway.

We thank Michela Mapelli, Sara Rastello, Giuliano Iorio and all the members of
the DEMOBLACK group for the useful discussions. This work was supported by NSF
Grant AST-2009916 to The University of North Carolina at Chapel Hill and
Carnegie Mellon University. The $N$-body simulations of young star clusters
adopted here are an open-data product of the European Research Council (ERC)
Consolidator grant DEMOBLACK, contract no. 770017 (PI: Michela Mapelli). CR
acknowledges support from a Charles E. Kaufman Foundation New Investigator
Research Grant, an Alfred P. Sloan Research Fellowship, and a David and Lucile
Packard Foundation Fellowship. The Flatiron Institute is supported by the
Simons Foundation.

## References

  * Abt (1983) Abt, H. A. 1983, ARA&A, 21, 343, doi: 10.1146/annurev.aa.21.090183.002015
  * Ballone et al. (2020) Ballone, A., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 496, 49, doi: 10.1093/mnras/staa1383
  * Ballone et al. (2021) Ballone, A., Torniamenti, S., Mapelli, M., et al. 2021, MNRAS, 501, 2920, doi: 10.1093/mnras/staa3763
  * Breivik et al. (2017) Breivik, K., Chatterjee, S., & Larson, S. L. 2017, ApJ, 850, L13, doi: 10.3847/2041-8213/aa97d5
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Chakrabarti et al. (2022) Chakrabarti, S., Williams, T. G., Druard, C., & et al. 2022, Monthly Notices of the Royal Astronomical Society, 520, 1815, doi: 10.1093/mnras/stab579
  * Chawla et al. (2022) Chawla, C., Chatterjee, S., Breivik, K., et al. 2022, ApJ, 931, 107, doi: 10.3847/1538-4357/ac60a5
  * Corral & Santana (2016) Corral, A., & Santana, R. 2016, International Journal of Modern Physics D, 25, 1644001, doi: 10.1142/S0218271816440017
  * Di Carlo et al. (2019) Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947, doi: 10.1093/mnras/stz1453
  * Di Carlo et al. (2020a) Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020a, MNRAS, 497, 1043, doi: 10.1093/mnras/staa1997
  * Di Carlo et al. (2020b) Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020b, MNRAS, 498, 495, doi: 10.1093/mnras/staa2286
  * Di Carlo et al. (2021) Di Carlo, U. N., Mapelli, M., Pasquato, M., et al. 2021, MNRAS, 507, 5132, doi: 10.1093/mnras/stab2390
  * Drout et al. (2012) Drout, M. R., Massey, P., & Meynet, G. 2012, ApJ, 750, 97, doi: 10.1088/0004-637X/750/2/97
  * El-Badry et al. (2023a) El-Badry, K., Rix, H.-W., Quataert, E., et al. 2023a, MNRAS, 518, 1057, doi: 10.1093/mnras/stac3140
  * El-Badry et al. (2023b) El-Badry, K., Rix, H.-W., Cendes, Y., et al. 2023b, MNRAS, 521, 4323, doi: 10.1093/mnras/stad799
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91, doi: 10.1088/0004-637X/749/1/91
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Garmany et al. (1982) Garmany, C., Conti, P., & Massey, P. 1982, The Astrophysical Journal, 263, L33, doi: 10.1086/183356
  * Giacobbo & Mapelli (2018) Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011, doi: 10.1093/mnras/sty1999
  * Giacobbo & Mapelli (2019) —. 2019, MNRAS, 482, 2234, doi: 10.1093/mnras/sty2848
  * Giacobbo & Mapelli (2020) —. 2020, ApJ, 891, 141, doi: 10.3847/1538-4357/ab7335
  * Giacobbo et al. (2018) Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959, doi: 10.1093/mnras/stx2933
  * Giesers et al. (2019) Giesers, B., Bastian, N., Evans, C. J., & et al. 2019, Monthly Notices of the Royal Astronomical Society, 482, 2407, doi: 10.1093/mnras/sty2798
  * Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., & et al. 2018, Monthly Notices of the Royal Astronomical Society, 475, L15, doi: 10.1093/mnrasl/sly016
  * Goldberg & Mazeh (1994) Goldberg, D., & Mazeh, T. 1994, A&A, 282, 801
  * Goodwin & Whitworth (2004) Goodwin, S. P., & Whitworth, A. P. 2004, A&A, 413, 929, doi: 10.1051/0004-6361:20031529
  * Hafen et al. (2022) Hafen, Z., Stern, J., Bullock, J., et al. 2022, MNRAS, 514, 5056, doi: 10.1093/mnras/stac1603
  * Hopkins (2015) Hopkins, P. F. 2015, MNRAS, 450, 53, doi: 10.1093/mnras/stv195
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Kroupa (2001) Kroupa, P. 2001, Monthly Notices of the Royal Astronomical Society, 322, 231
  * Küpper et al. (2011) Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300, doi: 10.1111/j.1365-2966.2011.19412.x
  * Lada & Lada (2003a) Lada, C. J., & Lada, E. A. 2003a, Annual Review of Astronomy and Astrophysics, 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Lada & Lada (2003b) —. 2003b, Annual Review of Astronomy and Astrophysics, 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Licquia & Newman (2015) Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96, doi: 10.1088/0004-637X/806/1/96
  * Mapelli et al. (2017) Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422, doi: 10.1093/mnras/stx2123
  * Mapelli et al. (2020) Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76, doi: 10.3847/1538-4357/ab584d
  * Marks & Kroupa (2012) Marks, M., & Kroupa, P. 2012, A&A, 543, A8, doi: 10.1051/0004-6361/201118231
  * Mazeh et al. (1992) Mazeh, T., Goldberg, D., Duquennoy, A., & Mayor, M. 1992, ApJ, 401, 265, doi: 10.1086/172058
  * McClintock et al. (2006) McClintock, J. E., Shafee, R., Narayan, R., et al. 2006, The Astrophysical Journal, 652, 518, doi: 10.1086/507561
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224, doi: 10.1103/PhysRev.136.B1224
  * Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L., & Gieles, M. 2010, Annual Review of Astronomy and Astrophysics, 48, 431, doi: 10.1146/annurev-astro-081309-130834
  * Portegies Zwart et al. (1997) Portegies Zwart, S. F., Verbunt, F., & Ergma, E. 1997, A&A, 321, 207, doi: 10.48550/arXiv.astro-ph/9701037
  * Rastello et al. (2023) Rastello, S., Iorio, G., Mapelli, M., et al. 2023, arXiv
  * Reed (2003) Reed, B. C. 2003, AJ, 125, 2531, doi: 10.1086/374771
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, Annual Review of Astronomy and Astrophysics, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Tanikawa et al. (2023) Tanikawa, A., Cary, S., Shikauchi, M., Wang, L., & Fujii, M. S. 2023, arXiv e-prints, arXiv:2303.05743, doi: 10.48550/arXiv.2303.05743
  * Torniamenti et al. (2021) Torniamenti, S., Ballone, A., Mapelli, M., et al. 2021, MNRAS, 507, 2253, doi: 10.1093/mnras/stab2238
  * Torniamenti et al. (2022) Torniamenti, S., Rastello, S., Mapelli, M., et al. 2022, MNRAS, 517, 2953, doi: 10.1093/mnras/stac2841
  * Trani et al. (2022) Trani, A. A., Rastello, S., Di Carlo, U. N., et al. 2022, MNRAS, 511, 1362, doi: 10.1093/mnras/stac122
  * Wang et al. (2015) Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070, doi: 10.1093/mnras/stv817
  * Wang et al. (2016) —. 2016, MNRAS, 458, 1450, doi: 10.1093/mnras/stw274
  * Wetzel et al. (2016) Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, The Astrophysical Journal, 827, L23, doi: 10.3847/2041-8205/827/2/L23