#  Microscopic solutions for vortex clustering in two-band type-1.5
superconductor

Igor Timoshuk  Department of Physics, The Royal Institute of Technology,
Stockholm SE-10691, Sweden  Wallenberg Initiative Materials Science for
Sustainability, Department of Physics, The Royal Institute of Technology,
Stockholm SE-10691, Sweden  Egor Babaev  Department of Physics, The Royal
Institute of Technology, Stockholm SE-10691, Sweden  Wallenberg Initiative
Materials Science for Sustainability, Department of Physics, The Royal
Institute of Technology, Stockholm SE-10691, Sweden

###### Abstract

It has been proposed that two-band superconductors exhibit a distinct phase
characterized by two correlation lengths, where one is smaller and the other
larger than the magnetic field penetration length. This regime was coined
type-1.5 superconductivity, with a number of unconventional properties, such
as vortex clustering. However a fully microscopic solution for vortex clusters
has remained challenging due to computational complexities beyond
quasiclassical models. In this work, we present numerical solutions obtained
in a fully self-consistent two-band Bogoliubov-de-Gennes model. We show the
presence of discrepant correlation lengths leading to vortex clustering in
two-band superconductors.

The original work by Ginzburg and Landau introduced the concept of coherence
length $\xi^{GL}$ [1] and classified superconductors by a single Ginzburg-
Landau parameter: the ratio $\kappa^{GL}=\lambda/\xi^{GL}$ of two fundamental
length scales: the magnetic field penetration length $\lambda$ and the
coherence length $\xi^{GL}$. The latter is a fundamental length scale that
governs the asymptotic behavior of the modulus of superconducting gap
$|\Delta|e^{i\theta}$ or equivalently, up to a different prefactor [2],
Ginzburg-Landau order parameter field $|\Psi|e^{i\theta}$. The existence of
this order parameter and, hence, the fundamental length, $\xi^{GL}$ is
guaranteed by the fact that a superconductor breaks local $U(1)$ symmetry [1].
Within Ginzburg-Landau’s (GL) theory, the superconductor allows repulsively
interacting vortices for $\kappa^{GL}>1/\sqrt{2}$, called type-2
superconductivity [3]. The vortices form lattices when the magnetic field is
larger than the first critical magnetic field $H_{c1}$ and smaller than the
second critical magnetic field $H_{c2}$ [3]. These critical magnetic fields
were introduced in the series of experimental works by Shubnikov et al. [4,
5]; hence the state forming $H_{c1}<H<H_{c2}$ also referred to as Shubnikov’s
phase. In what follows, we will absorb the factor $\sqrt{2}$ into the
definition of coherence length $\xi\equiv\sqrt{2}\xi^{GL}$. Hence, the GL
criterion for type-2 superconductivity in these notations is
$\kappa\equiv\sqrt{2}\kappa^{GL}>1$. Each vortex has a current-carrying area
of radius $\lambda$ around its core. At the vortex core, the modulus of the
order parameter is suppressed. The exact definition of $\xi$ is the
characteristic exponent of the gap decay far away from the center of the
vortex core. At low temperatures, the overall size of the vortex core is
smaller than $\xi$ in the simplest models [6]. Note also that in strongly
type-2 superconductors, the long-range asymptotic of a vortex is affected by
nonlinearities [7]. When $\kappa<1$, in an ordinary Ginzburg-Landau theory,
the vortex is too energetically expensive and thermodynamically unstable. At
$\kappa=1$ vortices do not interact [8, 9, 10] 111Going beyond Ginzburg-Landau
theory, at low temperatures, one obtains microscopic corrections that result
in small nonmonotonic inter-vortex forces in a narrow range o
$\kappa^{GL}\approx 1.1$ [43, 44, 45, 46, 47, 48]. It should be noted that it
is a small microscopic model-specific effect beyond Ginzburg-Landau theory,
i.e. it cannot be described by carrying Ginzburg-Landau expansion to higher
order in $\tau=(1-T/T_{c})$ [49]..

Ginzburg and Landau’s work [1] classified single-component superconductors.
Today, many of the superconducting states of interest break multiple
symmetries, for example, featuring a breakdown of time-reversal symmetry or
translation symmetry or nematicity. Hence, they require description in terms
of multiple order parameter fields $|\Psi_{i}|e^{i\theta_{i}}$ and hence must
be characterized by multiple coherence lengths. It was pointed out in [12]
that in multicomponent systems, a new regime is possible where some coherence
lengths are shorter than the magnetic field penetration length, and some are
larger: $\xi_{1}<\xi_{2}<...\lambda<\xi_{n}<\xi_{n+1}...$. This regime was
termed type-1.5 in [13]. The vortex excitations there can be viewed as
composite objects. Namely, they are bound states of elementary vortices with
phase winding only in one of the components $\oint dl\nabla\theta_{n}=\pm
2\pi$. Such elementary constituents each carry a fraction of flux quantum and
are much more energetically expensive than integer-flux vortices [14]. Then,
in an external field, the system is expected to form composite vortices where
all components have phase winding around the common core, which can also be
viewed as bound states of fractional vortices where fractions add up to one
flux quantum. The fractional vortices with phase winding only in the single
band have been recently experimentally observed [15]. When some of the
coherence lengths are larger than the magnetic field penetration length, the
density-density interaction results in long-range attractive intervortex
forces, while the magnetic- and current- interaction gives short-range
repulsion [12, 16, 17]. Consequently, in a low magnetic field, such a system
exhibits vortex clustering and phase separation in vortex droplets and the
Meissner domains [12, 16, 17, 18, 19]. The concept of type-1.5
superconductivity was also generalized to other systems beyond
superconductivity, such as the typology of Quantum Hall systems [20] and
neutron stars [21].

In the above, we emphasized the case of multiple broken symmetries. The more
nontrivial case for typology is represented by the commonly occurring
superconducting materials with multiband electronic structure. Such systems
have multiple superconducting gaps forming on different bands
$|\Delta_{\alpha}|e^{i\theta_{\alpha}}$ but where interband Josephson
interaction breaks symmetry down to single $U(1)$. In that case, symmetry does
not guarantee multiple correlation lengths, even when there are multiple
bands. Nonetheless, in the simplest Landau theories, breaking a symmetry
explicitly does not prohibit the existence of extra coherence length. It was
discussed at the level of two-band Ginzburg-Landau theory in [16, 17, 18, 22]
that in multiband $U(1)$ systems, multiple coherence lengths arise and are
associated with different linear combinations of the gaps fields 222Note that
in general, the inclusion of anisotropies creates multiple scales in the
magnetic field penetration lengths too [50]. However, the justification of
multiple coherence lengths is nontrivial: the two-band Ginzburg-Landau model
is an expansion in multiple small parameters associated with multiple small
gaps and corresponding gradients. Such an expansion is not always justified as
it is based not on a small parameter guaranteed by symmetry but depends on the
structure of the intercomponent interaction [24], such as strength and
presence or absence of frustration. The conditions and parameter range where
two coherence lengths occur in two-band $U(1)$ systems were studied in
microscopic quasi-classical Eilenberger formalism in [25, 24], confirming, at
the level of quasiclassical theory the existence of length scale hierarchy
$\xi_{1}<\sqrt{2}\lambda<\xi_{2}$ in two-band systems that break only a single
symmetry. In the simplest two-band models, it requires having weak interband
coupling [25, 24]. These works also calculated asymptotic intervortex forces
in two-band Eilenberger formalism. Solutions for vortex clusters in the
type-1.5 regime were obtained in several microscopically derived Ginzburg-
Landau models [24, 26, 27]. However, to date, no solutions for vortex clusters
in the type-1.5 regime were demonstrated in microscopic models. Vortex
clusters were observed experimentally in several multiband systems and
attributed to type-1.5 physics in [13, 28, 29, 30, 31, 32, 33, 34, 35].

A microscopic approach that retains even the shortest-length-scales physics is
the Bogoliubov-de-Gennes (BdG) formalism [36]. Within this formalism, fully
microscopic solutions, including self-consistent calculation of magnetic
field, were obtained for an isolated Abrikosov vortex in [37, 6, 38]. However,
obtaining the solutions for vortex clusters is significantly more challenging
as one cannot rely on an axially symmetric Anzats. Here, we report solutions
of vortex clusters in the fully self-consistent numerical treatment of the
multiband Bogoliubov–de- Gennes model, including a self-consistent solution
for the magnetic field.

The two-band BdG model, that we consider, is defined on a two-dimensional
square lattice, described by the mean-field Hamiltonian

$\begin{gathered}\mathcal{H}=-\sum_{\sigma\alpha}\sum_{<ij>}e^{\mathrm{i}qA_{ij}}c^{\dagger}_{i\sigma\alpha}c_{j\sigma\alpha}+\\\
+\sum_{i\alpha}\left(\Delta_{i\alpha}c^{\dagger}_{\uparrow
i\alpha}c^{\dagger}_{\downarrow
i\alpha}+{\rm{H.c.}}\right)+\frac{1}{2}F_{m}\,.\end{gathered}$ (1)

Here $<ij>$ denotes all nearest neighbor pairs, $c_{i\sigma\alpha}$ is the
fermionic annihilation operator at position $i$, with spin $\sigma$
($\sigma\in\\{\uparrow,\downarrow\\}$) and band index $\alpha$
($\alpha\in\\{1,2\\}$) and ${\rm{H.c.}}$ denotes Hermitian conjugation. The
phase factor $\exp{\mathrm{i}qA_{ij}}$ accounts for the interaction with the
magnetic vector potential $\bm{A}$ through Peierls substitution [39, 40],
$\frac{1}{2}F_{m}$ is the magnetic field energy density.

$A_{ij}=\int_{j}^{i}\bm{A}\cdot\differential{\bm{\ell}}\,,$ (2)

the gap fields are defined as

$\Delta_{i\alpha}=\sum_{\beta}V_{\alpha\beta}\expectationvalue{c_{\uparrow
i\beta}c_{\downarrow i\beta}}\,,$ (3)

where $V_{\alpha\beta}=V_{\beta\alpha}^{*}$ stands for quadratic interaction
term, and the fermionic current is

$J_{ij}=-2q\sum_{\alpha\sigma}\mathrm{Im}\left(\langle
c_{i\alpha\sigma}^{\dagger}c_{j\alpha\sigma}\rangle
e^{\mathrm{i}qA_{ij}}\right)\,.$ (4)

and discrete version of Maxwell’s equation
$\nabla\times\nabla\times\bm{A}=\bm{J}$ determines the connection between
$A_{ij}$ and $J_{ij}$.

The free energy associated with the tight binding Hamiltonian Eq. (1), may be
expressed as

$\begin{gathered}F_{H}=\sum_{i}\bm{\Delta}_{i}^{\dagger}V^{-1}\bm{\Delta}_{i}-k_{\rm{B}}T\Tr\ln\quantity(\mathrm{e}^{-\beta\mathcal{H}}+1)+\\\
+\frac{1}{2}\sum_{\mathrm{plaquettes}}B^{2},\end{gathered}$ (5)

where the magnetic field $B=\nabla\times\bm{A}$ is defined on plaquettes.

The self-consistency equations Eq. (3), Eq. (4), along with the Maxwell
equation, are solved numerically using iterative scheme, described in [41].
Two independent codes were used to validate the solutions. New values are
obtained for the vector potential and the gaps during each iteration, using
Eq. (3). They can be calculated by obtaining the eigenvectors $c_{\sigma
i\alpha}$ by direct diagonalization of the Hamiltonian or using the Chebyshev
spectral expansion scheme. The key results of the paper are obtained by
Graphic Processing Units (GPU)-based exact matrix diagonalization, which is
used for free energy calculation Eq. (5). In addition, we report approximate
solutions for a larger vortex cluster, obtained using the approximate
Chebyshev spectral expansion method for larger vortex clusters. The iteration
procedure stops when the convergence criteria $\absolutevalue{\delta
p/\left(p+\epsilon\right)}<\varepsilon$ is achieved for each of the parameters
$\Delta_{1}$, $\Delta_{2}$, and $A$ simultaneously.

Below, we report microscopic solutions for vortex clusters in the BdG model.
We demonstrate (i) the existence of multiple correlation lengths in the
microscopic solution, (ii) that these length scales form the required
hierarchy: $\xi_{1}<\lambda<\xi_{2}$, (iii) the intervortex interaction
potential has a minimum at a finite distance (iv) a multi-quanta vortex
separates into a bound state of single-quanta vortices, forming a cluster.

We study a square sample with linear size $L$ and open boundary conditions.
The crucial aspect is to avoid mesoscopic effects on vortex physics. For that
reason, the vortices are initiated with an initial guess far away from sample
boundaries, and the sample is chosen to be significantly large to avoid vortex
escape due to boundary attraction. For the same reason, only the regime with
moderate disparity of the length scales could be studied to have all the
characteristic length scales much smaller than the grid size and, at the same
time, significantly larger than lattice spacing. The calculations based on the
Chebyshev approximation method were performed on a grid with $L=64$. For free
energy calculations, we used the exact diagonalization method with double
precision, so the system size is decreased to $L=48$ sites. Since the
quasiclassical analysis [25] suggests that interband coupling should be very
weak to have well-defined multiple correlation lengths, we analyze the
Hamiltonian Eq. (1) with the following numerical parameters $q=0.6$,
$V_{11}=2.8$, $V_{22}=2.2$, $V_{12}=0.01$, and $T=0.264$. We use the the
convergence criteria $\epsilon=10^{-8}$, $\varepsilon=10^{-6}$.

First, we analyze the structure of a single vortex. Single-vortex states were
calculated using the Chebyshev spectral expansion approximation. From these
solutions, the asymptotics for magnetic field and order parameter correlation
length were obtained FIG. 1.

In the presence of even a small inter-band Josephson coupling $V_{12}$,
coherence lengths are affected quantitatively and qualitatively. Calculations
in two-band Ginzburg-Landau [16, 17] and Eilenberger [25, 24] formalisms
predict that away from a vortex, each gap field approaches its asymptotic
value $\left|\Delta_{1,2}^{\mathrm{u}}\right|$ with two length scales. In
these models, the gap asymptotic is given by a combination of two modified
Bessel functions:

$\displaystyle\left|\Delta_{1}\left(r\right)\right|=\left|\Delta^{\mathrm{u}}_{1}\right|-q_{1}\cos\Theta
K_{0}\left(r/\xi_{1}\right)+q_{2}\sin\Theta K_{0}\left(r/\xi_{1}\right)$ (6)
$\displaystyle\left|\Delta_{2}\left(r\right)\right|=\left|\Delta^{\mathrm{u}}_{2}\right|-q_{1}\sin\Theta
K_{0}\left(r/\xi_{1}\right)-q_{2}\cos\Theta K_{0}\left(r/\xi_{2}\right)$

This implies that due to interband coupling, the coherence lengths $\xi_{1,2}$
are associated with the linear combination of the gap fields rather than
individual bands.

The crossection on FIG. 1, shows two length scales in the gap function. From
the solution, we can see that even a weak interband coupling forces the same
long-range asymptotic on both gap functions despite very different behavior of
$\Delta_{1}$ near the origin. The solutions can be fitted with Eq. (6), which
gives us $\xi_{1}\approx 1.27$, $\xi_{2}\approx 4.53$, $q_{1}\approx 0.043$,
$q_{2}\approx-0.48$, $\Theta\approx 0.495\pi$, meanwile fit for
$B=q_{b}e^{-x/\lambda}\sqrt{\lambda/x}$ gives $\lambda\approx 2.25$ and
$q_{b}=0.14$. Therefore, it shows that $\xi_{1}<\lambda<\xi_{2}$, which
verifies, that the system is in the type-1.5 regime for these coupling
constants.

Figure 1: The absolute values for superconducting gap and magnetic field for
single vortex, showing three distinct length scales. Left - asymptotic of
single-vortex solution. Dots represent numerical data for the gap amplitudes
$\Delta_{1,2}$ and magnetic field $B$; lines show the fit for the function Eq.
(6). Right - heatmaps for $\left|\Delta_{1}\right|$, $\left|\Delta_{2}\right|$
and $\left|B\right|$. Dot lines show the cross-section line, represented on
the left figure.

Upon establishing the type-1.5 hierarchy of the length scales, we demonstrate
that the system forms vortex clusters due to the competition of long-range
attractive core-core interaction set by coherence length $\xi_{2}$ and
shorter-range current-current and magnetic interaction set by the magnetic
field penetration length $\lambda$. In the next step, we analyze the
interaction energy between two vortices. To ensure that the results are not
artifacts of numerical grid pinning, the calculations are performed for a
variety of initial conditions. The total free energy of the system Eq. (5) is
orders of magnitude higher than the intervortex interaction energy, so the
exact diagonalization method and double precision were used to ensure the
accuracy of the energy calculations.

If the distance between the vortices in the initial guess is large enough, the
interaction force will not move the vortex from its initial position due to
the exponential falloff of interaction, and the presence of an underlying
lattice. This on the one hand, allows us to calculate directly the energy of
vortex interaction on different distances away from the minimum of the
potential FIG. 2, on the other hand the grid pinning mean that we slightly
overstimate the interaction energy, compared to the analytic expression (8).
The interaction energy of vortex pair $E_{\mathrm{int}}\left(r\right)$ is
calculated as follows: We calculate numerically the energy of a vortex pair
put the positions $\textbf{r}_{1}$ and $\textbf{r}_{2}$:
$F_{\mathrm{p}}\left(L,\textbf{r}_{1},\textbf{r}_{2}\right)$ from this we
subtract (i) the energy of the system in the absence of vortices
$F_{\mathrm{u}}\left(L\right)$, and (ii) the energies of the solutions for
single vortices, which are calculated at the same positions as vortices in the
pair $\textbf{r}_{1,2}$, in order to diminish the finite-size and
discretization effects.

$E_{\mathrm{int}}\left(r\right)=F_{\mathrm{p}}\left(L,\textbf{r}_{1},\textbf{r}_{2}\right)+F_{\mathrm{u}}\left(L\right)-F_{\mathrm{v}}\left(L,\textbf{r}_{1}\right)-F_{\mathrm{v}}\left(L,\textbf{r}_{2}\right)$
(7)

where $r=\absolutevalue{\textbf{r}_{1}-\textbf{r}_{2}}$.

The long-range asymptotic form for the vortex-vortex interaction energy in
type-1.5 superconductor calculated in continuum two-band Ginzburg-Landau [16,
24] and Eilenberger [25] formalisms has the form

$V\left(r\right)=\alpha\left(q_{b}^{2}K_{o}\left(r/\lambda\right)-q_{2}^{2}K_{o}\left(r/\xi_{2}\right)-q_{3}^{2}K_{o}\left(r/\xi_{1}\right)\right)$
(8)

where $\alpha$ is a positive constant 333In more general models the prefactors
can have additional parameters [17] and more terms can be present [22, 50].

The two dominant terms with $\lambda$ and $\xi_{2}$ in this expression give a
minimum of the interaction potential at a certain intervortex distance that
depends on the competing coherence and magnetic field penetration length
scales. We can compare this approximate expression with the results of our
BdG-based calculations. First, since we are interested only in long-range
forces, we omit the part with the shortest length scale
$~{}K_{o}\left(r/\xi_{1}\right)$ from Eq. (6). In the above, we extracted the
coherence and magnetic field penetration lengths from our solutions. Using
these lengths, we fit the results from Eq. (7) using Eq. (8) and obtain
$\alpha$. Although, strictly speaking, the form Eq. (8) is derived using
different continuum models, that equation approximately fits the calculated
intervortex potential. It also shows that the long-range attraction is
dominated by density-density interaction between the extended vortex cores and
short-range repulsion is dominated by current-current and magnetic
interactions.

Figure 2: Left - vortex pair interaction energy Eq. (7), expressed in single-
vortex energy units $E_{\mathrm{v}}=F_{\mathrm{v}}-F_{\mathrm{u}}$ vs relative
distance $r/L$ between the vortex cores. The leftmost data point corresponds
to the calculated minimum of the interaction energy. The dashed line is the
analytic fit based on - Eq. (8). Right - $\left|\Delta_{1}\right|$,
$\left|\Delta_{2}\right|$ and $\left|B\right|$ for the minimal energy
solution. There is a noticeable core overlap in $\Delta_{2}$.

Finally, an approximate solution for a larger vortex cluster is obtained. The
exact diagonalization methods are computationally expensive, and larger system
sizes are necessary to study multi-vortex clusters, so the Chebyshev spectral
expansion method was used. The system Eq. (1) is evaluated for a single multi-
quanta vortex solution as an initial conditions. The giant vortex is not a
stable configuration and it decays into a vortex cluster. On FIG. 3 is shown
the solution for an initial guess of the giant vortex with charge 3.

Figure 3: Cluster made of 3 vortices. The panels show the absolute values
$\absolutevalue{\Delta_{1,2}}$ and phases $\varphi_{1,2}$ for gap
distributions, magnetic field $B$ and current $\absolutevalue{J}$ from Eq. (4)
in the system and tunneling conductance below $\sigma\left(0\right)$ and above
$\sigma\left(V_{1}\right)$ the gap. The initial state of the system was a
single vortex with charge 3. There is a visible overlap of the vortex cores in
the second component. The distance between the cores for each pair of vortices
is the same as for the minimal energy solution on FIG. 2

In conclusion, we demonstrated, in a fully microscopic formalism, the
superconductivity beyond type-1/type-2 dichotomy in two-band superconductors.
Such superconductors break only a single symmetry due to interband coupling.
Nonetheless, when the interband coupling is weak, the obtained vortex
solutions clearly show the effects of multiple correlation lengths. The
numerical solutions show that these correlation lengths are hybridized: i.e.
associated with different linear combinations of the gap fields. For weak
interband coupling, we find hierarchy of the length scales
$\xi_{1}<\lambda<\xi_{2}$ that leads to attractive intervortex interaction at
large separation due to core-core overlap. It does not exclude large disparity
of coherence lengths at stronger interband coupling, for example in the cases
of frustrated interband coupling or proximity to a phase transition into a
superconducting state with different symmetry [22, 26]. The Bogoliubov-de-
Gennes formalism allows us also to calculate signatures of the vortex clusters
in Scanning Tunneling Microscopy. In such a probe, a vortex cluster in the
considered microscopic model, can appear as a group of vortices with
individual small cores yet having a significant attractive interaction. This
is because the large coherence length is associated with an asymptotic tail
where both gaps are suppressed at a long range but with a very small
amplitude.

###### Acknowledgements.

We thank A. Talkachov for insightful discussions and for sharing his exact
diagonalization CUDA code. This work was supported by the Swedish Research
Council Grants 2022-04763, by Olle Engkvists Stiftelse and the Wallenberg
Initiative Materials Science for Sustainability (WISE) funded by the Knut and
Alice Wallenberg Foundation.

## References

  * Landau and Ginzburg [1950] L. Landau and V. Ginzburg, Zh. Eksp. Teor. Fiz 20, 546 (1950).
  * Gor’kov [1959] L. P. Gor’kov, Sov. Phys. JETP 9, 1364 (1959).
  * Abrikosov [1957] A. A. Abrikosov, Sov. Phys.-JETP (Engl. Transl.);(United States) 5 (1957).
  * Rjabinin and Shubnikow [1935] J. N. Rjabinin and L. Shubnikow, Nature 135, 581 (1935).
  * Shubnikov  _et al._ [1937] L. Shubnikov, V. Khotkevich, Y. D. Shepelev, and Y. N. Ryabinin, Zh. Eksp. Teor. Fiz 7, 221 (1937).
  * Gygi and Schluter [1990a] F. Gygi and M. Schluter, Phys. Rev. Lett. 65, 1820 (1990a).
  * Plohr [1981] B. Plohr, Journal of Mathematical Physics 22, 2184 (1981).
  * Kramer [1971] L. Kramer, Phys. Rev. B 3, 3821 (1971).
  * Saint-James  _et al._ [1969] D. Saint-James, G. Sarma, and E. J. Thomas,  _TYPE-II SUPERCONDUCTIVITY._ , Tech. Rep. (CEN, Saclay, France, 1969).
  * Bogomol’nyi and Vainshtein [1976] E. Bogomol’nyi and A. Vainshtein, Sov. J. Nucl. Phys.(Engl. Transl.);(United States) 23 (1976).
  * Note [1] Going beyond Ginzburg-Landau theory, at low temperatures, one obtains microscopic corrections that result in small nonmonotonic inter-vortex forces in a narrow range o $\kappa^{GL}\approx 1.1$ [43, 44, 45, 46, 47, 48]. It should be noted that it is a small microscopic model-specific effect beyond Ginzburg-Landau theory, i.e. it cannot be described by carrying Ginzburg-Landau expansion to higher order in $\tau=(1-T/T_{c})$ [49].
  * Babaev and Speight [2005] E. Babaev and M. Speight, Phys. Rev. B 72, 180502 (2005).
  * Moshchalkov  _et al._ [2009] V. Moshchalkov, M. Menghini, T. Nishio, Q. H. Chen, A. V. Silhanek, V. H. Dao, L. F. Chibotaru, N. D. Zhigadlo, and J. Karpinski, Phys. Rev. Lett. 102, 117001 (2009).
  * Babaev [2002] E. Babaev, Physical review letters 89, 067001 (2002).
  * Iguchi  _et al._ [2023] Y. Iguchi, R. A. Shi, K. Kihou, C.-H. Lee, M. Barkman, A. L. Benfenati, V. Grinenko, E. Babaev, and K. A. Moler, Science 380, 1244 (2023).
  * Babaev  _et al._ [2010] E. Babaev, J. Carlström, and M. Speight, Phys. Rev. Lett. 105, 067003 (2010).
  * Carlström  _et al._ [2011] J. Carlström, E. Babaev, and M. Speight, Phys. Rev. B 83, 174509 (2011).
  * Carlström  _et al._ [2011] J. Carlström, J. Garaud, and E. Babaev, Phys. Rev. B 84, 134515 (2011).
  * Díaz-Méndez  _et al._ [2017] R. Díaz-Méndez, F. Mezzacapo, W. Lechner, F. Cinti, E. Babaev, and G. Pupillo, Phys. Rev. Lett. 118, 067001 (2017).
  * Parameswaran  _et al._ [2012] S. Parameswaran, S. Kivelson, E. Rezayi, S. Simon, S. Sondhi, and B. Spivak, Phys. Rev. B 85, 241307 (2012).
  * Wood and Graber [2022] T. S. Wood and V. Graber, Universe 8, 228 (2022).
  * Carlström  _et al._ [2011] J. Carlström, J. Garaud, and E. Babaev, Phys. Rev. B 84, 134518 (2011).
  * Note [2] Note that in general, the inclusion of anisotropies creates multiple scales in the magnetic field penetration lengths too [50].
  * Silaev and Babaev [2012] M. Silaev and E. Babaev, Phys. Rev. B 85, 134514 (2012).
  * Silaev and Babaev [2011] M. Silaev and E. Babaev, Phys. Rev. B 84, 094515 (2011).
  * Garaud  _et al._ [2018] J. Garaud, A. Corticelli, M. Silaev, and E. Babaev, Phys. Rev. B 98, 014520 (2018).
  * Agterberg  _et al._ [2014] D. F. Agterberg, E. Babaev, and J. Garaud, Physical Review B 90, 064509 (2014).
  * Nishio  _et al._ [2010] T. Nishio, V. H. Dao, Q. Chen, L. F. Chibotaru, K. Kadowaki, and V. V. Moshchalkov, Phys. Rev. B 81, 020506 (2010).
  * Ray  _et al._ [2014] S. J. Ray, A. S. Gibbs, S. J. Bending, P. J. Curran, E. Babaev, C. Baines, A. P. Mackenzie, and S. L. Lee, Phys. Rev. B 89, 094504 (2014).
  * Biswas  _et al._ [2020] P. Biswas, F. N. Rybakov, R. Singh, S. Mukherjee, N. Parzyk, G. Balakrishnan, M. Lees, C. Dewhurst, E. Babaev, A. Hillier, _et al._ , Phys. Rev. B 102, 144523 (2020).
  * Bolotina  _et al._ [2022] N. Bolotina, O. Khrykina, A. Azarevich, N. Y. Shitsevalova, V. Filipov, S. Y. Gavrilkin, K. Mitsen, V. Voronov, and N. Sluchanko, Phys. Rev. B 105, 054511 (2022).
  * Zhang  _et al._ [2022a] A.-L. Zhang, L.-X. Gao, J.-Y. He, V. B. Filipov, S. Cao, Q.-L. Xiao, and J.-Y. Ge, Science China Physics, Mechanics & Astronomy 65, 297412 (2022a).
  * Bolotina  _et al._ [2023] N. Bolotina, O. Khrykina, A. Azarevich, N. Y. Shitsevalova, V. Filipov, S. Y. Gavrilkin, A. Y. Tsvetkov, S. Gabáni, K. Flachbart, V. Voronov, _et al._ , Solid State Sciences 142, 107245 (2023).
  * Curran  _et al._ [2023] P. Curran, S. Bending, A. Gibbs, and A. Mackenzie, Scientific Reports 13, 12652 (2023).
  * Zhang  _et al._ [2022b] A.-L. Zhang, L.-X. Gao, J.-Y. He, V. B. Filipov, S. Cao, Q.-L. Xiao, and J.-Y. Ge, Science China Physics, Mechanics & Astronomy 65, 297412 (2022b).
  * de Gennes [2018] P. de Gennes,  _Superconductivity of metals and alloys_ (CRC press, 2018).
  * Gygi and Schluter [1990b] F. Gygi and M. Schluter, Phys. Rev. B 41, 822 (1990b).
  * Gygi and Schlüter [1991] F. Gygi and M. Schlüter, Phys. Rev. B 43, 7609 (1991).
  * Peierls [1933] R. Peierls, Zeitschrift für Physik 80, 763 (1933).
  * Feynman  _et al._ [2011] R. P. Feynman, R. B. Leighton, and M. Sands,  _The Feynman lectures on physics, Vol. I: The new millennium edition: mainly mechanics, radiation, and heat_ , Vol. 1 (Basic books, 2011).
  * Benfenati  _et al._ [2023] A. Benfenati, M. Barkman, and E. Babaev, Phys. Rev. B 107, 094503 (2023).
  * Note [3] In more general models the prefactors can have additional parameters [17] and more terms can be present [22, 50].
  * Jacobs [1973] A. Jacobs, J. Low Temp. Phys. 10, 137 (1973).
  * Leung and Jacobs [1973] M. C. Leung and A. E. Jacobs, J. Low Temp. Phys. 11, 395 (1973).
  * Eilenberger and Büttner [1969] G. Eilenberger and H. Büttner, Zeitschrift fur Physik 224, 335 (1969).
  * Weber  _et al._ [1987] H. Weber, M. Botlo, F. Sauerzopf, H. Wiesinger, and U. Klein, Japanese journal of applied physics 26, 917 (1987).
  * Klein [1987] U. Klein, J. Low Temp. Phys. 69, 1 (1987).
  * Miranović and Machida [2003] P. Miranović and K. Machida, Physical Review B 67, 092506 (2003).
  * Rybakov and Babaev [2021] F. N. Rybakov and E. Babaev, Journal of Mathematical Physics 62, 121901 (2021).
  * Winyard  _et al._ [2019] T. Winyard, M. Silaev, and E. Babaev, Phys. Rev. B 99, 064509 (2019).


# Thermodynamics to infer the astrophysics of binary black hole mergers

Patrick Hu  Department Physics and Astronomy, Vanderbilt University, 2301
Vanderbilt Place, Nashville, TN, 37235, USA  Karan Jani  Department Physics
and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN,
37235, USA<EMAIL_ADDRESS>Kelly Holley-Bockelmann  Department
Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place,
Nashville, TN, 37235, USA  Gregorio Carullo  Dipartimento di Fisica “Enrico
Fermi”, Università di Pisa, Pisa I-56127, Italy  INFN sezione di Pisa, Pisa
I-56127, Italy

###### Abstract

We introduce the Merger Entropy Index ($\mathcal{I}_{\mathrm{BBH}}$), a new
parameter to measure the efficiency of entropy transfer for any generic binary
black hole merger in General Relativity. We find that
$\mathcal{I}_{\mathrm{BBH}}$ is bounded between an asymptotic maximum and
minimum. For the observed population of mergers detected by LIGO [1] and Virgo
[2], we find that $\mathcal{I}_{\mathrm{BBH}}$ is $\lesssim 30\%$ of its
theoretical maximum. By imposing the thermodynamical consistency between the
pre- and post-merger states through $\mathcal{I}_{\mathrm{BBH}}$, we showcase
BRAHMA – a novel framework to infer the properties and astrophysical
implications of gravitational-wave detections. For GW190521 – the heaviest
confirmed binary black hole merger observed so far – our framework rules out
high mass-ratio, negative effective inspiral spin, and electromagnetic
counterpart claims. Furthermore, our analysis provides an independent
confirmation that GW190521 belongs to a separate population.

_Entropy transfer_.— A binary black hole (BBH) merger [3], regardless of its
astrophysical origins, is an irreversible process. From the Second Law of
Thermodynamics, the merger inevitably leads to an increase of entropy,
$\Delta{S}_{\mathrm{BBH}}$, in the Universe:

$\Delta{S}_{\mathrm{BBH}}\equiv S_{\mathrm{f}}-(S_{1}+S_{2})\geq 0,$ (1)

where $S_{\mathrm{f}}$ is the entropy of the remnant black hole formed after
merger, and $S_{1}$, $S_{2}$ are the entropies of the two pre-merger black
holes. According to the laws of black hole mechanics [4], this change in
entropy is a function of the mass and spin magnitudes of black holes before
$(m_{1,2},|\boldsymbol{\chi_{1,2}}|)$ and after merger
$(m_{\mathrm{f}},|\boldsymbol{\chi_{\mathrm{f}}}|)$. For the family of Kerr
black holes [5], the entropy can be measured as:

$S_{i}=\frac{2\pi
G}{c\hbar}m_{i}^{2}\left(1+\sqrt{1-|\boldsymbol{\chi_{i}}|^{2}}\right)\,.$ (2)

We construct a mass-invariant ratio, $\epsilon_{S}$, to quantify the entropy
transfer during the coalescence of BBH from its pre- to post-merger stage, and
note that it must be bracketed as follows:

$\epsilon_{S,\mathrm{min}}\leq\epsilon_{S}\equiv\frac{{S}_{\mathrm{f}}}{{S}_{\mathrm{1}}+{S}_{\mathrm{2}}}\leq\epsilon_{S,\mathrm{max}}\,,$
(3)

where $\epsilon_{S,\mathrm{max}}$ and $\epsilon_{S,\mathrm{min}}=1$ are the
global maximum and minimum of entropy transfer across the BBH parameter space
$\Lambda(q,\boldsymbol{\chi_{1}},\boldsymbol{\chi_{2}})$ with
$q=m_{2}/m_{1}\leq 1$ .

_Global Maximum for BBH Entropy_.— By inserting Eq. (2) in Eq. (3), we find
that entropy transfer has a leading order term equal to $(1+q)^{2}/(1+q^{2})$.
Therefore, the entropy transfer will be larger for equal-mass BBH mergers
$(q{=}1)$. Instead, it is in the extreme mass-ratio limit $(q\to 0)$ of BBH
mergers that $\epsilon_{S}\to 1$.

To acquire an understanding of the dependence of $\epsilon_{S}$ on
$\Lambda(q,\boldsymbol{\chi_{1}},\boldsymbol{\chi_{2}})$, and obtain an
estimate of $\epsilon_{S,\mathrm{max}}$, it is useful to first specialise to
mergers of comparable masses $(q{\simeq}1)$ and spins
$(|\boldsymbol{\chi_{1,2}}|\simeq\chi)$. In this case, the entropy transfer
can be approximated as:

$\epsilon_{S}\approx
2\left(1-E_{\mathrm{rad}}\right)^{2}\left(\frac{1+\sqrt{1-|\boldsymbol{J}-\boldsymbol{J_{\mathrm{rad}}}|^{2}}}{1+\sqrt{1-\chi^{2}}}\right)\,,$
(4)

where
$\boldsymbol{J}=\boldsymbol{\chi_{1}}+\boldsymbol{\chi_{2}}+\boldsymbol{L}$ is
the total angular momentum in the pre-merger system, $\boldsymbol{L}$ is the
corresponding orbital angular momentum, and $E_{\mathrm{rad}}$ and
$\boldsymbol{J_{\mathrm{rad}}}$ are the energy and angular momentum radiated
to infinity via gravitational waves. The radiated quantities and angular
momentum are expressed here in dimensionless units, and normalized with binary
mass $(M=m_{1}+m_{2})$.

Based on Eq. (4), if a BBH coalescence radiates more energy, it would lead to
a smaller change in entropy. If there is no residual angular momentum after
the merger, it would lead to a higher change in entropy. First noted by
Hawking [6], the theoretical maximum energy that can be radiated in the
comparable mass and spin case is $E_{\mathrm{rad}}=1/2$. A minimum entropy
transfer would require a binary with two maximally spinning black holes
$(\chi\to 1)$ to radiate all the angular momentum
$(\boldsymbol{{J}_{\mathrm{rad}}}~{}{\sim}~{}\boldsymbol{J})$. In this case,
Eq. (4) approaches $\epsilon_{S,\mathrm{min}}$. However, there exist no
configurations in the $\Lambda$ hyperspace satisfying such a conditions. In
fact, numerical relativity simulations find that BBH systems are unable to
radiate more than ${\sim 10\%}$ rest mass [7, 8, 9]. Thus, for near equal-mass
BBH systems, the entropy transfer is always greater than the lower-bound
imposed by the Second Law.

For generic quasi-circular BBH systems, $E_{\mathrm{rad}}$ and
$\boldsymbol{J_{\mathrm{rad}}}$ are strongly dependent on the mass-ratio $q$,
spin vectors $\boldsymbol{\chi_{1,2}}$ and their orientation with respect to
$\boldsymbol{L}$. The radiated quantities dictates the mass and spin of the
black hole formed after the merger [10], and consequently, the final entropy
$S_{\mathrm{f}}$. To find the global maximum for quasi-circular BBHs, we
numerically solve Eq. (3) across the $\Lambda$ hyperspace. We utilize the
phenomenological fits of $E_{\mathrm{rad}}$ and
$\boldsymbol{J_{\mathrm{rad}}}$ [11, 12, 13]. These fits are calibrated on
numerical relativity simulations of BBH mergers [7, 9], which incorporate the
apparent horizon framework [14]. We find good agreement on
$\epsilon_{S,\mathrm{max}}$ between the two different fits employed,
$\epsilon_{S,\mathrm{max}}=3.6235$ from [11], and
$\epsilon_{S,\mathrm{max}}=3.6342$ from [13]. As expected,
$\epsilon_{S,\mathrm{max}}$ occurs for a system of equal-mass, maximally
spinning black holes, aligned in the direction opposite that of the orbital
angular momentum.

For the special case of head-on collisions $(L\sim 0)$ of equal-mass BBHs, the
entropy transfer can, in principle, be even higher. Suppose a collision occurs
between two maximally spinning black holes that are anti-aligned with each
other $(J\sim 0)$. Provided the collision radiates no energy
$(E_{\mathrm{rad}}~{}{\sim}~{}0)$, only then, Eq. (4) would approach a
theoretical upper-bound of $\epsilon_{S,\mathrm{max}}=4$. However, based on
numerical relativity simulations, we can infer that such collisions would emit
${\gtrsim}1.35\%$ of rest mass. Thus, a subset of highly eccentric, maximally
spinning, equal-mass mergers can approach an entropy transfer in the range
$3.62\lesssim\epsilon_{S,\mathrm{max}}\lesssim 3.91$.

Figure 1: Global behavior of Merger Entropy Index
$(\mathcal{I}_{\mathrm{BBH}})$ within the intrinsic BBH parameter space,
$\Lambda$. For demonstration, this figure only shows non-precessing binaries
with equal spin vectors.

_Merger Entropy Index_.— By using the bounds of $\epsilon_{S}$, we construct
the Merger Entropy Index as,

$\boxed{\mathcal{I}_{\mathrm{BBH}}\equiv\frac{\pi}{9}\frac{\Delta{S}_{\mathrm{BBH}}}{{S}_{\mathrm{1}}+{S}_{\mathrm{2}}}}\,\,.$
(5)

This index measures the efficiency of entropy transfer during a BBH
coalescence. All BBH systems in General Relativity will satisfy
$\mathcal{I}_{\mathrm{BBH}}\in(0,1]$. By construction, the index is mass-
invariant and is thus applicable to a wide range of astrophysical scenarios.
For all astrophysical binaries, we anticipate
$\mathcal{I}_{\mathrm{BBH}}\lesssim 0.9$.

Fig. 1 shows the global behavior of Merger Entropy Index within a slice in the
hyperspace $\Lambda$, the intrinsic BBH parameters. We employ RIT fits [11] to
compute $\mathcal{I}_{\mathrm{BBH}}$, which are valid for quasi-circular to
moderately eccentric BBHs. If these BBHs are non-spinning, we find that
$\mathcal{I}_{\mathrm{BBH}}\propto q/(1+q)^{2}\equiv\eta$, the symmetric mass-
ratio.

For near equal-mass BBHs, the index has a rather strong dependence on the
orientation of spin vectors. When the effective inspiral spin parameter
$(\chi_{\mathrm{eff}})$ [15] is negative, $\mathcal{I}_{\mathrm{BBH}}$
exponentially approaches to 0.9. On the contrary, when $\chi_{\mathrm{eff}}$
is positive, $\mathcal{I}_{\mathrm{BBH}}$ has very limited variance of
${\sim}20\%$. This is consistent regardless of the mass-ratio of the system.
For high mass-ratio BBHs with high primary BH spin, the index exponentially
scales to the lower-bound of 0. We note that BBHs with high in-plane spins is
that they tend to have high $\mathcal{I}_{\mathrm{BBH}}$.

Figure 2: Cumulative distribution of Merger Entropy Index for a uniform prior
within the intrinsic BBH parameter space $\Lambda$ (red curve). The dotted
green curve is for region of $\Lambda$ derived from the observed population of
BBH by LIGO-Virgo [16]

In Fig. 2, we show the cumulative distribution of $\mathcal{I}_{\mathrm{BBH}}$
across the intrinsic BBH parameter space $\Lambda$. The red curve is induced
by a uniform distribution within mass-ratios $1/q\in[1,20]$, spin magnitudes
$\chi_{1,2}\in[0,0.99]$ and cosine of the spin orientations
$\cos{\theta_{1,2}}\in[-1,1]$. This induced distribution is a good
representation of the default prior on the Merger Entropy Index for quasi-
circular BBH mergers in General Relativity. When assuming the uniform prior,
$\mathcal{I}_{\mathrm{BBH}}\lesssim 0.5$ at 99% credibility, greatly
restricting its allowed range. This is consistent with Fig. 1, which shows
that higher values of $\mathcal{I}_{\mathrm{BBH}}$ are confined in a small
domain of $\Lambda$.

The dotted green curve in Figure 3 shows a distribution of
$\mathcal{I}_{\mathrm{BBH}}$ drawn out of an astrophysical prior. This prior
is derived directly from the constraints on the masses and spins for the
observed population of BBH mergers in GWTC-3 [16] from the LIGO-Virgo
detectors. In particular, we use the distribution $p(q)\propto q^{\beta}$ from
PowerLaw+Peak Model to fetch the mass-ratio, and the Gaussian Spin Model for
the distribution of black hole spin magnitudes and orientations. For the
astrophysical prior, we find $\mathcal{I}_{\mathrm{BBH}}$ is steeply centered
around 0.2. There are no values $\gtrsim 0.3$ with 99% confidence interval.

From the distributions discussed above, we can infer the Merger Entropy Index
has an upper-bound as:

$0\leq\mathcal{I}_{\mathrm{BBH}}\lesssim\begin{cases}1.0~{}~{}~{}\texttt{Theoretical
Limits; }\\\ 0.5~{}~{}~{}\texttt{Uniform Prior;}\\\
0.3~{}~{}~{}\texttt{Astrophysical Prior}\end{cases}$ (6)

While the theoretical and uniform prior bounds of $\mathcal{I}_{\mathrm{BBH}}$
in Eq. (6) are independent of the BBH origins, the astrophysical bounds could
well vary between the supermassive and stellar mass BBH population, and would
be interesting to explore in future work. We found that KS-convergence can
confidently differentiate (with a probability of one part in billion) between
the LIGO-Virgo derived distribution of $\mathcal{I}_{\mathrm{BBH}}$ [16] from
the uniform prior distribution. This suggests $\mathcal{I}_{\mathrm{BBH}}$ can
be an effective tool to test astrophysical formation channels.

Figure 3: Schematic diagram of BRAHMA Framework. We use a BBH detection
(wavelet at the top) whose pre- and post-merger states have been independently
measured. From this information we compute the Merger Entropy Index matrix. By
imposing the bounds of $\mathcal{I}_{\mathrm{BBH}}$ stated in Equation (6, we
filter parts of the inspiral posterior distribution (rows in the figure) which
pass this criterion with $\geq\rho\%$.  Figure 4: New constraints on
GW190521’s mass-ratio and effective spins using BRAHMA Left: Constraints using
the pre- and post-merger measurement from set-1, i.e. the LIGO-Virgo-KAGRA
Collaboration (LVK). The black curve refers to the published distribution from
the SEOBNRv4PHM model. The dashed green line shows the constraints we obtain
in BRAHMA by imposing theoretical bounds of $\mathcal{I}_{\mathrm{BBH}}$. The
filled green contour is the constraint from BRAHMA for a uniform prior bound
of $\mathcal{I}_{\mathrm{BBH}}$. Right: Constraints using pre- and post-merger
measurement from set-2, i.e. the Open Gravitational-Wave Catalog (OGC). The
black curve refers to the published distribution from the IMRPhenomXPHM model.
The red dashed and red contour refers to the constraints obtained from BRAHMA
by imposing theoretical and uniform prior bounds of
$\mathcal{I}_{\mathrm{BBH}}$ on this set. The two-dimensional curves show the
90% credible regions.

_Novel GW Framework_.— By adopting the Merger Entropy Index, we introduce a
generalized framework called BRAHMA (BinaRy blAck Holes Merger entropy
Analysis). The schematic flow of BRAHMA is shown in Fig. 3. Every BBH
detection whose pre- $(S_{1},S_{2})$ and post-merger $(S_{\mathrm{f}})$
entropies have been independently measured (e.g. because models with different
physical content have been used to infer them) have to satisfy the upper and
lower bounds of the theoretical limits of $\mathcal{I}_{\mathrm{BBH}}$ from
Eq. (6). This is under the assumption that the observed BBH population is
consistent with the predictions of General Relativity, which by now has a
strong experimental support [17, 18, 19]. With the current generation of
gravitational-wave data analysis tools, the pre-merger entropy will typically
be measured by conducting a Bayesian analysis over the intrinsic $\Lambda$
(masses, spins) and extrinsic $\Theta$ (distance, orientation, sky-location)
BBH parameter space with inspiral-merger-ringdown (IMR) waveform models (green
wavelet in Fig. 3). For the post-merger entropy, our framework strictly
requires an independent measurement of the black hole ringdown (pink wavelet
in Fig. 3), which is achieved by an independent measurement of the post-merger
properties [17, 18, 19]. See Ref. [20] for thermodynamical constraints on the
information emission process of LIGO-Virgo BHs using post-merger measurements.

From the posterior distributions of pre-merger entropies (with size $1\times
m$) and post-merger entropies ($1\times n$), we compute a Merger Entropy Index
matrix $(m\times n)$. Assuming the detected BBH is consistent with General
Relativity, every element in the matrix has to satisfy (pass) theoretical
bounds stated in Equation (6). For illustration, we show them as grey boxes of
$\mathcal{I}_{\mathrm{BBH}}$ in the matrix in Figure 3. The matrix may contain
some failed elements, i.e they do not obey the $\mathcal{I}_{\mathrm{BBH}}$
bounds (shown in the Fig. 3 as red boxes). For every row in the
$\mathcal{I}_{\mathrm{BBH}}$ matrix, we compute the the fraction of passed
columns to that of the total elements in the column $(n)$. If this fraction is
less than the threshold $\rho$, we discard this row vector and filter the pre-
merger distribution to a new length $1\times(m-i)$. As a default setting in
our framework, we fix $\rho=90\%$. BRAHMA directly filters the parameter space
constrained by the IMR model. All the BBH inspiral parameters, including
extrinsic $(\Theta)$ ones, are now further constrained than the original
output from the Bayesaian analysis, even though BRAHMA does not explicitly
require any direct information on most of these parameters (e.g. luminosity
distance, coalescence phase). Naturally, the constraints on BBH parameters are
tighter if the bounds on the Merger Entropy Index are set from a uniform or
astrophysical prior, compared to the most general theoretical limit.

The framework is built on the idea that for certain systems, post-merger
analyses provide a more robust measurement of the end state properties than
full-signal analyses. This situation can happen when full-signal models do not
correctly represent complete solutions of the Einstein equations (i.e. they
are ”unfaithful” to accurate numerical simulation of BBHs), or when they lack
the inclusion of certain physical effects. The post-merger stage is generally
simpler to model regardless of the pre-merger history of the BBH system, and
that the agnostic post-merger models employed in our analysis do not rely on
numerical calibration to quasi-circular mergers. For example, eccentricity
effects are much less relevant in the post-merger due to fast binary
circularisation in the late stages of the pre-merger dynamics, see Fig. 4 of
[21]. On the contrary, pre-merger models employed by the LVK collaboration
currently cannot incorporate astrophysically interesting scenarios such as
eccentric orbits, dynamical captures and high-mass ratio precessing binaries.
Some progress in some of these directions has been recently obtained using
Effective One Body models for non-precessing eccentric systems and eccentric
dynamical captures [22, 23]. However, such templates currently lack the
inclusion of BHs spins and higher harmonics.

_Applications to GW190521_.— To demonstrate the application of our BRAHMA
framework, we apply it to the parameters measurements of GW190521. This event
has been shown to suffer from systematic uncertainties of its pre-merger
properties [24, 25, 26, 23], thus falling in the category of systems to which
our framework is especially relevant. Detected during the third observing run
of the LIGO and Virgo detectors, this event marks the first direct measurement
of an intermediate-mass black hole [27], though its astrophysical formation
remains so far mysterious [28]. Such massive BBHs have an advantage of an
independent measurement of the remnant mass and spin
$(m_{\mathrm{f}},\chi_{\mathrm{f}})$ from the ringdown stage [29]. Several
studies have performed a full Bayesian inference on the pre- and post-merger
properties of GW190521 using a variety of gravitational waveform models and
numerical relativity simulations [28, 30, 31, 32, 33]. From these data
products, we chose two fully independent set of measurements to test BRAHMA:

  1. set-1:

We choose the waveform model SEOBNRv4PHM [34] to measure the pre-merger
entropies $(S_{1},S_{2})$. We do so using its posteriors on detector-frame
masses $(m_{1},m_{2})$ and spin magnitudes $(\chi_{1},\chi_{2})$, available
from the Gravitational-Wave Open Science Center (GWOSC) maintained by the
LIGO-Virgo-KAGRA Collaboration [35]. These pre-merger constraints were
obtained with the Bayesian parameter scheme RIFT [36]. For the post-merger
entropy $(S_{\mathrm{f}})$, we use the posteriors on the remnant mass and spin
$(m_{\mathrm{f}},\chi_{\mathrm{f}})$ obtained through the pyRing analysis
framework [37] with the Kerr_222 model [38, 28]. This template incorporates
the $\ell=m=2$ quasi-normal modes up to the second overtone ($n=0,1,2$), and
is applied starting at the peak of the waveform (a proxy for the merger).
Among the post-merger models employed in [28], it provides the most precise
measurements of ringdown, and thus that of remnant mass and spin.

  2. set-2:

From the Open Gravitational Catalog (OGC) [39], which is an independent
analysis on the open data, we get the pre-merger estimates for the
IMRPhenomXPHM waveform model [40]. Their analysis is performed using the
Bayesian code pyCBC Inference [41]. The post-merger estimates are obtained
from a Kerr template [42] incorporating the $(\ell,m,n)=\\{(2,2,0),(3,3,0)\\}$
modes, as there is claimed evidence for both these ringdown modes in GW190521
[43].

| set-1: | set-2:
---|---|---
| LVK | With | OGC | With
Parameter | SEOB-PHM | BRAHMA | IMR-XPHM | BRAHMA
Lower-bound on $q$ | $0.32$ | $0.39$ | $0.22$ | $0.37$
Lower-bound on $\chi_{\mathrm{eff}}$ | $-0.28$ | $-0.12$ | $-0.60$ | $+0.06$
Probability($z\leq$ AGN) | $10\%$ | $7\%$ | $41\%$ |  $5\%$
Table 1: Source properties for GW190521 with $\mathcal{I}_{\mathrm{BBH}}$
framework and imposing uniform prior. For mass-ratio $q$ and effective
inspiral spin $\chi_{\mathrm{eff}}$, we provide the lower-bound at the 90%
confidence interval. The third row is the cumulative probability for the
inferred redshift $z$ to be equal or less than the one from the claimed AGN
flare, J1249+349 [44].

As shown in the two panels of Fig. 4-5, the inference of pre-merger properties
by these two sets is significantly different (compare the black curves in
left-right panels). However, both the sets provide support for (i) an
asymmetric mass distribution $(q\lesssim 0.3)$, (ii) negative effective
inspiral spins $(\chi_{\mathrm{eff}}\lesssim-0.2)$ and (iii) a source redshift
that overlaps with an AGN flare, J1249+349 [44]. These features carry strong
implications for the astrophysical origins of GW190521. By applying our BRAHMA
framework, we find that both set-1 and set-2 no longer provides evidence for
the features (i)-(iii). Table 1 summarizes the overall improvement in pre-
merger properties of GW190521 using BRAHMA.

Figure 5: New constraints on the redshift of GW190521 using the BRAHMA
framework. The color coding on both the panels is same as Figure 4. The dotted
grey line refers to the redshift of the electromagnetic candidate observed by
the ZTF collaboration [44].

After applying the theoretical prior of $\mathcal{I}_{\mathrm{BBH}}$ in
BRAHMA, set-1 (LVK:SEOB-PHM) no longer shows the bimodality in
$q-\chi_{\mathrm{eff}}$ parameter space. For set-2 (OGC:IMR-XPHM), the
bimodality still remains, however, the support for
$\chi_{\mathrm{eff}}\lesssim 0$ is significantly lower (see dashed lines in
Fig. 4). These bounds also have significant impact on the source redshift for
set-2, dropping its cumulative probability at $z\leq 0.438$ (electromagnetic
candidate claim) by a factor of ${\sim}2$, confirming GW190521 is farther than
the claimed electromagnetic candidate.

When the uniform prior bounds are adopted in BRAHMA, we see a significant
reduction in the 90% confidence intervals in both the sets. The most notable
changes appear in our constraint from set-2 on the effective inspiral spin of
GW190521, which suggests the parameter is positive with 90% confidence. Also,
for the same set the cumulative probability for source redshift at $z\leq
0.438$ drops to ${\sim}5\%$ (a factor of 8). Filtered posteriors from both
sets no longer show any bimodality, and in fact, converge to similar a post-
merger distribution after applying our $\mathcal{I}_{\mathrm{BBH}}$ framework.
Moreover, the distribution from these two sets now seems consistent with other
waveform models that have analyzed GW190521 [28], suggesting a possible
improvement in waveform systematic with $\mathcal{I}_{\mathrm{BBH}}$.

By imposing the astrophysical prior in BRAHMA, we find that 99.9% (95%) of
pre-merger distribution is filtered out from set-2 (set-1). This provides an
additional confirmation that GW190521 is at the tail end of the astrophysical
population derived from all other LIGO-Virgo black holes [16].

_Conclusion_.— We introduce a new thermodynamical formalism into
gravitational-wave astrophysics. For the current generation detectors, we find
our formalism is particularly powerful for probing GW190521-like lower-range
intermediate-mass black holes $({\sim}10^{2-3})~{}M_{\odot}$ [45, 46]. Our
framework is applicable to BBH lying in regions of the parameter space where
full numerical simulations are not extensively available and for which the
pre- and post-merger states can be independently measured, for all mass
ranges. Our formalism will be even more germane for next-generation detectors
Cosmic Explorer [47], Einstein Telescope [48], upcoming space-mission LISA
[49] and potential deci-Hertz experiments [50, 51], as we expect the post-
merger measurements to be more often and louder. Future work will explore the
implications of the Merger Entropy Index $(\mathcal{I}_{\mathrm{BBH}})$ and
the BRAHMA framework on BBH demographics, formation channels and gravitational
waveform systematics.

##  I acknowledgments

We thank Maximiliano Isi and Nathan Johnson-McDaniel for useful comments.
P.H.’s research was supported by the Immersion Initiative at Vanderbilt
University. K.J and K.H.B.’s research was supported by the Stevenson Chair
endowment funds at Vanderbilt University.

This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. This research
has made use of data, software and/or web tools obtained from the
Gravitational Wave Open Science Center (https://www.gw-openscience.org/), a
service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo
Collaboration. LIGO Laboratory and Advanced LIGO are funded by the United
States National Science Foundation (NSF) as well as the Science and Technology
Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS),
and the State of Niedersachsen/Germany for support of the construction of
Advanced LIGO and construction and operation of the GEO600 detector.
Additional support for Advanced LIGO was provided by the Australian Research
Council. Virgo is funded, through the European Gravitational Observatory
(EGO), by the French Centre National de Recherche Scientifique (CNRS), the
Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef,
with contributions by institutions from Belgium, Germany, Greece, Hungary,
Ireland, Japan, Monaco, Poland, Portugal, Spain.

## References

  * LIGO Scientific Collaboration [2015] LIGO Scientific Collaboration, Classical and Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Virgo Collaboration [2015] Virgo Collaboration, Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * LIGO Scientific Collaboration and Virgo Collaboration [2016] LIGO Scientific Collaboration and Virgo Collaboration, Physical Review Letters 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Bardeen  _et al._ [1973] J. M. Bardeen, B. Carter, and S. W. Hawking, Communications in Mathematical Physics 31, 161 (1973).
  * Bekenstein [1973] J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * Hawking [1971] S. W. Hawking, Phys. Rev. Lett. 26, 1344 (1971).
  * Mroué  _et al._ [2013] A. H. Mroué, M. A. Scheel, B. Szilágyi, H. P. Pfeiffer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, S. Ossokine, N. W. Taylor, A. Zenginoğlu, L. T. Buchman, T. Chu, E. Foley, M. Giesler, R. Owen, and S. A. Teukolsky, Physical Review Letters 111, 241104 (2013).
  * Jani  _et al._ [2016] K. Jani _et al._ , Classical and Quantum Gravity 33, 204001 (2016).
  * Healy  _et al._ [2017] J. Healy, C. O. Lousto, Y. Zlochower, and M. Campanelli, Classical and Quantum Gravity 34, 224001 (2017).
  * Barausse and Rezzolla [2009] E. Barausse and L. Rezzolla, Astrophys.J. 704, L40 (2009), arXiv:0904.2577 [gr-qc] .
  * Healy  _et al._ [2014] J. Healy, C. O. Lousto, and Y. Zlochower, Phys.Rev. D90, 104004 (2014), arXiv:1406.7295 [gr-qc] .
  * Johnson-McDaniel  _et al._ [2016] N. K. Johnson-McDaniel, A. Gupta, D. Keitel, P. Ajith, O. Birnholtz, F. Ohme, and S. Husa, LIGO DCC , LIGO-T1600168-v6 (2016).
  * Varma  _et al._ [2019] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. Research 1, 033015 (2019).
  * Dreyer  _et al._ [2003] O. Dreyer, B. Krishnan, D. Shoemaker, and E. Schnetter, Phys. Rev. D 67, 024018 (2003).
  * Pan  _et al._ [2011] Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 84, 124052 (2011).
  * The LIGO Scientific Collaboration  _et al._ [2021a] The LIGO Scientific Collaboration, The Virgo Collaboration, and The KAGRA Scientific Collaboration, arXiv e-prints , arXiv:2111.03634 (2021a), arXiv:2111.03634 [astro-ph.HE] .
  * Carullo [2021] G. Carullo, Phys. Rev. D 103, 124043 (2021), arXiv:2102.05939 [gr-qc] .
  * Isi  _et al._ [2021] M. Isi, W. M. Farr, M. Giesler, M. A. Scheel, and S. A. Teukolsky, Physical Review Letters 127 (2021), 10.1103/physrevlett.127.011103.
  * Abbott  _et al._ [2021a] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. Adhikari, V. Adya, C. Affeldt, and et al., Physical Review D 103 (2021a), 10.1103/physrevd.103.122002.
  * Carullo  _et al._ [2021] G. Carullo, D. Laghi, J. Veitch, and W. Del Pozzo, Phys. Rev. Lett. 126, 161102 (2021), arXiv:2103.06167 [gr-qc] .
  * Huerta  _et al._ [2019] E. A. Huerta _et al._ , Phys. Rev. D 100, 064003 (2019), arXiv:1901.07038 [gr-qc] .
  * Nagar  _et al._ [2021] A. Nagar, P. Rettegno, R. Gamba, and S. Bernuzzi, Phys. Rev. D 103, 064013 (2021), arXiv:2009.12857 [gr-qc] .
  * Gamba  _et al._ [2021] R. Gamba, M. Breschi, G. Carullo, P. Rettegno, S. Albanesi, S. Bernuzzi, and A. Nagar, (2021), arXiv:2106.05575 [gr-qc] .
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw, P. D. Lasky, E. Thrane, and J. C. Bustillo, Astrophys. J. Lett. 903, L5 (2020), arXiv:2009.04771 [astro-ph.HE] .
  * Bustillo  _et al._ [2021a] J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forné, and J. A. Font, Phys. Rev. Lett. 126, 201101 (2021a), arXiv:2009.01066 [gr-qc] .
  * Bustillo  _et al._ [2021b] J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forné, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, Phys. Rev. Lett. 126, 081101 (2021b), arXiv:2009.05376 [gr-qc] .
  * LIGO Scientific Collaboration and Virgo Collaboration [2020] LIGO Scientific Collaboration and Virgo Collaboration, Phys. Rev. Lett. 125, 101102 (2020), arXiv:2009.01075 [gr-qc] .
  * Abbott  _et al._ [2020] R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. B. Adya, C. Affeldt, M. Agathos, and et al., The Astrophysical Journal 900, L13 (2020).
  * Carullo  _et al._ [2019] G. Carullo, W. Del Pozzo, and J. Veitch, Physical Review D 99 (2019), 10.1103/physrevd.99.123029.
  * Romero-Shaw  _et al._ [2020] I. Romero-Shaw, P. D. Lasky, E. Thrane, and J. Calderón Bustillo, Astrophys. J. Lett 903, L5 (2020), arXiv:2009.04771 [astro-ph.HE] .
  * Nitz and Capano [2021] A. H. Nitz and C. D. Capano, The Astrophysical Journal 907, L9 (2021).
  * Estellés  _et al._ [2021] H. Estellés, S. Husa, M. Colleoni, M. Mateu-Lucena, M. de Lluc Planas, C. García-Quirós, D. Keitel, A. Ramos-Buades, A. K. Mehta, A. Buonanno, and S. Ossokine, arXiv e-prints , arXiv:2105.06360 (2021), arXiv:2105.06360 [gr-qc] .
  * Gayathri  _et al._ [2020] V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. Lousto, and R. O’Shaughnessy, arXiv e-prints , arXiv:2009.05461 (2020), arXiv:2009.05461 [astro-ph.HE] .
  * Ossokine  _et al._ [2020] S. Ossokine, A. Buonanno, S. Marsat, R. Cotesta, S. Babak, T. Dietrich, R. Haas, I. Hinder, H. P. Pfeiffer, M. Pürrer, and et al., Physical Review D 102 (2020), 10.1103/physrevd.102.044055.
  * Abbott  _et al._ [2021b] R. Abbott _et al._ , SoftwareX 13, 100658 (2021b).
  * Lange  _et al._ [2018] J. Lange, R. O’Shaughnessy, and M. Rizzo, “Rapid and accurate parameter inference for coalescing, precessing compact binaries,” (2018), arXiv:1805.10457 [gr-qc] .
  * [37] G. C. et al., “pyRing: a time-domain ringdown analysis python package,” https://git.ligo.org/lscsoft/pyring.
  * Giesler  _et al._ [2019] M. Giesler, M. Isi, M. A. Scheel, and S. Teukolsky, Phys. Rev. X 9, 041060 (2019), arXiv:1903.08284 [gr-qc] .
  * Nitz  _et al._ [2021] A. H. Nitz, C. D. Capano, S. Kumar, Y.-F. Wang, S. Kastha, M. Schäfer, R. Dhurkunde, and M. Cabero, arXiv e-prints , arXiv:2105.09151 (2021), arXiv:2105.09151 [astro-ph.HE] .
  * Pratten  _et al._ [2021] G. Pratten, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, M. Mateu-Lucena, R. Jaume, M. Haney, D. Keitel, J. E. Thompson, and et al., Physical Review D 103 (2021), 10.1103/physrevd.103.104056.
  * Biwer  _et al._ [2019] C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, Publications of the Astronomical Society of the Pacific 131, 024503 (2019).
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021c), arXiv:2010.14529 [gr-qc] .
  * Capano  _et al._ [2021] C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz, A. B. Nielsen, and B. Krishnan, arXiv e-prints , arXiv:2105.05238 (2021), arXiv:2105.05238 [gr-qc] .
  * Graham  _et al._ [2020] M. Graham, K. Ford, B. McKernan, N. Ross, D. Stern, K. Burdge, M. Coughlin, S. Djorgovski, A. Drake, D. Duev, and et al., Physical Review Letters 124 (2020), 10.1103/physrevlett.124.251102.
  * Jani  _et al._ [2020] K. Jani, D. Shoemaker, and C. Cutler, Nature Astronomy 4, 260 (2020), arXiv:1908.04985 [gr-qc] .
  * The LIGO Scientific Collaboration  _et al._ [2021b] The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration, arXiv e-prints , arXiv:2105.15120 (2021b), arXiv:2105.15120 [astro-ph.HE] .
  * Reitze  _et al._ [2019] D. Reitze _et al._ , arXiv e-prints , arXiv:1907.04833 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ ,  _Proceedings, 14th Workshop on Gravitational wave data analysis (GWDAW-14): Rome, Italy, January 26-29, 2010_ , Class. Quant. Grav. 27, 194002 (2010).
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ , arXiv e-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Sato  _et al._ [2017] S. Sato _et al._ , Journal of Physics: Conference Series 840, 012010 (2017).
  * Jani and Loeb [2021] K. Jani and A. Loeb, Journal of Cosmology and Astroparticle Physics 2021, 044 (2021).