∎

11institutetext: Antonios Tsokaros 22institutetext: University of Illinois at
Urbana-Champaign, Loomis Laboratory, 1110 W Green St, Urbana, IL 61801, USA.
Tel.: +1-217-3004508
Fax : +1-217-3334898
22email<EMAIL_ADDRESS>33institutetext: Kōji Uryū 44institutetext:
University of the Ryukyus, Department of Physics, Senbaru, Nishihara, Okinawa
903-0213, Japan.
Tel. & Fax: +81988958521
44email<EMAIL_ADDRESS>
# Methods for relativistic self-gravitating fluids: From binary neutron stars
to black hole-disks and magnetized rotating neutron stars. ††thanks:
A.Tsokaros is supported by National Science Foundation Grant No. PHY-1662211
and the National Aeronautics and Space Administration (NASA) Grant No.
80NSSC17K0070 to the University of Illinois at Urbana-Champaign. K. Uryū is
supported by JSPS Grant-in-Aid for Scientific Research (C) 15K05085 and
18K03624 to the University of Ryukyus.

Antonios Tsokaros  Kōji Uryū

(Received: date / Accepted: date)

###### Abstract

The cataclysmic observations of event GW170817, first as gravitational waves
along the inspiral motion of two neutron stars, then as a short $\gamma$-ray
burst, and later as a kilonova, launched the era of multimessenger astronomy,
and played a pivotal role in furthering our understanding on a number of
longstanding questions. Numerical modelling of such multimessenger sources is
an important tool to understand the physics of compact objects and, more
generally, the physics of matter under extreme conditions. In this review we
present a unified view of various techniques used to obtain equilibrium and
quasiequilibrium solutions for three astrophysically relevant relativistic,
self-gravitating fluid systems: Binary neutron stars, black hole-disks, and
magnetized rotating neutron stars. These solutions are necessary not only for
modeling such compact objects, but equally important, for providing self-
consistent initial data in numerical relativity simulations. Instead of
presenting the full details of the formulations and numerical algorithms, we
focus on painting the broadbrush picture of the methods developed to address
these problems, and facilitate future work in the area.

This article is dedicated to the memory of Yoshiharu Eriguchi, whose work on
self-gravitating fluids has been a source of inspiration.

##  1 Introduction

Although the majority of neutron stars are observed as isolated pulsars that
emit electromagnetic radiation from their magnetic poles Hewish1968 ;
lyne_graham-smith_2012 , a small percentage of them appear in binary systems.
In our own Galaxy the first binary neutron star (BNS) system, PSR B1913+16,
was discovered in 1974 HulseTaylor74  while currently 19 such binaries are
known Zhu2020 ). Despite the small number of BNS currently observed their
strong gravitational and electromagnetic interactions constitutes them as
ideal probes for relativistic astrophysical phenomena. In particular:

  1. (i)

During their inspiral and merger they produce gravitational waves an important
feature of the general theory of relativity. Fundamental tests related to
gravity in an otherwise inaccessible strong-field regime are therefore
possible MillerColeman2016 ; Abbott2019 .

  2. (ii)

The late inspiral as well as the merger process can yield important
information about the masses and radii of the component stars and in effect
about their equation of state (EOS) Yagi2013b ; lorimer:lr ; Lattimer2016 ;
Ozel2016 .

  3. (iii)

They can be the progenitors of short duration ($<2\;\rm s$) $\gamma$-ray
bursts (sGRBs) Eichler89 ; Narayan92 ; Fong2013 .

  4. (iv)

They constitute premier sites for the production of rare heavy elements like
gold and platinum Lattimer74 ; Lattimer76 ; Symbalisty82 , through a rapid
neutron-capture process (or r-process), in which neutrons are captured by
lighter nuclei like iron in a dense neutron-rich environment Burbidge1957 ;
Cameron1957 .

In a single stroke the first detection by LIGO/Virgo of a BNS approximately
$40\;\rm Mpc$ away in the Galaxy NGC 4993, together with follow-up detections
by the Fermi, INTEGRAL spacecrafts and other telescopes, addressed all of the
problems above. Event GW170817 showed the inspiral of two neutron stars
GW170817prl  which was followed by a sGRB, named GRB170817A $1.7\;\rm s$ later
Goldstein2017 ; Savchenko2017 ; Abbott2017d . Lastly a luminous optical
counterpart, named AT2017gfo, was also observed 11 hours after the merger
Coulter2017 ; Soares-Santos2017 ; Arcavi2017 ; Diaz2017 ; Abbott2017c . This
kilonova Metzger:2010  transient was powered by the radioactive decay of heavy
neutron-rich elements created in the expanding merger ejecta Li:1998  and
served as a direct probe of the astrophysical origin of the heaviest elements
in the Universe (see recent review by Metzger Metzger2019 ). Despite the
scarcity of detections like event GW170817, it is apparent that the synthesis
of gravity, electromagnetism, and microphysical processes in binary black
hole-neutron star (BHNS) or BNS mergers (binary black holes are not expected
to produce electromagnetic radiation due to the absence of matter) is
astrophysically very fruitful ShibataTaniguchilrr-2011-6 ; Berger2013b ;
Baiotti2016 ; Paschalidis2016 ; Shibata2019a ; Meszaros2019 ; Friedman2020 ;
Foucart2020 ; Ciolfi2020 ; Radice2020 .

The background spacetime in which many of the complex phenomena mentioned
above develop, is the strong gravity regime. Depending on the question at
hand, or the timescale, one can consider two categories of problems. In the
first category there is a dominant compact object (either single or composite)
that sets up the spacetime and the rest of the system is evolving either
without affecting the gravitational sector or by perturbing it. An example of
such system is a black hole surrounded by a non self-gravitating disk
Abramowicz2011 . A widely separated BNS where the neutron stars are treated as
point masses, also falls within this class of problems Blanchet06 . The
complement of the first category constitutes of all systems where the
gravitational interactions between the various components are equally
significant, cannot be isolated, and vary with time. In this category one
necessarily ends up with a dynamical spacetime where the self-gravity of every
component is important. Analytical work for these kind of systems is limited,
and full 3+1 (spatial and temporal) computational calculations are the main
tool of investigation. In this endeavor the Einstein’s equations are playing
the leading role since both electromagnetic and nuclear interactions are
acting within the unknown spacetime. Therefore numerical simulations that try
to explain complex phenomena behind BHNS or BNS mergers inevitably have to be
able to integrate in a stable way the equations of general relativity, a
second order system of partial differential equations. In addition in order to
incorporate magnetic fields one needs Maxwell equations, and finally the
energy-momentum and radiation transport equations for the evolution of matter
and radiation.

Any differential equation needs initial values in order to be integrated and
the aforementioned dynamical system is no exception. Even more, due to the
complexity of the problem there are two main difficulties. The first one is
mathematical and is associated with the fact that the Einstein dynamical
system itself, therefore its initial values too, is not trivially posed.
Afterall the Einstein equations for the pair $(M,g_{\alpha\beta})$, where $M$
is a 4-dimensional manifold, and $g_{\alpha\beta}$ a Lorentzian metric, are
invariant under diffeomorphisms of $M$ and the associated isometries of
$g_{\alpha\beta}$. In addition even if a formulation is mathematically well
posed it does not mean that it will be numerically stable. The second
difficulty is rooted into the fact that the initial value functions must
represent the physical system under consideration which is a complicated task
in the general theory of relativity. Basic concepts like the mass, angular
momentum or center of mass are not trivially defined as they are in Newtonian
mechanics. In other words it is not always clear what kind of assumptions one
has to make in order to get a snapshot of the system under consideration which
will be physically meaningful. The subject of numerical relativity and
relativistic hydrodynamics has grown considerably since the pioneering binary
black hole simulations Pretorius:2005gq ; Campanelli06 ; Baker:2005vv  and
details can be found in a number of textbooks Wilson2003 ; Alcubierre:2008 ;
Bona2009 ; Baumgarte2010 ; gourgoulhon20123+1 ; Friedman2012 ;
Rezzolla_book:2013 , while recent reviews Lehner2014 ; Tichy2016rev ;
Paschalidis2017b ; Duez2019 ; Radice2020  present different aspects of it.

In this review we will touch upon the initial value problem and its numerical
implementation of three astrophysically relevant self-gravitating systems that
include matter: BNSs, black hole-disks (BHD), and magnetized rotating neutron
stars (MRNSs). For reasons mentioned above the modeling of a neutron star in a
binary or single setting as well as its magnetic field is very important not
only for gravitational wave astronomy but also multimessenger astrophysics. In
addition since the most promising scenario for the existence of a sGRB is a
black hole surrounded by a massive disk Eichler89 ; Narayan92 ; Shibata99d ;
Baiotti08 ; Anderson2007 ; Liu:2008xy ; Bernuzzi2011 ; Paschalidis2014 ;
Ruiz2016  the study of such systems is also well-motivated. Furthermore BHD
are ubiquitous in the Universe and self-gravity can be important at certain
times during their evolution. Although BHD have been intensely studied in the
past Frank2002 ; Kato2008 ; Abramowicz2011 ; Rezzolla_book:2013 , here we will
focus only on general relativistic self-gravitating disks that are not covered
in the bibliography, where a fixed Kerr black hole is assumed. Apart from the
BNS initial value problem which is discussed extensively in the bibliography,
BHD and MRNS are still missing a concrete exposition, therefore our effort
here is to close this gap and offer a unified approach to the subject. Also we
do not discuss the initial value problem of binary black holes, BHNSs or
rotating neutron stars (RNS) since it is sufficiently covered in the textbooks
and the reviews mentioned above. Due to limitations in space, our discussion
will focus only on the salient features and methods of obtaining self-
gravitating solutions of generic BNS, BHD, and MRNS (quasi)equilibria, without
getting into the details of those calculations. We will narrow our exposition
only in full general relativistic methods. Neither we will address the
important subject of stability and evolution in general for these compact
objects that inevitably will make this review grow manyfold. Our main goal is
to offer a bird’s eye view of the subject, while at the same time motivate
further research in the field of numerical general relativistic solutions.

In this review Greek indices are taken to run from 0 to 3 while Latin indices
from 1 to 3. We use a signature $(-,+,+,+)$ for the spacetime line element,
and a system of units in which $c=G=M_{\odot}=1$ (unless explicitly shown).
The list of acronyms adopted in this paper are listed below:

BNS: | Binary neutron star |  | sGRB: | Short $\gamma$-ray burst
---|---|---|---|---
RNS: | Rotating neutron star |  | EOS: | Equation of state
MRNS: | Magnetized rotating neutron star |  | ADM: | Arnowitt-Deser-Misner
BHD: | Black hole-disk |  | KEH: | Komatsu-Eriguchi-Hachisu
MHD: | Magnetohydrodynamics |  | IWM: | Isenberg-Wilson-Mathews
IMHD: | Ideal magnetohydrodynamics |  |  |

##  2 Binary neutron stars

In general relativity contrary to Newtonian gravity binary compact objects
evolve by the emission of gravitational waves. When two neutron stars are
separated by a distance much larger than their radii one can approximate them
as point particles Blanchet2013 . At that stage the neutron stars evolve in an
adiabatic manner with the gravitational wave timescale being larger than the
orbital period $P_{\rm orb}$,

$\frac{t_{\rm gw}}{P_{\rm
orb}}\approx\left(\frac{r}{6M}\right)^{5/2}\left(\frac{M}{4\mu}\right).$ (1)

Here $M$ is the total mass of the binary and $\mu$ its reduced mass ($\mu=M/4$
for equal mass binaries). When the separation becomes smaller, the
gravitational wave timescale decreases, as does the orbital period. When the
two become comparable, the radial velocity of the neutron stars increases
significantly and the adiabatic approximation breaks down. For a typical
system this happens at $35\;\rm km$ or $1\;\rm KHz$ gravitational wave
frequency with the orbital period being approximately $2\;\rm ms$. The neutron
stars merge shortly afterwards. At the intermediate stage when on one hand the
two neutron stars are not too far (distances less than $60\;\rm km$), but on
the other not too close, (distances greater than $35\;\rm km$) the system can
be approximated as stationary in the corotating frame and finite size effects
are important. The combined Eistein-Euler system needs to be resolved for an
accurate representation of the system. The methods described below aim at that
stage. The solutions obtained can be used as accurate initial values for
performing full general relativistic simulations and study the late inspiral,
merger and postmerger at the nonlinear regime, but they can also be used on
their own for gaining important information regarding this stage in the
evolution of the binary.

###  2.1 Isenberg-Wilson-Mathews formulation

One of the pillars in constructing binary neutron star initial data is the so-
called Isenberg-Wilson-Mathews (IWM) formulation Isenberg08 ; Wilson89 ;
Wilson95 ; Wilson96 . Three waveless formulations were proposed by Isenberg in
1978 Isenberg08  which simplify computations of astrophysical systems of
compact objects by decoupling the gravitational wave part and an “induced”
part of strong gravity (which may be associated with the matter source terms).
Isenberg never implemented the formulation into a numerical code, but later
Wilson and Mathews Wilson89 ; Wilson95 ; Wilson96  implemented one of the
waveless formulations, in which the spatial metric is assumed to be
conformally flat, for the evolution of BNS.

Since a binary system is inherently nonaxisymmetric, all 10 metric components
in the spacetime line element are necessary. A convenient way to express that
is by the use of the 3+1 form

$ds^{2}=-\alpha^{2}dt^{2}\,+\,\gamma_{ij}(dx^{i}+\beta^{i}dt)(dx^{j}+\beta^{j}dt)\,,$
(2)

where $\alpha,\beta^{i}$ are the lapse function and shift vector respectively,
and $\gamma_{ij}$ the spatial metric. A conformal 3-metric,
$\tilde{\gamma}_{ij}$, and a conformal traceless extrinsic curvature,
$\tilde{A}_{ij}$, are introduced through Lichnerowicz44 ; Nakamura1994

$\gamma_{ij}=\psi^{4}\tilde{\gamma}_{ij},\qquad\mbox{and}\qquad
A_{ij}=\psi^{4}\tilde{A}_{ij},$ (3)

where $\psi$ is the confomal factor. In the IWM formulation, the spatial
metric is conformally flat $\tilde{\gamma}_{ij}=f_{ij}$ ($f_{ij}=\delta_{ij}$
in Cartesian coordinates), the slicing is maximal $K=0$, and the time
derivatives of the conformal metric vanish
$\partial_{t}\tilde{\gamma}_{ij}=0$. Under such assumptions, part of the
Einstein system reduces to 5 elliptic equations for the conformal factor
$\psi$, the lapse $\alpha$ and the shift $\beta^{i}$. These are the
Hamiltonian constraint, the spatial trace of $\partial_{t}K_{ij}$, and the
momentum constraint. In the IWM formulation the tracefree part of the
extrinsic curvature is written in terms of the lapse, the shift, and the
conformal factor as

$A^{ij}=K^{ij}=\frac{\psi^{-4}}{2\alpha}\left({\raise 4.30554pt\hbox{${}^{\
\circ}$}}\\!\\!\\!\\!\\!D^{i}\beta^{j}+{\raise 4.30554pt\hbox{${}^{\
\circ}$}}\\!\\!\\!\\!\\!D^{j}\beta^{i}-\frac{2}{3}f^{ij}{\raise
4.30554pt\hbox{${}^{\ \circ}$}}\\!\\!\\!\\!\\!D_{k}\beta^{k}\right)\,,$ (4)

where ${\raise 4.30554pt\hbox{${}^{\ \circ}$}}\\!\\!\\!\\!\\!D$ is the
covariant derivative with respect to the flat metric, ${\raise
4.30554pt\hbox{${}^{\ \circ}$}}\\!\\!\\!\\!\\!D_{c}f_{ab}=0$ (in Cartesian
coordinates ${\raise 4.30554pt\hbox{${}^{\
\circ}$}}\\!\\!\\!\\!\\!D_{i}=\partial_{i}$). This scheme was first used in an
evolutionary study of BNS by Wilson et al. Wilson95 ; Wilson96 . Soon after,
it was realized that the IWM formulation was even more useful for the
construction of accurate initial data sets for full numerical relativity
simulations. The IWM formulation is still used in simulations of binary
neutron stars or binary white dwarf mergers in order to incorporate (part of)
relativistic gravity by replacing the Newtonian gravitational potential with
the above metric potentials Bauswein2010 ; Bauswein2014 . The connection
between gravity and matter can be accomplished in many ways and the IWM
formalism is one of them. A better formulation will be presented later in
Section 2.9 Bonazzola:2003dm ; Shibata04a .

###  2.2 Mass, angular momentum and the first law for binary systems

Two important characteristics of a BNS system are its mass-energy content, and
its angular momentum. In an asymptotically flat spacetime like the one
representing an isolated BNS a definition for a global mass was presented by
Arnowitt-Deser-Misner (ADM) Arnowitt1960 ; Arnowitt62unfindable . It is now
called the ADM mass

$M=\frac{1}{16\pi}\oint_{S_{{}_{\infty}}}(f^{ai}f^{bk}-f^{ki}f^{ab}){\raise
4.30554pt\hbox{${}^{\ \circ}$}}\\!\\!\\!\\!\\!D_{k}\gamma_{ab}dS_{i},$ (5)

where the integral is performed on a sphere whose radius tends at infinity.
Similarly the ADM angular momentum associated with a rotational Killing vector
$\phi^{i}$ is

$J=-\frac{1}{8\pi}\oint_{S_{{}_{\infty}}}\pi^{i}_{\ j}\phi^{j}dS_{i},$ (6)

were $\pi^{ij}=-(K^{ij}-\gamma^{ij}K)\sqrt{\gamma}$ the momentum conjugate to
$\gamma_{ij}$. Having a spacetime metric in the form of Eq. (2), one can
rewrite Eqs. (5) and (6) in terms of $\alpha,\psi,\beta^{i},K_{ij}$ (see e.g.
gourgoulhon20123+1 ) and therefore have a measure of the mass and angular
momentum in terms of the 3+1 variables.

According to the first law of thermodynamics for binary systems by Friedman et
al. Friedman01a , if one assumes a spatial geometry $\Sigma_{t}$ that is
conformally flat, neighboring equilibria of asymptotically flat spacetimes
with a helical Killing vector satisfy

$\delta M=\Omega\delta J+\int_{\Sigma_{t}}[\bar{T}\Delta dS+\bar{\mu}\Delta
dM_{0}+V^{\alpha}\Delta
d\mathcal{C}_{\alpha}]+\sum_{i}\frac{1}{8\pi}\kappa_{i}\delta A_{i}.$ (7)

Here $M$ and $J$ are the ADM mass and angular momentum of the spacetime while
$\Omega$ is the orbital angular velocity; $\bar{T}$ and $\bar{\mu}$ are the
redshifted temperature and chemical potential; $dM_{0}$ is the baryon mass of
a fluid element; $d\mathcal{C}_{\alpha}$ is related to the circulation of a
fluid element, and $V^{\alpha}$ is the velocity with respect to the corotating
frame; $\kappa_{i},\ A_{i}$ are the surface gravity, and the areas of black
holes. For isentropic fluids, dynamical evolution conserves the baryon mass,
entropy, and vorticity of each fluid element, and thus the first law yields

$\delta M=\Omega\delta J.$ (8)

Eq. (7) implies that a natural measure to characterize the spin of a neutron
star in a binary setting is its circulation in a similar manner to the way
rest mass characterizes the mass.

###  2.3 Equations for perfect fluids

In this review the stress-energy tensor for the matter will be described by a
perfect fluid with 4-velocity $u^{\alpha}$, pressure $p$ and total energy
density $\epsilon$,

$T^{\alpha\beta}_{\rm\scalebox{0.6}{M}}=(\epsilon+p)u^{\alpha}u^{\beta}+pg^{\alpha\beta}\,.$
(9)

Its conservation leads to

$\nabla_{\beta}T^{\alpha\beta}_{\rm\scalebox{0.6}{M}}=\rho[u^{\beta}\nabla_{\beta}(hu^{\alpha})+\nabla^{\alpha}h]+hu^{\alpha}\nabla_{\beta}(\rho
u^{\beta})-\rho T\nabla^{\alpha}s=0,$ (10)

where the first law of thermodynamics $dh=Tds+dp/\rho$ has been used for the
pressure gradient. Here $\rho,\,h,\,s$ are the rest-mass density, specific
enthalpy ($h:=(\epsilon+p)/\rho$), and specific entropy respectively.

The inspiral of two cold neutron stars can be considered as an isentropic
process therefore the last term in Eq. (10) can be set to zero. In addition
the system conserves rest mass, $\nabla_{\alpha}(\rho u^{\alpha})=0$,
therefore one arrives at the relativistic Euler equation,
$u^{\beta}\nabla_{\beta}(hu_{\alpha})+\nabla_{\alpha}h=0$.

Decomposing the 4-velocity with respect to the corotating observer as

$u^{\alpha}=u^{t}(k^{\alpha}+V^{\alpha}),$ (11)

with

$k^{\alpha}=t^{\alpha}+\Omega\phi^{\alpha}$ (12)

being the helical vector, one can rewrite the conservation of rest mass as

$\mathcal{L}_{\bf k}(\rho u^{t})+\nabla_{\alpha}(\rho u^{t}V^{\alpha})=0,$
(13)

and the Euler equation as,

$\mathcal{L}_{\bf
k}(hu_{\alpha})+\nabla_{\alpha}\left(\frac{h}{u^{t}}+hu_{\beta}V^{\beta}\right)+V^{\beta}\omega_{\alpha\beta}=0\,,$
(14)

where
$\omega_{\alpha\beta}=\nabla_{\beta}(hu_{\alpha})-\nabla_{\alpha}(hu_{\beta})$
is the vorticity tensor. Here $\mathcal{L}_{\bf k}$ is the Lie derivative
along $k^{\alpha}$, and vector $\phi^{\alpha}$ generates rotations in the
$\phi$ direction. Eqs. (13) and (14) are valid for any type of flow field. In
the next sections we will describe 3 types of flows: corotating, irrotational,
and arbitrary spinning.

###  2.4 Corotating binary neutron star initial data

The first BNS quasiequilibrium initial data have been constructed by Baumgarte
et al. Baumgarte97 ; Baumgarte98b  using the IWM formalism. Neglecting
deviations from a strictly periodic circular orbit and assuming the two stars
to be corotating, the fluid 4-velocity is proportional to the helical vector
$k^{\alpha}$, as in a single rotating neutron star, which is assumed to be a
Killing vector – a time translation symmetry in a rotating frame. Here
$\Omega$ in Eq. (12) is the orbital angular velocity of the binary system
which is also the spin angular velocity of the corotating component stars.
Similarly to the single rotating star case the Euler equations can be
integrated to yield a first integral. Equation (14) with $V^{\alpha}=0$ and
$\mathcal{L}_{\bf k}(hu_{\alpha})=0$ results to

$\frac{h}{u^{t}}=C\,,$ (15)

where $C$ is a constant. This first integral can be used to compute the single
unknown thermodynamic variable (for example $h$). The component $u^{t}$ can be
found from the normalization of the 4-velocity and the gravitational
variables, $u^{t}=1/\sqrt{\alpha^{2}-\omega^{i}\omega_{i}}$, where
$\omega^{i}=\beta^{i}+\Omega\phi^{i}$ is the so-called corotating shift
vector. There are two constants that appear in the equations, which are the
orbital angular velocity $\Omega$ and the constant $C$ that has the meaning of
the injection energy Thorne1967 . These are evaluated from two conditions, one
is to fix $\rho$ at the center of the star, and the other to fix the
separation. 111 Sometimes the number of constants is augmented to 3 or 4, and
additional conditions are needed. This results to a system of 3 or 4 non-
linear equations that needs to be solved at every iteration. The whole
procedure above is repeated until the differences in all computed quantities
between two subsequent iteration steps drop below a certain threshold error.

Baumgarte et al. Baumgarte97 ; Baumgarte98b  built quasiequilibrium binary
sequences, i.e. sequences of solutions of constant rest mass and decreasing
separation, which approximate evolutionary trajectories of neutron star
binaries undergoing slow inspiral via the generation of gravitational
radiation. By construction these sequences maintained corotation (i.e. the
spin frequency of the neutron stars increases as they become closer) which is
not realistic, since tidal torques due to realistic viscosity mechanisms are
not strong enough to synchronize the neutron star spin. The conserved quantity
in a BNS system besides the rest mass is circulation Kochanek92 ; Bildsten92 ;
Friedman01a . Nevertheless, constructing corotating sequences was a very
important first step in understanding close binary dynamics. They found that
the maximum density of the component stars decreases as they approach each
other and become tidally deformed 222The behavior of maximum density as the
binary system merges was the subject of intense debate which stemmed from
relativistic numerical simulations Wilson95 ; Wilson96 ; Mathews97  in which
it has been noted that as the stars approach each other their interior maximum
density increases. As a consequence, depending on the EOS, BNS would collapse
individually toward black holes prior to merger. At that time many authors
have disputed the finding of density increase Lai96 ; Wiseman97 ; Brady97 ;
Shibata98b ; Baumgarte98c , while Flanagan Flanagan:1998c  pointed to an
inconsistency in the solution of the shift vector employed in Wilson95 ;
Wilson96 ; Mathews97  and was responsible for this erroneous behavior
Mathews00 . . At the same time, the mass that a given maximum density can
support increases as the stars approach each other. The authors found that the
maximum allowed mass of neutron stars in quasiequilibrium binaries increases
with decreasing separation. These effects are larger for a smaller polytropic
index and hence a stiffer EOS. In this review we use the term stiff in two
ways. If the EOS is a simple polytrope, then a stiff EOS has a higher value of
polytropic exponent $\Gamma$. On the other hand if we have a realistic EOS (as
those represented in piecewise poytropic form Read:2009a ) a stiff EOS
signifies that a more massive neutron star can be supported333A lot of times
this means that for a given mass, stiff EOSs result to larger radius.
Confusion may arise since for simple polytropes (in terms of normalized mass
and radius), a higher $\Gamma$ result to neutron stars with smaller radius
(for a given mass), or to a smaller maximum mass, although the maximum
compactness $M/R$ becomes higher..

Baumgarte et al. computed the binding energy of the system
$E_{b}=M-M_{\infty}$, where $M$ is the ADM mass at a given separation, and
$M_{\infty}=M_{1}+M_{2}$ the sum of the two gravitational masses of the two
component stars at infinite separation. In a corotating binary as the stars
approach each other, the angular velocity increases while the binding energy
decreases as in Newtonian gravity. When the separation becomes sufficiently
small, finite size effects start playing an important role, and in particular
the (positive) rotational energy of the stars reduces the negative binding
energy. If the EOS is sufficiently stiff ($\Gamma\gtrapprox 2$ so that the
moment of inertia is sufficiently large) the binding energy goes through a
minimum and then increases again prior to contact. This minimum in energy
which coincides with the minimum in the angular momentum of the system (see
Eq.(8) Friedman01a ), approximately signifies the innermost stable circular
orbit (ISCO) beyond which a rapid plunge and merger occurs. By locating the
turning points in their total energy versus separation curves, the authors
identified the onset of orbital instability, and the orbital parameters at
that critical radius finding no evidence of destabilization Wilson95 ;
Wilson96 ; Mathews98 . For stiff EOSs the stars are more extended and such
minimum in the energy is absent i.e. the binding energy is monotonically
decreasing and merger happens at some minimum distance at the inner Lagrange
point.

Corotating BNSs have also been calculated by Marronetti et al. Marronetti98
using the IWM scheme similarly to Baumgarte97 , where the focus was in the
direct determination of the ISCO. The authors presented equal mass corotating
binaries from a larger separation, where the Newtonian approximation is valid,
to a smaller separation close to ISCO. Further analysis of conformally flat
corotating configurations has been presented by Miller et al. Miller:2003vc ,
who explicitly showed that if one takes into account the spin energy of the
neutron stars, the effective binding energy (i.e. the binding energy of an
effective irrotational BNS) no longer exhibits a minimum (this was in
agreement with irrotational quasiequilibrium BNS sequences as we describe in
Sec 2.5). The authors made a complete comparison of the quasiequilibrium
corotating sequences with full general relativistic simulations and analyzed
the conformal flatness assumption, and the assumption of the 4-velocity being
proportional to the Killing vector field $k^{\alpha}$ in a BNS evolution. It
turned out that the violation of the assumption regarding the timelike helical
Killing vector field was an order of magnitude larger than the violation of
the assumption of conformal flatness. In particular, for corotating BNSs at
separations less than $47M_{0}$ (which was slightly more than 6 neutron star
radii), where $M_{0}$ is the neutron star rest mass, the quasiequilibrium
solutions violate the Einstein field equations at the order of $10\%$. On the
other hand the conformal flatness assumption was violated at the order of
$1\%$.

###  2.5 Irrotational binary neutron star initial data

As mentioned earlier corotation is not maintained during the evolution of a
BNS due to their small viscosity Bildsten92 ; Kochanek92 ; Miller:2003vc .
Instead the conserved quantity is circulation, and it is expected that most
BNSs can be approximated as irrotational, i.e. having no vorticity in the
inertial frame. In this case the stars are counterotating in the rotating
frame and therefore the fluid velocity is nonzero there. This means that in
contrast to the corotating case, the continuity equation is nontrivial and
leads to an equation for the divergence of the fluid velocity. A first
formulation that treated the case of nonsynchronized binaries was presented by
Bonazzola et al. Bonazzola97  and Asada Asada1998 . Later Shibata Shibata98
and Teukolsky Teukolsky98  presented a simplified formulation introducing the
relativistic velocity potential. Since for irrotational BNSs the three
formalisms are equivalent (see Appendix A in Gourgoulhon01 ), with the latter
two being much simpler, and becoming the workhorse in all subsequent numerical
implementations, we will limit our discussion to those. As in the corotating
case, we assume that the system preserves its properties under the action of
the helical Killing vector (12) (stationarity property), i.e.

$\mathcal{L}_{\bf k}(u_{\alpha})=\mathcal{L}_{\bf k}(h)=0.$ (16)

For an irrotational flow the vorticity tensor $\omega_{\alpha\beta}$ is by
definition zero which implies that the 4-velocity can be derived from a
potential $\Phi$

$hu_{\alpha}=\nabla_{\alpha}\Phi,$ (17)

which together with the conservation of rest mass, Eq. (13), will yield an
elliptic equation for $\Phi$,
$\nabla_{\alpha}(\frac{\rho}{h}\nabla^{\alpha}\Phi)=0$, or in 3+1 form under
helical symmetry (16), (13)

$D_{i}(\alpha\rho u^{t}V^{i})=0,\qquad\mbox{with}\qquad
V^{i}=\frac{D^{i}\Phi}{hu^{t}}-\omega^{i}\,.$ (18)

Here $D_{i}$ is the 3D covariant derivative associated with $\gamma_{ij}$. A
boundary condition for this Poisson type of equation can be derived by setting
the fluid velocity in the corotating frame $V^{\alpha}$ at the stellar
surface, $r=R_{s}$, to be tangent to the surface itself,
$V^{\alpha}\nabla_{\alpha}p|_{r=R_{s}}=0$, which yields

$(D^{i}\Phi-hu^{t}\omega^{i})D_{i}\rho|_{r=R_{s}}=0.$ (19)

Having determined the enthalpy vector (canonical momentum), $hu_{\alpha}$, one
needs to calculate the enthalpy $h$ in order to have a complete solution of
the fluid equations for a barotropic EOS. Equation (14) under the assumptions
of Eqs. (16) yields the first integral

$\frac{h}{u^{t}}+V^{\beta}\nabla_{\beta}\Phi={\rm const}.$ (20)

We emphasize here that in order to arrive to the first integral Eq. (20), and
the elliptic equation for $\Phi$, Eq. (18), both the first (in particular
$\mathcal{L}_{\bf k}u^{t}=0$) and the second parts of Eq. (16) are used. In
other words both the thermodynamic profile and the velocity field need to
respect the approximation of stationarity in the corotating frame.

The irrotational formulation by Shibata Shibata98  and Teukolsky Teukolsky98
together with the IWM formalism Isenberg08 ; Wilson95 ; Wilson96  has been
used by many groups around the world in order to obtain realistic BNS initial
data. The first such calculation has been done by Bonazzola, et al.
Bonazzola98b  which explicitly showed that the maximum density of constant
rest-mass irrotational binaries decreases with respect to its value at
infinite separation, as the binaries are approaching each other. As expected
the decrease is not as pronounced as the corresponding decrease of the
corotating binaries. This can be explained by the fact that the latter have
considerable spin at close distances which leads to a larger decrease in their
central density, as in a sequence of single rotating neutron stars that starts
at the spherical limit and evolves towards the mass-shedding limit. These
first results have soon been confirmed by Marronetti et al. Marronetti99  as
well as by Uryū & Eriguchi Uryu00 . The picture that emerged from these early
BNS initial data studies Bonazzola98b ; Gourgoulhon-etal-2000:2ns-initial-data
; Marronetti99 ; Uryu00 ; Uryu00a ; Usui2000 ; Taniguchi02b ; Faber:2002zn ;
Taniguchi03 ; Usui2002 ; Bejger04  can be summarized as follows: 1)
Irrotational binary neutron star systems have been found to be dynamically
stable during the inspiraling stage until Roche lobe overflow starts at the L1
point as in the Newtonian cases. The quasiequilibrium sequences showed no
increase of the maximum energy density during the inspiraling phase as a
result of gravitational wave emission. 2) In general, equal mass corotating
BNS sequences terminate at contact point, while all the other types of
sequences (corotating nonequal mass BNSs, irrotational equal mass, and
irrotational nonequal mass binaries) end at the mass shedding point. 3)
Turning points (defining the ISCO) in the binding energy and angular momentum
of constant rest-mass sequences are found for corotating binary systems for
$\Gamma\geq 2$ and for irrotational ones for $\Gamma\geq 2.5$. It was
conjectured, that if the turning points are found for the Newtonian binaries
for some adiabatic index $\Gamma$, they should exist in the general
relativistic computations for the same value of $\Gamma$. 4) For different
mass binary systems a turning point may not always exist. 5) The deformation
of the star is determined by the orbital separation and the mass ratio and is
not affected much by its compactness. 6) The decrease of the central energy
density depends on the compactness of the star and not on that of its
companion.

The first three groups that computed irrotational BNS quasiequlibrium initial
data used the same formalism for the gravitational field and hydrodynamics but
completely different numerical techniques to solve these equations. In
particular Bonazzola et al. Bonazzola98b ; Gourgoulhon-etal-2000:2ns-initial-
data  used a multi-domain pseudo spectral method with several coordinate
systems one of which is fitted to the stellar surface. The implementation has
been done within the lorene code lorene_web  and produced solutions with
excellent accuracy with only limited amount of resources. Marronetti et al.
Marronetti99  on the other hand, used a Cartesian coordinate system and a
finite difference method. They used domain decomposition techniques which
provided a natural way of code parallelization, that reduced the processing
time enormously. For the conservation of rest-mass density equation, they
splitted the elliptic equation for the potential field into homogeneous and
inhomogeneous parts, with the homogeneous field satisfying the boundary
condition. While formally correct, this solution was difficult to implement
accurately and the results strongly depended on grid resolutions (see Sec. 2.7
for recent work on this issue). In the works by Uryū & Eriguchi Uryu00
spherical coordinates and finite differences were used while the Poisson
equations were inverted using a Green’s function approach. This was the first
application of the Komatsu-Eriguchi-Hachisu (KEH) method to BNSs in general
relativity. Separate coordinates for the gravitational and fluid equations
have been employed (surface fitted coordinates) that enabled the accurate
determination of the neutron star surface.

A systematic study of quasiequilibrium BNSs has been performed by Taniguchi &
Shibata Taniguchi2010  where a large number of systems with different mass
ratios $q=M_{1}/M_{2}$ ($0.71\leq q\leq 1$), total masses, and EOSs has been
studied under the IWM formulation using the lorene code. By constructing a
large number of BNS sequences or constant rest mass, the authors investigated
the behavior of the binding energy and total angular momentum, the endpoint of
such sequences, and the orbital angular velocity as a function of time. They
found that for piecewise polytropic EOSs the change in stellar radius at fixed
core stiffness makes the orbital angular velocity at the mass-shedding limit
vary widely, while the change in stiffness of the core EOS at fixed stellar
radius does not change the orbital angular velocity at mass-shedding
significantly. Since the less massive star in an unequal-mass binary is
tidally deformed by the companion more massive star and starts shedding mass
at larger separation than that for the equal-mass case, the orbital angular
velocity at the closest separation decreases as we decrease the mass ratio. On
the other hand the orbital angular velocity at the mass-shedding limit
increases as the neutron star mass increases, since a more massive star
becomes more compact and more difficult to be tidally disrupted for the same
EOS. This implies that BNS with massive stars need to come closer than those
with less massive stars for reaching the mass-shedding limit. The authors
derive an empirical formula for the orbital angular velocity at the mass-
shedding limit using a Newtonian argument

$M\Omega_{\rm ms}=0.27C_{1}^{3/2}\left(1+\frac{1}{q}\right)^{3/2}q^{1/2}$ (21)

where $C_{1}=M_{1}/R_{1}$ the compactness of the first neutron star.

In all the initial data above, the IWM formulation as well as helical symmetry
for the matter (Eq.(16)) is assumed. Any sequence constructed from such
initial data, satisfies Eq. (8). Stated differently, a sequence whose entropy
and rest mass is kept constant, and the flow is either corotational or
circulation conserving, approximately inspirals as a result of gravitational
wave emission. However, because such solutions represent circular orbits of
the IWM spacetime, their sequences deviate from the realistic inspiraling
orbits in full general relativity. When such initial data are used for precise
numerical relativity simulations of inspiraling BNS (typically for highly
accurate gravitational wave extraction), an approaching velocity for the
component stars is added to minimize the deviation from the inspiraling orbit
(see section 2.8).

###  2.6 Spinning binary neutron star initial data

One of the most important characteristic of a neutron star is its rotational
frequency, which in isolation has been observed to be as high as $716\;{\rm
Hz}$, corresponding to a period of $1.4\;{\rm ms}$ for PSR J1748-2446ad
Hessels2006 . For the BNS systems currently known in the Galaxy Tauris2017 ;
ZhuX2018  the rotational frequencies are typically smaller. The neutron star
in the system J1807-2500B has a period of $4.2\;{\rm ms}$ while systems
J1946+2052 Stovall2018 , J1757-1854 Cameron2018 , J0737-3039A Kramer2006  have
periods $16.96,\ 21.50,\ 22.70\;{\rm ms}$ respectively. Although the majority
of the BNS simulations are based on formulations that assume the neutron stars
to be irrotational this is no longer true when the spins are as high as
approximately 10 times the orbital period at merger ($\sim 3\;\rm ms$). In
other words for accurate gravitational wave analysis we cannot ignore neutron
star spins that are $30\;\rm ms$ or less. According to ZhuX2018  this will be
the case for binaries J1946+2052, J1757-1854, J0737-3039A which at merger will
have periods of $18.23,\ 27.09,\ 27.17\;{\rm ms}$, respectively. Note that
typical spin down rate observed for neutron stars is around
$10^{-13}-10^{-17}\mbox{s s}^{-1}$, where $10^{-15}$s s${}^{-1}=30$
ms$/10^{6}$ yr. Also, let’s not forget that event GW170817 GW170817prl  was
unable to rule out high spin priors and thus two sets of data (for low and
high spins) were consistent with the observations.

####  2.6.1 Earlier formulations and calculations

The first attempt to address the neutron star spin in a general relativistic
binary setting was by Marronetti & Shapiro Marronetti03 . Despite the problems
with their approach, it laid the foundations for the more recent advances in
this topic by Tichy Tichy11  that we will discuss in the next section. The
central points of Marronetti03  are the use of the generalized Bernoulli law

$\mathcal{L}_{\mathbf{u}}(hu_{\alpha}k^{\alpha})=0$ (22)

together with a decomposition

$v^{i}=v_{\rm\scriptscriptstyle RS}^{i}+v_{\rm\scriptscriptstyle
RI}^{i}+\epsilon^{i}_{\ jk}\Omega^{j}x^{k}$ (23)

for the fluid velocity. Here $v^{i}=u^{i}/u^{t}$ and Eq. (23) decomposes the
velocity with respect to the inertial frame $v^{i}$, into an orbital velocity
(last term) and a velocity with respect to the rotating frame (other two
terms). The latter is further written as the sum of a solenoidal
$v_{\rm\scriptscriptstyle RS}^{i}$, and an irrotational part
$v_{\rm\scriptscriptstyle RI}^{i}$, i.e.

$\partial_{i}v_{\rm\scriptscriptstyle
RS}^{i}=0,\qquad\epsilon_{kji}\partial^{j}v_{\rm\scriptscriptstyle RI}^{i}=0.$
(24)

Here $\epsilon_{kji}$ is the three-dimensional Levi-Civita tensor. Equation
(22) holds if $k^{\alpha}$ is a Killing vector, and a symmetry vector for the
stress-energy tensor. It generalizes the relativistic Bernoulli law for
stationary flows, i.e. that the enthalpy per unit rest mass is constant along
the flow lines, $hu_{t}=\rm const$, and one can show that in the cases of
corotating or irrotational flows the more stringent relation

$hu_{\alpha}k^{\alpha}={\rm const}\qquad\mbox{everywhere inside the star}$
(25)

holds. The authors make their first assumption at this point by taking Eq.
(25) to hold for the arbitrary spinning binaries too. For the components of
the fluid velocity in Eq. (23) they assume

$v_{\rm\scriptscriptstyle RS}^{i}=(a-1)\epsilon^{i}_{\
jk}\Omega^{j}(x^{k}-x_{0}^{k}),\qquad v_{\rm\scriptscriptstyle
RI}^{i}=\partial^{i}\Phi.$ (26)

The first equation above assigns a uniform angular velocity of magnitude
$(a-1)\Omega$ (a fraction of the orbital angular velocity) to the neutron star
companions in the rotating frame ($x^{k},\,x^{k}_{0}$ are the coordinate
position vector and the position of maximum baryonic density respectively),
while the second expresses the fact that every irrotational vector field is a
gradient of a scalar. The free parameter $a$ controls the approximate spin of
the neutron star $\Omega_{s}^{i}=a\Omega^{i}$ as seen by a distant observer in
the inertial frame. Since neither $a$ nor $\Omega_{s}^{i}$ have a strict
physical meaning, in order to measure the spin of each neutron star the
authors introduce the concept of circulation,

$\mathcal{C}=\oint_{c}hu_{\alpha}dx^{\alpha},$ (27)

which according to the Kelvin-Helmholtz theorem Carter1979_conservation_laws
is conserved for isentropic flows along any closed path $c$ on hypersurfaces
of constant proper time. Carter Carter1979_conservation_laws  has shown that
conservation laws like Eq. (27) can hold beyond hypersurfaces of constant
proper time, and Marronetti & Shapiro Marronetti03  assume constant
circulation along a sequence of orbits. This assumption is consistent with
neglecting the radial velocity of the fluid as the binary moves to closer
separations. In reality both the radial velocity is nonzero and the
circulation will slightly change as the stars get closer together. For an
irrotational flow the closed integral (27) over a gradient yields zero
circulation, while for a corotating binary $\mathcal{C}$ is monotonically
increasing from zero at infinite separation to a finite value (but relatively
small with respect to the maximum possible, see Sec 2.7) all the way to
merger. The irrotational part of the velocity, $v_{\rm\scriptscriptstyle
RI}^{i}$ in (26), leads to a Poisson-type of equation, similar to Eq. (18),
which the authors solved following the methods of Marronetti et al.
Marronetti99 . Sequences of different values of circulation have been
constructed and the authors reported on a spin-up effect for all cases
examined, except for the sequence with the largest value of compactness.

A different approach to construct spinning BNS initial data has been
introduced by Baumgarte & Shapiro Baumgarte:2009 ; Baumgarte:2009e . In an
effort to reduce the Euler equations to an elliptic problem they took the
divergence of Eq. (14) and derived a Poisson-type of equation for the
auxiliary variable

$H=\frac{h}{u^{t}}+V^{i}\hat{u}_{i},$ (28)

where $\hat{u}_{i}=\gamma_{i}^{\ \alpha}hu_{\alpha}$ the projected enthalpy
current. Two generalizations for the enthalpy current have been introduced:

$\displaystyle U1)\quad\hat{u}_{i}$ $\displaystyle=$ $\displaystyle
D_{i}\Phi+\eta\Omega\phi_{i},$ (29) $\displaystyle U2)\quad\hat{u}_{i}$
$\displaystyle=$ $\displaystyle D_{i}\Phi+\eta
hu^{t}(\beta_{i}+\Omega\phi_{i})$ (30)

where $\eta$ is a parameter. Decomposition ($U1$) does not lead to a
corotating flow when $\eta=1$, while ($U2$) does. For both decompositions when
$\eta=0$ one gets an irrotational flow. Although the introduction of
decompositions $(U1)$ or $(U2)$ were in the correct dirrection (see Tichy11
and Tsokaros2018 ) the main problem of the approach in Baumgarte:2009  was the
effort to reduce the Euler equation into an elliptic problem. As it was
identified by Gourgoulhon Baumgarte:2009e , by doing so, one solves only a
subset of the equilibrium equations in general, but not all of them. In
particular for the Euler equation to be satisfied, both its divergence and its
curl has to vanish. For the limiting cases of corotational or irrotational
fluid flow the curl of these equations does vanish identically, but this is
not necessarily the case for intermediate circulation.

####  2.6.2 Tichy’s formulation for spinning binary neutron star initial data

A formulation that remedies the problems of Marronetti03 ; Baumgarte:2009  has
been proposed by Tichy Tichy11 . The main ingredients of his approach was the
decomposition of the spatialy projected enthalpy current

$\hat{u}_{i}=D_{i}\Phi+s_{i},$ (31)

as in the decomposition $(U1)$ of Baumgarte:2009 . Here $s^{i}$ is the neutron
star spin, which in principle can have any form. The second main ingredient in
Tichy’s formulation is a first integral of the Euler equation as in Eq. (20).
Although this first integral has the same form as the one in the irrotational
case in reality it is different because the velocity with respect to the
corotating observer is now

$V^{i}=\frac{D^{i}\Phi+s^{i}}{hu^{t}}-\omega^{i},$ (32)

i.e. it depends on the choice of the input spin $s^{i}$. Tichy derived this
first integral by using the following assumptions: (A1) Helical symmetry for
the metric $\mathcal{L}_{\mathbf{k}}g_{\alpha\beta}=0$. (A2) Helical symmetry
for the fluid thermodynamic variables like $h$ or $\rho$. (A3) Helical
symmetry for the irrotational part of the velocity $\gamma^{\nu}_{\
i}\mathcal{L}_{\mathbf{k}}(\nabla_{\nu}\Phi)=0$. (A4) Spin is constant along
$\nabla\Phi/(hu^{t})$ which is parallel to the worldline of the star center,
$\gamma^{\nu}_{\ i}\mathcal{L}_{\nabla\Phi/(hu^{t})}(s_{\nu})=0$. (A5) Second
order terms on spin are neglected $\gamma^{\nu}_{\ i}\mathcal{L}_{{\bf
s}/(hu^{t})}(s_{i})=0$. In particular, Tichy argued that for the spinning case
one cannot assume helical symmetry on the 4-velocity
$\mathcal{L}_{\mathbf{k}}u_{\alpha}=0$, but only on the irrotational part,
assumption (A3) Tichy2016rev . In such case, from the relation
$u^{t}=g^{t\alpha}u_{\alpha}$, and assumption (A1), one has that
$\mathcal{L}_{\mathbf{k}}u^{t}\neq 0$, which implies that the conservation of
rest mass in the form of Eq. (18) (first part) will not hold. Fortunately the
assumption $\mathcal{L}_{\mathbf{k}}u_{\alpha}\neq 0$ is not necessary, and
the same set of equations can still be derived. We will come back to this
point in the following section.

Equation (31) generalizes the irrotational condition Eq. (17) with the
inclusion of a spin vector $s^{i}$. Although this vector in principle can be
chosen arbitrarily, it was shown Tichy11 ; Tichy12  that a choice

$s^{i}=\Omega_{s}^{a}\phi_{s(a)}^{i}$ (33)

minimizes differential rotation and leads to a negligible shear. Here
$\Omega_{s}^{a}$ are the three components of the spin angular velocity and
$\phi_{s(a)}^{i}$ are the rotation vectors along the neutron star’s three axes
(the subscript $(a)$ denotes a different vector). One important difference
between the works of Tichy Tichy11  and Marronetti & Shapiro Marronetti03  is
the fact that the latter decomposes $u^{i}/u^{t}=v^{i}$ in irrotational and
rotational parts while the former performs this decomposition on the spatial
enthalpy current $\hat{u}_{i}$. This choice is better first because $v^{i}$ is
not a spatial vector while $\hat{u}_{i}$ is, and second because it leads to
equations that have the correct limit in the irrotational case while the ones
in Marronetti03  do not.

With the help of Eqs. (20), (31), (32), the conservation of rest-mass density,
Eq. (18), as well as the rest of the IWM equations for the gravitational
fields, one can obtain a quasiequilibrium solution for spinning BNS following
the same steps as in the irrotational case. Spinning equilibria with their
spins aligned with the orbital angular momentum have been presented in Tichy12
using a polytropic EOS ($\Gamma=2$), and the sgrid code Tichy:2009  that
employs pseudospectral methods. Those sequences conserve the rest-mass of the
neutron star while keeping the spin parameter $\Omega_{s}$ the same. Although
a better description for BNS quasiequilibria conserve rest mass and
circulation Marronetti03 ; Tsokaros2018 , the two descriptions give similar
results Tsokaros2018 . A large suite of results using the sgrid code have been
presented in Dietrich:2015b , where different spins, eccentricities, mass
ratios, compactions, and EOSs have been explored. The authors produced highly
eccentric and eccentricity-reduced initial data (see Section 2.8), as well as
unequal-mass binaries with mass ratios $q\approx 2$. In addition they
constructed binaries with arbitrary spins, misaligned with the orbital angular
momentum, and studied precession. In all cases the dimensionless spins of the
neutron stars were $\lesssim 0.16$. Recent upgrades of the sgrid code
Tichy2019  that involved a new grid structure, the use of different
coordinates, as well as a reformulation of the equations for the conformal
factor and the velocity potential, enabled the same group to reach spins all
the way to the mass-shedding limit ($\sim 0.59$ for a $\Gamma=2$ polytrope)
and compactions up to $C\sim 0.28$.

###  2.7 Other developments

A completely different method to compute BNS initial data with arbitrary spin
has also been presented by Tsatsin and Marronetti Tsatsin2013 . Their method
does not look for solutions of the Hamiltonian and momentum constraints: their
satisfaction is only asymptotic with binary separation. Instead it consists of
a variant of metric superposition that addresses two common problems, large
stellar shape oscillations and orbital eccentricities. It reduces the former
to variations of the order of $1\%$ and offers great control over orbital
eccentricities. The authors also found that these initial data sets possess
less junk radiation than that found in standard quasiequilibria. Techniques
based on superposition and the conformal thin sandwich formulation of the
constraint equations have been used for generic initial data calculations in
East2012d , while spinning initial data based on superposition of irrotational
and corotating flows have been employed in Kastaun2013 .

Another approach to BNS initial data has been presented in Tsokaros2015 ;
Tsokaros2018  which employs the cocal library that has been used in the past
to succesfully compute a great variety of quasiequilibria, from RNS and quark
stars (axisymmetric or triaxial) Huang08 ; Uryu2016a ; Uryu2016b ; Uryu2017 ;
Zhou2017xhf ; Zhou2019hyy , binary black holes Uryu2012 ; Tsokaros2012 ;
Uryu:2012b , MRNSs with mixed poloidal and toroidal magnetic fields Uryu2014 ;
Uryu2019  (see section 4), as well as self-gravitating BHDs Tsokaros2018a
(see section 3). The main characteristics of the cocal code is the use of
finite differences and a Green’s function approach as first developed in
Uryu00  for neutron stars and in Tsokaros2007  for black holes to achieve a
convergent solution through a Picard type of iteration. The field equations
are solved in spherical coordinates in multiple patches and a smooth solution
is obtained everywhere through boundary surface integrals. For the fluid
equations surface fitted coordinates are being implemented that allow accurate
representation of the neutron star surface which is important in order to
impose boundary conditions. Comparison both with the pseudospecral code kadath
Uryu:2012b  as well as with lorene Tsokaros2015 ; Tsokaros2016  showed
excellent agreement, with the small initial differences being leveled out in
the first couple of iterations of an evolution. The cocal code can produce
accurate initial data for a variety of neutron or quark star EOSs in a
piecewise or tabulated form. It has been employed for the construction of the
highest compactness BNS system ($C=0.34$ with a total mass of
$M=7.90M_{\odot}$) in quasiequilibrium to date, using a causal and
compressible EOS Tsokaros:2019lnx .

The first spinning quasiequilibrium sequences with a nuclear EOS have been
presented by Tsokaros et al. Tsokaros2015 . The authors used the formulation
by Tichy Tichy11  but started from different assumptions. In particular using
Eqs. (31) and (32) one can rewrite the Euler Eq. (14) as

$\gamma_{i}^{\ \alpha}[\mathcal{L}_{\bf k}(hu_{\alpha})+\mathcal{L}_{\bf
V}(s_{\alpha})]+D_{i}\left(\frac{h}{u^{t}}+V^{j}D_{j}\Phi\right)=0\,.$ (34)

Instead of using assumptions (A1)-(A5) Tichy11 , they assumed (A1), (A2), and

${\rm(B1)}:\
\mathcal{L}_{\mathbf{k}}\hat{u}_{i}=0,\qquad\mbox{and}\qquad{\rm(B2)}:\
\mathcal{L}_{\bf V}s_{i}=0.$ (35)

Assumption (B1) enforces helical symmetry on the fluid velocity, consistent
with the helical symmetry assumption on the spacetime (A1), and more
importantly, consistent with the helical symmetry on the thermodynamic
variables (A2). Assumption (B1) is also necessary for taking
$\mathcal{L}_{\mathbf{k}}u^{t}=0$ and deriving the elliptic equation for the
fluid potential $\Phi$. On the other hand, assumption (B2) roughly expresses
the fact that the spin does not change significantly with respect to the
corotating observer. Although helical symmetry is increasingly more accurate
as the binary separation becomes larger (irrespective of the spin), closer
binaries with rapidly spinning neutron stars (with ms or smaller periods) will
not satisfy the helical symmetry assumptions, and in fact will not satisfy any
of (A1)-(A5), (B1),(B2).

BNS initial data within the IWM formulation have also been presented by Tacik
et al. Tacik15 ; Tacik16  based on the spells elliptic solver libraries
developed by Pfeiffer et al. Pfeiffer:2002wt ; Foucart2008 . These methods
employ pseudospectral techniques that have been applied from the same group
with great success in the binary black hole problem Cook:2004kt ;
Pfeiffer:2004qz ; Caudill:2006hw ; Boyle:2006ne . The authors implement the
spinning formulation by Tichy Tichy11  and introduce a new diagnostic for
measuring the neutron star spin in a BNS setting which is based on their work
on binary black holes Lovelace2008c . In particular they construct the
quasilocal spin angular momentum Brown1993 ; Ashtekar01a ; Ashtekar03a

$J_{\rm ql}=\frac{1}{8\pi}\oint_{S}K_{i}^{\ j}\zeta^{i}dS_{j},$ (36)

where $\zeta^{i}$ is an approximate Killing vector of the spacetime444Note
that this expression is essentially identical to the ADM angular momentum
calculated at infinity Eq. (6).. For spacetimes with axisymmetry, as for
example a rotating black hole, $\zeta^{i}$ should be chosen as the rotational
Killing vector $\phi^{i}$. For binary systems where no such symmetry exists
one needs to construct $\zeta^{i}$ through a minimization principle that
results to an eigenvalue problem Lovelace2008c . The surface of integration
$S$ is the apparent horizon for the case of a black hole, while for neutron
stars although there is no clear choice, the stellar surface is the most
natural one. The authors performed this integration both on the stellar
surface (which is calculated when a convergent solution is obtained) as well
as on spheres of increased radius outside the star. They showed that the
calculated spin is independent of the precise choice of $S$ within an accuracy
of $1\%$. For irrotational BNSs they found a quasilocal spin residue of $\sim
10^{-4}$. Highly spinning BNS initial data with $J_{\rm ql}/M_{1}^{2}\sim
0.43$, where $M_{1}$ is the ADM mass of a single star at infinite separation,
were presented, as well as precessing binaries. The spin and orbital
precession of the stars were well described by the post-Newtonian
approximation. In addition the authors implemented eccentricity control
algorithms Pfeiffer2007  for its reduction targeted for highly accurate
waveform production (see section 2.8).

Different spin measures for BNSs were examined in detail in Tsokaros et al.
Tsokaros2018 . The authors examined the relationship between the circulation
$\mathcal{C}$, dimensionless spin $J_{\rm ql}/M_{1}^{2}$, and spin parameter
$\Omega_{s}$ in Eq. (33). For the quasilocal spin, Eq. (36), the authors
assumed $\zeta^{i}$ to be the rotational vector around each star’s center
$\phi_{s}^{i}$. Realistic spinning sequences that conserve both rest mass and
circulation were presented and compared with corotating and irrotational (zero
circulation) ones. Regarding the corotating sequence, it was found that the
circulation and the dimensionless spin of the neutron stars close to the ISCO
(turning point of the $E_{b}(\Omega)$ curve) are much smaller than the
possible maximums at the mass-shedding limit. Practically this means that
corotating neutron stars at the ISCO can be considered as slowly rotating.
Similarly to isolated slowly uniformly rotating neutron stars, the circulation
and the dimensionless spin of the component stars increase linearly with the
angular velocity as the binary proceeds to closer separations.

For a spinning sequence of constant rest mass and circulation, the binding
energy (or angular momentum) as a function of angular velocity typically
follows a curve parallel to the irrotational sequence, but shifted to higher
energy (or angular momentum). This is expected since now the system contains
the spin rotational energy which remains approximately constant along the
sequence. For moderate spins (i.e. spins smaller than the corotating value at
the ISCO), $E_{b}(\Omega)$ or $J(\Omega)$ curves of a spinning sequence
exhibit an intersection with the corotating curve at some angular velocity
$\Omega_{i}$. For $\Omega<\Omega_{i}$ the binding energy (or angular momentum)
of the spinning sequence is larger than the corresponding one from the
corotating and irrotational binaries (with the same rest mass) due to the
excess of rotational energy. On the other hand for $\Omega>\Omega_{i}$
(binaries in closer separations) the binding energy (angular momentum) of the
spinning sequence although larger than the corresponding irrotational one, it
is smaller than the corotating. The dimensionless spin of each star, $J_{\rm
ql}/M_{1}^{2}$, is nearly constant for larger separations and exhibits an
increase of the order of $10-15\%$ as ones moves closer to the ISCO. For
closer separations, the choice $\zeta^{i}=\phi^{i}_{s}$ is less accurate and
one needs to compute the approximate Killing vector for the quasilocal spin
Lovelace2008c . In addition the rotational parameter $\Omega_{s}$ was found to
be approximately constant along a constant rest mass and circulation sequence.

Motivated by the circulation expression for single stars and corotating
binaries

$\mathcal{C}=\oint_{c}hu^{t}\gamma_{ij}(\beta^{i}+\Omega\phi^{i})dx^{j}\,,$
(37)

the authors proposed a modification for the fluid velocity decomposition Eq.
(31) according to

$\hat{u}_{i}=D_{i}\Phi+hu^{t}s_{i}\,,$ (38)

which apart from the shift term in Eq. (37), resembles the circulation of a
single rotating star. This decomposition is similar to the choice $(U2)$ in
Baumgarte & Shapiro Baumgarte:2009 . The velocity with respect to the
corotating observer, Eq. (32), is now modified as

$V^{i}=\frac{D^{i}\Phi}{hu^{t}}-(\omega^{i}-s^{i})\,,$ (39)

and new equations of hydrostatic equilibrium are derived. Notice that Eq. (39)
resembles the corresponding velocity of irrotational binaries Eq. (18), with
the spin vector $s^{i}$ modifying the corotating shift $\omega^{i}$. Similar
modification is present in the boundary condition, Eq. (19). Using Eq. (38)
the authors computed new sequences of constant rest mass and circulation and
found that for a fixed angular velocity, the binding energy (and angular
momentum) is now larger than the corresponding one using the decomposition of
Eq. (31). This shows that the way the 4-velocity is decomposed can influence
important quantities like the energy of the system or its angular momentum.

The irrotational and spinning formulations involve the solution of an elliptic
equation for the velocity potential $\Phi$, which typically is performed on
the so-called surface fitted coordinates for greater accuracy. These are fluid
coordinates that track the position of the surface of the star at every
iteration and are used to impose the boundary condition Eq. (19) for the
irrotational case (or a similar one for the spinning case). In order to avoid
such complication, Tsao et al. Tsao2020  developed a new technique by
employing the source term method proposed by Towers Towers2018 , where the
boundary condition is treated as a jump condition and is incorporated as
additional source terms in the Poisson equation. If the domain of the star is
denoted by $Q^{+}$, its exterior by $Q^{-}$, and its boundary by $\partial
Q^{+}$, Tsao et al. considered the boundary value problem

$\displaystyle\nabla^{2}\Phi$ $\displaystyle=S^{+}(x),\quad x\in Q^{+},$
$\displaystyle\nabla^{2}\Phi=S^{-}(x),\quad x\in Q^{-},$ (40)
$\displaystyle\left[\frac{\partial\Phi}{\partial n}\right]$
$\displaystyle=a(x),\quad x\in\partial Q^{+},$ $\displaystyle[\Phi]=b(x),\quad
x\in\partial Q^{+},$ (41)

which is a generalization of the irrotational/spinning one for the velocity
potential $\Phi$. Here we denote by $[\Phi]\equiv\Phi^{+}-\Phi^{-}=b(x)$ and
$\left[\frac{\partial\Phi}{\partial
n}\right]\equiv\left(\frac{\partial\Phi}{\partial
n}\right)^{+}-\left(\frac{\partial\Phi}{\partial n}\right)^{-}=a(x)$, with
$\Phi^{+},\left(\frac{\partial\Phi}{\partial n}\right)^{+}$ being the solution
$\Phi$ and its normal derivative evaluated at the interior of $Q^{+}$, while
$\Phi^{-},\left(\frac{\partial\Phi}{\partial n}\right)^{-}$ the same
quantities at the exterior $Q^{-}$.

The source term method converts the boundary conditions on $\partial Q^{+}$ to
jump conditions that can be absorbed into the sources. The generalized Poisson
equation that is solved in the extended domain $Q=Q^{+}\cup Q^{-}$ becomes

$\nabla^{2}\Phi=\nabla^{2}(bH)-H\nabla^{2}b-\left(a-\frac{\partial b}{\partial
n}\right)|\nabla\rho|\delta(\rho)+S,$ (42)

where $\rho$ is the rest-mass density,
$S(x)=S^{+}H(\rho(x))+S^{-}(1-H(\rho(x)))$ with $H(\rho)$ the Heaviside
function, and $\delta(\rho)$ the Dirac delta function. The authors presented a
comparison between the solution of Eq. (42) on a Cartesian grid, with the one
obtained by cocal on surface fitted coordinates, and found excellent agreement
with a maximum difference of $1.4\%$. The source term method can be used in
other problems where nonsmooth solutions are present, as for example in MRNSs.

More recently spinning BNS solutions have been presented by Papenfort et al.
Papenfort2021  using the KADATH spectral library kadath ; Grandclement09 . The
authors used Tichy’s formulation Tichy11  to produce a large suite of highly
spinning, and asymmetric systems including one with mass ratio $q=0.455$ with
the primary having dimensionless spin of $\sim 0.6$. Eccentricity reduced
techniques (as described in the next section) have also been implemented.

###  2.8 Eccentricity in binary neutron stars

All BNS initial data computations that we have mentioned are based on the IWM
formulation, whose equations hold even without the existence of a helical
Killing vector. Since gravitational radiation reaction and the accompanying
approaching velocity are neglected, quasicircular binary initial data of this
kind exhibit a small but nonzero deviation from strict circularity of the
order of $\sim 0.01$ Despite its small value, the residual eccentricity
becomes problematic in the evolution of a binary when accurate gravitational
wave analysis is needed in order to determine various parameters such as the
tidal deformability Favata2014 ; Yagi2014c ; Wade2014 , which is crucial for
constraining the neutron star EOS.

Kyutoku et al. Kyutoku2014  presented a method to further reduce the orbital
eccentricity by using a similar methodology employed for binary black holes
Pfeiffer2007 . The main idea is to incorporate an approaching velocity term in
the formulation that can be adjusted through sequential evolutions until the
eccentricity is reduced. The difference between binary black holes and BNSs is
that this approaching velocity is incorporated on the apparent horizon
boundary conditions in the black hole case, while for neutron stars it is
through the helical Killing vector that controls the hydrodynamic equations.
The procedure starts with the computation of quasicircular initial data and a
modified symmetry vector for the hydrodynamical fields as

$k^{\alpha}=t^{\alpha}+\Omega\phi^{\alpha}+v\frac{r}{r_{0}}(\partial_{r})^{\alpha}\,,$
(43)

instead of Eq. (12). Here $v$ is the radial velocity, and $r_{0}$ the
separation from the coordinate origin. Note that the extra approaching
velocity term is a conformal Killing vector of the flat metric, and hence does
not affect the gravitational field equations. The quasicircular initial data
($v=0$) are then evolved for a sufficient time interval whose duration on one
hand has to be long enough to include more than one modulation eccentricity
cycles, and on the other, short enough to avoid strong influence of a long-
term secular evolution. The authors assumed an interval of $[0.5P,3P]$, where
$P=2\pi/\Omega$ is the initial orbital period. To measure the eccentricity
they locate the star center $x^{i}_{{}_{\rm NS}}=(x_{{}_{\rm NS}},y_{{}_{\rm
NS}},0)$ through the maximum of the conserved rest-mass density, $\rho\alpha
u^{t}\sqrt{\gamma}$, and compute the coordinate orbital separation and orbital
phase as

$\displaystyle d(t)$ $\displaystyle=$ $\displaystyle\sqrt{(x_{{}_{\rm
NS,1}}-x_{{}_{\rm NS,2}})^{2}+(y_{{}_{\rm NS,1}}-y_{{}_{\rm NS,2}})^{2}},$
(44) $\displaystyle\phi(t)$ $\displaystyle=$
$\displaystyle\tan^{-1}\left(\frac{y_{{}_{\rm NS,1}}-y_{{}_{\rm
NS,2}}}{x_{{}_{\rm NS,1}}-x_{{}_{\rm NS,2}}}\right)+2\pi N,$ (45)

where $N$ is the number of orbits (in order to make $\phi$ continuous) and
indices $1,2$ refer to the two neutron stars. Following Pfeiffer2007 ;
Boyle:2007ft , Kyutoku et al. assume

$\frac{d\Omega}{dt}=A_{0}+A_{1}t+B\cos(\omega t+\phi_{0})\,,$ (46)

where $\\{A_{0},A_{1},B,\omega,\phi_{0}\\}$ are fitting parameters. The term
$A_{0}+A_{1}t$ is due to radiation reaction, while the term $B\cos(\omega
t+\phi_{0})$ represents the modulation from the nonzero eccentricity $e$. For
a perfectly circular orbit it should be $B=0$. For a Newtonian orbit with
small eccentricity ($e<<1$), it is $\Omega(t)\approx\Omega[1+2e\sin(\omega
t+\phi_{0})]$ which leads to the estimate $e\approx\frac{|B|}{2\omega\Omega}$.
Following Newtonian considerations Buonanno2011  the adjustments in the
orbital velocity and the approaching velocity turn out to be

$\Omega\rightarrow\Omega-\frac{B\omega\sin\phi_{0}}{4\Omega^{2}},\qquad\mbox{and}\qquad
v\rightarrow v+\frac{Bd\cos\phi_{0}}{4\Omega}.$ (47)

Using the updated values of $\Omega$ and $v$, the elliptic initial value
equations are then solved again, and the whole procedure is repeated until an
acceptable value of eccentricity ($e\lesssim 10^{-3}$) is obtained. Typically
this requires 3 such iterations. By calculating the gravitational waves in
every iteration the authors prove that this eccentricity reduction procedure
leads to gauge invariant results, i.e. the considerations based on the
coordinate orbital separation are not gauge artifacts.

The authors provide an alternative way to measure the eccentricity based on
the gravitational wave angular velocity $\Omega_{\rm gw}(t)$,

$e_{\rm gw}(t)=\frac{\Omega_{\rm gw}(t)-\Omega_{\rm gw,fit}(t)}{2\Omega_{\rm
gw,fit}(t)}$ (48)

where $\Omega_{\rm gw,fit}(t)$ a fourth order polynomial (5 fitting
constants). They showed that the results obtained with this method are similar
to the ones based on the orbital motion adding confidence about the
reliability of both methods.

Similar algorithms to remove the eccentricity were presented by Tacik et al.
Tacik15 ; Tacik16 , Dietrich et al. Dietrich:2015b , and Papenfort et al.
Papenfort2021 . Dietrich et al. Dietrich:2015b  compared two merger
simulations with the SLy EOS one from quasicircular initial data having
$e=1.241\times 10^{-2}$, and the other from eccentricity reduced initial data
having $e=8.7\times 10^{-4}$. They found that the phase difference
$\delta\phi_{22}$ oscillates in the range of $[-0.06,0.06]$ rad, i.e. it
produces an approximate dephasing of $\sim 0.12$ rad. The amplitude of the
non-eccentricity-reduced data oscillates around the eccentricity-reduced ones
by $5\%$ at early times and decrease as the system approaches merger. Note
that in this comparison the initial data are at a large distance of $\sim 10$
orbits before merger. At much closer distances of $\sim 3$ orbits, Tsokaros et
al. Tsokaros2016  found that the dephasing that comes by evolving
quasicircular initial data from the cocal and lorene initial value codes and
evolved by the WhiskyTHC Radice2013b ; Radice2013c  evolution code, can be 10
times larger. Therefore for accurate gravitational wave analysis one has to
take multiple factors into consideration including the truncation errors at
various resolutions.

When BNSs merge they will follow highly circular trajectories, since
gravitational radiation reaction circularizes the orbit Peters:1964 . Despite
of that, BNS mergers in eccentric orbits are still possible, either by
dynamical interactions in dense stellar regions, such as globular clusters
Oleary2009 ; Lee2010 ; Tsang2013 , or by exciting their eccentricity by, e.g.,
the Kozai mechanism in a hierarchical triple Antonini2012 ; Naoz2012 ;
Seto2013 ; Antognini2014 . Although eccentric BNSs have been constructed by
Gold et al. Gold2012 , as well by East & Pretorius East2012c  using different
approximations, the first consistent method to construct such initial data has
been presented by Moldenhauer et al. Moldenhauer2014 . Their method
generalizes the approximate helical Killing vector approach that is used to
solve the Euler equation in quasicircular binaries, to a pair of inscribed
helical symmetry vectors (one for each star), and allows them to provide
initial data for binary neutron stars with arbitrary eccentricity. The authors
assume each star center moves along a segment of an elliptic orbit at apoapsis
(which they take to be on the x-axis), and approximate this small orbital
segment by the circle inscribed into the elliptical orbit there. The radii of
the inscribed circles are $r_{c1,2}=(1-e)d_{1,2}$, where $e$ is the
eccentricity of the elliptical orbit, and $d_{1,2}=|x_{1,2}-x_{\rm cm}|$ are
the distances of the neutron stars from the center of mass. The centers of the
inscribed circles are at $x_{c1,2}=x_{\rm cm}+e(x_{1,2}-x_{\rm cm})$, thus the
approximate Killing vector (near each star) for the elliptical orbit is

$k_{\rm
ecc1,2}^{\alpha}=t^{\alpha}+\Omega[(x-x_{c1,2})y^{\alpha}-yx^{\alpha}]\,.$
(49)

In addition in order to accommodate for the energy loss due to gravitational
wave emission, a radial velocity can be added to Eq. (49) similar in spirit to
Eq. (43) Kyutoku2014 . Consistent initial data with eccentricities as large as
$e=0.5$ have been presented in Dietrich:2015b .

###  2.9 Non-conformally flat binary neutron star initial data

In an effort to correct the error coming from the choice of the conformally
flat three-geometry, two groups presented an improved formulation for initial
data starting from different viewpoints. Bonazzola et al. Bonazzola:2003dm
aimed to reformulate the whole 3+1 numerical relativity system (the Einstein
system), and instead of using a free-evolution scheme (like the Baumgarte-
Shapiro-Shibata-Nakamura formulation Shibata95 ; Baumgarte99 ) that involves
hyperbolic equations, they proposed to use a fully constrained evolution
method. In that formulation one only solves the maximum number of elliptic
equations and the minimum number of hyperbolic equations: the two wave
equations corresponding to the two degrees of freedom of the gravitational
field. To achieve this they assumed a decomposition of the conformal spatial
metric as

$\tilde{\gamma}_{ij}=f_{ij}+h_{ij},$ (50)

where $f_{ij}$ is the flat metric in the chart of the 3-dim hypersurface, and
$h_{ij}$ the components that need to be evaluated. The authors imposed a
condition for the determinants, $\det\tilde{\gamma}_{ij}=\det f_{ij}$, and
used maximal slicing and a generalization of the Dirac gauge to curvilinear
coordinates

${\raise 4.30554pt\hbox{${}^{\
\circ}$}}\\!\\!\\!\\!\\!D_{b}\tilde{\gamma}^{ab}=0.$ (51)

This gauge fixes the spatial coordinates $x^{i}$ in each hypersurface
$\Sigma_{t}$, and has been introduced by Dirac Dirac59a  as a way to fix the
coordinates in the Hamiltonian formulation of general relativity. The non-
conformal flat part of the metric $h_{ij}$ satisfies a wave-like equation
which in spherical coordinates and in the presence of the Dirac gauge can be
reduced to two scalar wave equations. This new evolution scheme in effect
determines a new initial data formulation also, where the quantities
$\\{\psi,\alpha,\beta^{i},h_{ij}\\}$ must be determined from a set of elliptic
equations.

At the same time Shibata et al. Shibata04a  proposed a formulation for BNSs
where the full Einstein system is solved similar to Bonazzola:2003dm . They
also used the Dirac gauge to derive elliptic equations for the non-conformal
part of the methric $h_{ij}$, but not in spherical coordinates. The authors
provide asymptotic falloff conditions for the lapse, the shift, the spatial
metric, and the extrinsic curvature which ensure the equality of ADM and Komar
masses, $M_{\scalebox{0.6}{K}}=M$, therefore generalizing the results by Beig
Beig1978 , and Ashtekar & Magnon-Ashtekar Ashtekar79a  for stationary systems.
The equality of the ADM and Komar masses is related to a virial relation for
equilibrium,

$\int_{\Sigma}x^{a}\gamma_{a}^{\ \alpha}\nabla_{\beta}T_{\alpha}^{\
\beta}\sqrt{-g}d^{3}x=0,$ (52)

and the first law, Eq. (7), that are used for evaluating the accuracy of the
numerical solutions including non-axisymmetric ones computed from the above
formulation.

The first quasiequilibrium sequences of irrotational BNS under the formulation
described above were calculated by Uryū et al. Uryu2006 . Together with the
assumptions of maximal slicing and the Dirac gauge (or spatially transverse
condition), the authors restrict the time derivative terms so that all
components of the field equations are elliptic, and hence that all metric
components, including the spatial metric have Coulomb-type falloff. In
particular for the conformal metric a condition
$\partial_{t}\tilde{\gamma}_{ij}=0$ is imposed while for the extrinsic
curvature and the fluid variables helical symmetry is assumed
$\mathcal{L}_{\bf k}K_{ij}=0$.

In order to impose conditions (51) and have a self-consistent iteration
scheme, an adjustment is necessary for the $h_{ij}$. This is accomplished by
introducing new gauge vector potentials $\xi^{a}$ as in Uryu:2009ye , (or
Uryu2016a  Eq. (29)-(32)) through the transformation

$\delta\gamma^{ab}\rightarrow\delta\gamma^{ab}\,-\,{\raise
4.30554pt\hbox{${}^{\ \circ}$}}\\!\\!\\!\\!\\!D^{a}\xi^{b}\,-\,{\raise
4.30554pt\hbox{${}^{\ \circ}$}}\\!\\!\\!\\!\\!D^{b}\xi^{a},$ (53)

which when combined with Eq. (51), yield another set of elliptic equations for
$\xi^{a}$. The augmented system for
$\\{\psi,\alpha\psi,\beta^{a},h_{ij},\xi^{a}\\}$ is then solved using a self-
consistent method. The authors named this new formulation as “waveless”
formulation due to the absence of gravitational waves in the constructed
spacetimes. The work of Uryū et al. Uryu2006  has two main conclusions. First
when one computes the binding energy of the system it is found that although
at large separations the results agree with the ones coming from post-
Newtonian or from using the IWM formulation, this is not true any more for
close orbits. In particular it was shown that the IWM formalism overestimates
$|E_{b}|$ by neglecting the nonconformally flat part of the conformal
geometry. The second important conclusion of Uryu2006  is that the IWM
formulation does not enforce circularity as accurately as the waveless
formulation even for large separations. In order to prove that, the authors
calculated the norm of $\mathcal{L}_{\bf k}K_{ij}$ which should be zero in
exact helical symmetry. They found that in the waveless formulation the number
is at least one order of magnitude less than in the IWM formulation and the
discrepancy is larger at larger separations.

These results have been extended in the sequel work of Uryū et al. Uryu:2009ye
where two formulations for nonconformally flat initial data are examined:
waveless and near-zone helically symmetric Yoshida06a . In each formulation,
the Einstein-Euler system, written in 3+1 form on an asymptotically flat
spacelike hypersurface, is exactly solved for all metric components, including
the spatially nonconformally flat potentials, and assuming an irrotational
flow. In the helically symmetric approach helical symmetry is imposed in the
near zone, from the center of mass to a radius $\sim\lambda=\pi/\Omega$, and
either the waveless formulation is applied outside, or the computational
domain is truncated at that radius. Here $\lambda$ is the wavelength of the
dominant mode (primarily $\ell=m=2$ quadrupole) of the gravitational waves
expected to be radiated from the system. As the compactness of the component
neutron stars increases, the behavior of the binding energy and angular
momentum of binary sequences more closely approximates that of point masses.
By using a variety of realistic EOSs in a piecewise representation the authors
show that this is true for the IWM and waveless/helically symmetric sequences,
but the effect is more pronounced for the IWM models than in the
waveless/helically symmetric ones. The behavior of the IWM sequence was
interpreted as the effacing of the tidal effects due to the spatially
conformally flat approximation. In effect this correction reflects an
overestimation in the IWM formulation by 1 cycle in the gravitational wave
phase during the last several orbits.

##  3 Self-gravitating black hole-disks.

BHD systems are omnipresent in astrophysics, from the core collapse of massive
stars Woosley1993 ; MacFadyen1999 , and the cores of active galactic nuclei
LyndenBell1969 ; Shakura1973 ; Paczynski1978 , to the merger of two compact
objects at least one of which is not a black hole Shibata99d ; Shibata06d ;
Etienne2007b ; Liu:2008xy ; Baiotti08 ; Anderson2007 ; Rezzolla:2010 ;
Lovelace2013 . In addition, stellar-mass binaries in active galactic nuclei
and quasars or massive black hole binaries in extended disks are also possible
scenarios with high astrophysical interest. Such systems constitute excellent
candidates of “multimessenger astronomy” since they will produce copious
amounts of electromagnetic and gravitational radiation, detectable by the
future Laser Interferometer Space Antenna (LISA) Amaro2017 .

In many of the cases above the spacetime can be approximated by the Kerr
solution, and one assumes that the fluid motion is not backreacting on the
gravitational sector. However it is possible that a massive black hole (or a
binary black hole system) is immersed in an exteded disk with mass comparable
or even greater than the black hole itself. In other words, there may be a
timeframe where such BHD spacetimes cannot be described by the Kerr metric and
the self-gravity of the disk needs to be taken into account. The disk
geometry, and density profile can instigate a number of instabilities
Papaloizou84 ; Papaloizou85 ; Blaes1985 ; Goldreich1986 ; Narayan1987 ;
Abramowicz1980  that can affect the gravitational wave signal detected by LISA
Korobkin2011 ; Korobkin2013 ; Mewes2016 ; Kiuchi2011b ; Wessel2020 .
Sufficiently massive disks can cause spin precession or even spin flips
Gergely:07 . At the same time gas accretion can appreciably change the mass
and the spin of the black hole(s) Gammie:2003qi . If a binary black hole is
found within a massive disk its gravitational pull, hydrodynamical drag, and
planetary migration can change the trajectories of the companions leaving an
imprint on the inspiral rate and the GW phase evolution Petrich89 ;
Barausse2007 ; Barausse2007a ; Yunes2011 . More recently, the announcement of
the high-mass binary black hole, GW190521 by aLIGO and VIRGO Abbott2020 ,
initiated a debate regarding the several possibilities for forming black holes
in the mass gap $\sim 50-120\ M_{\odot}$. One of the intriguing proposals
Shibata2021  was that the source of GW190521 is a stellar collapse of a very
massive star leading temporarily to a black hole with a massive disk that is
dynamically unstable to the one-armed spiral-shape deformation.

For all these reasons simulating self-gravitating BHDs is important in order
to understand the plethora of current and future observations. Only by
including self-gravity in full general relativity and tracking the
nonaxisymmetric perturbations that it may trigger can gravitational waves from
the disk be calculated reliably. The methods that we will describe below aim
towards calculating self-gravitating BHD (quasi)equilibria, that can serve as
self-consistent initial data for various evolutionary scenarios.

###  3.1 Formulation

Most of the work in BHDs as well as in MRNSs has been done under the
assumption of a circular, stationary, and axisymmetric spacetime. In these
spacetimes (used extensively in the construction of RNSs Friedman2012 ;
Paschalidis2017b ), the line element, Eq. (2) can be constrained to a smaller
number of unknown metric components. In particular a spacetime is stationary
and axisymmetric if there exist a timelike Killing vector $\xi^{\alpha}$ with
integral curves $b(t),\,t\in\mathbb{R}$, and a spacelike Killing vector
$\chi^{\alpha}$ with closed integral curves $c(\phi),\,\phi\in[0,2\pi]$. Here
$t,\,\phi$ are arbitrary parameters. Carter has shown Carter1970  that there
is no loss of generality if one further assumes that these two vector fields
commute,
$\xi^{\alpha}\nabla_{\alpha}\chi^{\beta}=\chi^{\alpha}\nabla_{\alpha}\xi^{\beta}$.
The commutation property can be used to promote parameters $t,\,\phi$ into
coordinates of the spacetime manifold. In particular one can choose
coordinates $\\{t,x^{1},x^{2},\phi\\}$ such that
$\xi^{\alpha}:=(\partial_{t})^{\alpha}$ and
$\chi^{\alpha}:=(\partial_{\phi})^{\alpha}$. In that case the 10 spacetime
metric components $g_{\mu\nu}$ will depend only on coordinates $x^{1},x^{2}$
and not on $t,\,\phi$.

If in addition the spacetime is circular, then Papapetrou66 ; Kundt66 ; Wald84

$\xi^{\mu}R_{\mu}^{\
[\alpha}\xi^{\beta}\chi^{\gamma]}=0,\qquad\chi^{\mu}R_{\mu}^{\
[\alpha}\xi^{\beta}\chi^{\gamma]}=0,$ (54)

where $R_{\alpha\beta}$ is the Ricci tensor, and square brackets denote
antisymmetrization. Eq. (54) guarantee that the 2-dimensional planes
orthogonal to $\xi^{\alpha}$ and $\chi^{\alpha}$ are integrable (we have
tacitly assumed that the spacetime is asymptotically flat, therefore a
rotational axis where $\chi^{\alpha}$ vanishes, must exist). Under such
assumptions the spacetime metric can be further simplified as

$ds^{2}=g_{tt}dt^{2}+2g_{t\phi}dtd\phi+g_{\phi\phi}d\phi^{2}+g_{\scriptscriptstyle
AB}dx^{\scriptscriptstyle A}dx^{\scriptscriptstyle B},$ (55)

where $A,B\in\\{x^{1},x^{2}\\}$, with only 6 nonzero components. A proper
choice of $x^{1},x^{2}$ will make the 2-dimensional metric
$g_{\scriptscriptstyle AB}$ diagonal (for example spherical coordinates
$\\{x^{1}=r,x^{2}=\theta\\}$) and since all 2-dimensional metrics are
conformally related, one more component can be eliminated. At the end, the
form of a circular, stationary and axisymmetric spacetime in the so-called
quasi-isotropic form is Bardeen71 ; Shibata2007

$ds^{2}=-\alpha^{2}dt^{2}+\psi^{4}[e^{2q}(dr^{2}+r^{2}d\theta^{2})+r^{2}\sin^{2}\theta(d\phi+\beta
dt)^{2}]$ (56)

where all four functions $\alpha,\psi,q,\beta$ depend on $r,\,\theta$ only.
From now on, in accordance with the coordinate notation of the metric, we will
denote the Killing vectors $\xi^{\alpha}$ and $\chi^{\alpha}$ as $t^{\alpha}$
and $\phi^{\alpha}$.

In the presence of nonzero sources in the Einstein’s equations, the Ricci
tensor in conditions (54) can be replaced by the stress-energy tensor, which
result into imposing the circularity property on them. In particular for a
perfect fluid, conditions (54) imply that

$u^{[\alpha}t^{\beta}\phi^{\gamma]}=0,$ (57)

i.e. the fluid 4-velocity belongs to the hyperplane spanned by $t^{\alpha}$
and $\phi^{\alpha}$, or

$u^{\alpha}=u^{t}(t^{\alpha}+\Omega\phi^{\alpha}),\qquad\mbox{with}\qquad\Omega:=\frac{u^{\phi}}{u^{t}}=\frac{d\phi}{dt}.$
(58)

Such kind of fluid flow, which is typically assumed in RNSs, will be employed
for the BHDs considered below. In this case, the rest mass conservation
equation, $\nabla_{\alpha}(\rho u^{\alpha})=0$, is identically satisfied, and
the equations of motion for the fluid
$(g_{\alpha\beta}+u_{\alpha}u_{\beta})\nabla_{\gamma}T^{\beta\gamma}=0$ yield

$\frac{\nabla_{i}p}{\epsilon+p}=-\nabla_{i}\ln(u_{t})+\frac{\Omega\nabla_{i}\ell}{1-\Omega\ell},$
(59)

where

$\displaystyle\displaystyle\ell:=-\frac{u_{\phi}}{u_{t}}=-\frac{g_{t\phi}+\Omega
g_{\phi\phi}}{g_{tt}+\Omega
g_{t\phi}}\quad,\qquad\Omega=\frac{u^{\phi}}{u^{t}}=-\frac{g_{t\phi}+\ell
g_{tt}}{g_{\phi\phi}+\ell g_{t\phi}},$ (60) $\displaystyle\displaystyle
u_{t}=\sqrt{\frac{g_{t\phi}^{2}-g_{tt}g_{\phi\phi}}{g_{\phi\phi}+2\ell
g_{t\phi}+\ell^{2}g_{tt}}}.$ (61)

Eq.(59) is integrable if a one-parameter EOS555We often use a barotropic EOS,
$p=p(\epsilon)$, or $p=p(\rho)$ and $\epsilon=\epsilon(\rho)$. and a rotation
law, $\ell=\ell(\Omega)$ (or $\Omega=\Omega(\ell)$), are prescribed. The same
equation can be written alternatively as

$\nabla_{\alpha}\ln\frac{h}{u^{t}}+u^{t}u_{\phi}\nabla_{\alpha}\Omega=0,$ (62)

where $h$ is the relativistic specific enthalpy. Similarly Eq.(62) is
integrable if a one-parameter EOS and the rotation law
$j:=u^{t}u_{\phi}=j(\Omega)$ (or $\Omega=\Omega(j)$) are prescribed. The first
integral of Eq.(59) or Eq. (62) defines the surfaces of constant pressure
(which correspond to the equipotential surfaces in Newtonian theory). When the
self-gravity of the disk is negligible, these surfaces are derived from the
Kerr metric components and the rotation law $\ell=\ell(\Omega)$. Such
equilibria, the so-called Polish doughnuts, have been computed by Fishbone &
Moncrief Fishbone76 , Abramowicz et al. Abramowicz78 , and Kozlowski et al.
Kozlowski1978 . They constitute the workhorse of most BHD studies where disk’s
gravity is neglected. They are also taken as the starting point for the
calculation of a self-gravitating disk around a black hole.

For the case where $\ell=\mbox{constant}$, and assuming a polytropic EOS
$p=k\rho^{\Gamma}$, Eq.(59) can be integrated as

$\rho=\left[\left(\frac{\Gamma-1}{\Gamma k}\right)\left(\frac{u_{\rm
in}-u_{t}}{u_{t}}\right)\right]^{\frac{1}{\Gamma-1}}.$ (63)

Here $u_{\rm in}=u_{t}(r_{\rm in},\theta=\pi/2)$ is the specific energy at the
inner point of the disk on the equatorial plane. For a radiation dominated
disk, $\Gamma=4/3$. Ignoring the self-gravity of the disk and assuming a Kerr
metric results in the determination of $u_{t}$ through Eq. (61), and therefore
the rest-mass density $\rho$. On the other hand if one assumes a general
differential rotation law of the form $\ell=\ell(\Omega)$, either of Eqs. (60)
is solved pointwise to calculate $\ell=\ell(x,y,z)$ (or
$\Omega=\Omega(x,y,z)$), which in turn is used to compute $u_{t}$ from Eq.
(61), and finally the rest-mass density is obtained from the first integral of
Eq. (59) or (62).

The inclusion of disk’s self-gravity is achieved by solving simultaneously the
Einstein equations in conjunction with Eq. (59) or Eq. (62) through an
iterative method. An assumption about the EOS, and the disk’s differential
rotation law as mentioned above, is required. Thus the procedure described in
the paragraph above for a non self-gravitating disk, is now repeated at every
iteration step, with the difference being the metric components that determine
$u_{t}$ are now the solutions of the Einstein system, and not the ones coming
from the Kerr spacetime.

###  3.2 Mass, angular momentum, and Smarr formula

In the case of axisymmetric systems the angular momentum can be expressed
covariantly (and therefore in a gauge invariant way) as a Komar integral

$J_{\scalebox{0.6}{K}}=\frac{1}{8\pi}\oint_{S_{{}_{\infty}}}\nabla^{\alpha}\phi^{\beta}dS_{\alpha\beta}=\underbrace{-\frac{1}{8\pi}\oint_{H}\nabla^{\alpha}\phi^{\beta}dS_{\alpha\beta}}_{J_{h}}+\underbrace{\int_{\Sigma_{t}}T^{\alpha}_{\
\beta}\phi^{\beta}dS_{\alpha}}_{J_{t}},$ (64)

where $S_{{}_{\infty}}$ denotes a sphere whose radius tends at infinity. From
Eq. (6) one can show that the Komar angular momentum coincides with the ADM
angular momentum, $J=J_{\scalebox{0.6}{K}}$. Using Stokes theorem, the surface
integral at infinity can be converted to a surface integral on the black hole
horizon $H$, and a volume integral on the spatial slice $\Sigma_{t}$. The
former can be identified with the black hole angular momentum $J_{h}$, while
the latter with the disk angular momentum $J_{t}$ Bardeen73b ;
Bardeen73_Stars_Disks_Black_Holes ; Carter73_BH_Equilibrium_States , which
with the help of
$\nabla_{\beta}\nabla_{\alpha}\phi^{\beta}=R_{\alpha\beta}\phi^{\beta}$ can be
expressed in terms of the stress-energy tensor $T_{\alpha\beta}$, Eq. (64).

In an analogous manner to the Komar angular momentum one can define in the
presence of the timelike Killing vector field $t^{\alpha}$ (stationary
spacetimes) the Komar mass as a surface integral at infinity

$M_{\scalebox{0.6}{K}}=-\frac{1}{4\pi}\oint_{S_{{}_{\infty}}}\nabla^{\alpha}t^{\beta}dS_{\alpha\beta}=\underbrace{\frac{1}{4\pi}\oint_{H}\nabla^{\alpha}t^{\beta}dS_{\alpha\beta}}_{M_{h}}+\underbrace{\int_{\Sigma_{t}}\left(T\delta^{\alpha}_{\
\beta}-2T^{\alpha}_{\ \beta}\right)t^{\beta}dS_{\alpha}}_{M_{t}},$ (65)

where $T=T^{\mu}_{\ \mu}$. Beig Beig1978 , Ashtekar & Magnon-Ashtekar
Ashtekar79a , and Shibata, et al.  Shibata04a  have proved that the Komar and
ADM masses are identical, $M=M_{\scalebox{0.6}{K}}$, for stationary systems.
The integral over the horizon can be identified as the black hole mass
$M_{h}$, while the volume integral as the gravitational mass of the disk
$M_{t}$. In the vacuum case the volume integral over the stress-energy tensor
$M_{t}=0$, and $M_{h}$ can be evaluated through the Smarr formula Smarr73a

$M_{h}=\frac{1}{4\pi}\kappa_{h}A_{h}+2\omega_{h}J_{h}\,.$ (66)

Here $\kappa_{h}$ is the surface gravity of the black hole, $A_{h}$ the area
of the horizon, and $\omega_{h}$ the frame-dragging at the horizon defined by
$\omega_{h}=-g_{t\phi}/g_{\phi\phi}$. For the Kerr black hole,
$\kappa_{h}=\sqrt{M^{2}-a^{2}}/(2Mr_{+})$, and $A_{h}=8\pi Mr_{+}$, with
$r_{+}=M+\sqrt{M^{2}-a^{2}}$ the radial coordinate of the event horizon in
Boyer-Lindquist coordinates. One important point already discussed in
Bardeen73_Stars_Disks_Black_Holes  is that in the presence of matter around
the black hole, $M_{h}$ is not a good choice for the gravitational mass of the
black hole and can lead to erroneous results (for example that
$J_{h}/M_{h}^{2}>1$). On the contrary, $J_{h}$ is always a good measure of the
black hole angular momentum.

###  3.3 Self-gravitating thin disk around black hole

The first general relativistic computation of a BHD system in which the self-
gravity of the disk was taken into account has been computed by Lanza
Lanza1992a  for a thin pressureless (dust) disk whose stress energy tensor is
$T_{\alpha\beta}=\sigma u_{\alpha}u_{\beta}$, $\sigma$ being its surface
energy density. In his method $\sigma$ was prescribed, and then the metric and
the rotation law were iterated until the field equations and the equation of
motion were satisfied. For computing the metric, the author used the Bardeen &
Wagoner formulation Bardeen71  in which the self-gravity of a massive
infinitesimally thin disk was incorporated as boundary conditions of the
Einstein’s equations at the equatorial plane. Although Lanza’s black hole-disk
model involves a significant simplification in computing the thin disk, its
qualitative behavior agrees well with that of the self-gravitating equilibrium
thick disk (having nonzero pressure) around a black hole.

The starting point of Lanza1992a  was the determination of a non self-
gravitating density profile $\rho$, Eq. (63), around a rotating black hole
Abramowicz78 ; Kozlowski1978  under the assumption of an $\ell={\rm const}$
law and the polytropic EOS for a radiation dominated gas ($\Gamma=4/3$). The
free quantities to be chosen are the specific angular momentum $\ell$, and the
inner point of the disk $r_{\rm in}$. During the iteration procedure to
incorporate the disk’s self-gravity, the author fixes the surface density
which is calculated as a quadrature of $\rho$ in the $\theta$ direction,
$\sigma(r)=\int_{0}^{\pi}\rho r\psi^{2}e^{q}d\theta$, taken from the above
non-self-gravitating thick disk. The projection of the divergence of the
stress-energy tensor orthogonal to the 4-velocity (Euler equation), when
$p=0$, leads to the geodesic equation

$\partial_{\alpha}g_{tt}+2\Omega\partial_{\alpha}g_{t\phi}+\Omega^{2}\partial_{\alpha}g_{\phi\phi}=0\,,$
(67)

which the author employs to calculate $\Omega$ on the equatorial plane. The
specific angular momentum profile, and the 4-velocity are then calculated from
Eqs. (60), (61), while the rest mass of the disk, $M_{0}$, is computed as a
volume integral over the fluid support.

Since the reference model of non-self-gravitating thick disk is composed of an
ideal gas plus radiation, the polytropic constant $k$ in Eq. (63) is a
function of the ratio $\delta=p_{\rm gas}/p$, $p$ being the total pressure. By
varying $\delta$, and thus $k$, solutions of BHDs are computed where it was
found that for fixed specific angular momentum $\ell$, the disk rest mass
$M_{0}$ increases with $\delta$, as the central density becomes larger. Also
for fixed ratio $\delta$, increasing $\ell$ results to larger disks and
therefore $M_{0}$ also increases.

Lanza computes BHD sequences of increasing rest mass, where the disk has the
same inner radius $r_{\rm in}$ and specific angular momentum $\ell$, while the
black hole has fixed angular momentum $J_{h}$ (aligned with the disk angular
momentum), and area $A_{h}$. He found that along this sequence:

  1. 1.

The ADM mass and angular momentum of the system are increasing.

  2. 2.

The mass of the black hole, given by the diagnostic $M_{h}$, is decreasing.

  3. 3.

The black hole surface gravity $\kappa_{h}$ is decreasing.

  4. 4.

The black hole horizon radius $r_{h}$ is decreasing.

  5. 5.

The black hole horizon angular velocity $\Omega_{h}$ is decreasing.

Finding (1) is not surprising; it is due to the increase of the gravitational
energy and angular momentum caused by the disk’s self-gravity. Finding (2) is
also expected since $M_{h}$ contains part of the binding energy of the system
Bardeen73_Stars_Disks_Black_Holes , that becomes increasingly negative as the
mass of the disk increases. In fact $M_{h}$ is not a good diagnostic for the
black hole mass, as clearly shown by Shibata Shibata2007 .

The fact that the surface gravity of the black hole decreases (finding (3))
can be explained as follows: As the gravitational field of the disk becomes
significant, a zero angular momentum observer located between the black hole
and the disk will feel the outward pull of the disk decreasing his physical
acceleration with respect to infinity, i.e. decreasing the effective gravity.
In the limit as one goes to the horizon this is decribed by $\kappa_{h}$666A
particle between the black hole and the disk will need less angular momentum
to stay in a Keplerian orbit Abramowicz1984 ; Chakrabarti1988 .. Because the
surface gravity and the radius of the horizon are related, $\kappa_{h}=4\pi
r_{h}/A_{h}$, and since $A_{h}$ is kept constant, the radius of the horizon
behaves analogously (finding (4)).

Will Will1974 ; Will1975  has pointed that in the presence of self-gravitating
matter one cannot make a clear distiction for the individual contributions of
the black hole and the disk. This makes possible for the black hole to have
zero angular momentum but non zero angular velocity or zero angular velocity
and negative angular momentum. This effect is due to the dragging of the
inertial frames by the external self-gravitating disk. In the case of a slowly
rotating black hole the metric function $-\beta$ has no longer a maximum at
the horizon but at the center of the disk and this maximum increases with
$M_{0}$. This explains finding (5).

Lanza also computes sequences of increasing disk mass with fixed $\Omega_{h}$,
and horizon radius $r_{h}$, for a slowly rotating black hole (small and
positive $\Omega_{h}$). Now the black hole’s angular momentum is decreasing
and becomes increasingly more negative (in order to keep $\Omega_{h}={\rm
const}$) consistent with finding (5). On the other hand, the apparent horizon
area (or the irreducible mass) is increasing, consistent with finding (3), and
the fact that $\kappa_{h}=4\pi r_{h}/A_{h}$.

###  3.4 Black hole-toroid in equilibrium

The first full calculation of a self-gravitating thick disk in general
relativity has been performed by Nishida & Eriguchi
NishidaEriguchi94:ToroidAroundBH . The starting point of their work was the
Bardeen formalism and its horizon boundary conditions. The new ingredients
were: 1) The stress energy tensor was assumed to be a perfect fluid, Eq. (9),
with a polytropic EOS $p=k\epsilon^{\Gamma}$777Notice that in most works
reported in this review the polytropic EOS is $p=k\rho^{\Gamma}$.. Contrary to
Lanza Lanza1992a , the authors assumed $p\neq 0$ and computed the hydrostatic
equilibrium of geometrically thick disks. 2) The differential rotation law for
the disk was assumed to be $u^{t}u_{\phi}=A^{2}(\Omega_{c}-\Omega)$ folllowing
the KEH Komatsu89 ; Komatsu89b  works on rotating stars. Here $A$ is an input
constant parameter that determines the degree of differential rotation and
$\Omega_{c}$ a constant that is evaluated during the iteration scheme. 3) The
numerical solution of the elliptic equations was performed using the KEH
method where the second order operator (a Laplacian) is inverted by employing
a Green’s function approach.

In general, each equilibrium solution of a BHD is obtained by specifying the
following: i) the black hole mass, ii) the black hole spin, iii) the EOS of
the toroid, iv) the rotation law of the toroid, v) the rest mass of the
toroid, vi) the position of the toroid relative to the black hole, vii) the
total angular momentum of the toroid. For a strict equilibrium, the directions
of black hole spin and the disk’s angular momentum should be aligned or
antialigned. Properties i), ii) are free parameters, and similarly functional
relations iii) and iv) can be chosen freely. On the other hand properties
v)-vii) determine a unique BHD model. In actual numerical computations
different authors use different ways to specify v)-vii). For example property
v) is controlled by specifying the maximum density $\rho_{\rm max}$ (or
$\epsilon_{\rm max}$). For properties vi) and vii) Nishida & Eriguchi choose
to fix the inner and the outer point of the disk.

The authors constructed sequences of BHDs around nonspinning as well as mildly
rotating black holes ($J_{h}/M_{h}^{2}\lessapprox 0.6$) for $\Gamma=2$, and
$\Gamma=5/3$ polytropes. In general they found that the disk plays the role of
an “anchor” since a larger angular momentum is needed to keep the same angular
velocity as the corresponding Kerr black hole. They paid special attention to
the cases with $\omega_{h}=0$, and confirmed previous results Lanza1992a ;
Will1974 ; Will1975  that even when $\ell\neq\rm const$ the self-gravity of
the torus induces negative angular momentum $J_{h}$ on the black hole.
Alternatively one can have a black hole with $J_{h}=0$, while $\omega_{h}\neq
0$. The authors find that black holes with $J_{h}=0$ have ratio of polar to
equatorial proper circumference $C_{p}/C_{e}$ always equal to unity
irrespective of the mass of the disk. This motivates them to define a
“nonrotating” black hole as one that has $J_{h}=0$ or $C_{p}/C_{e}=1$. Since
those equilibria will have $\omega_{h}\neq 0$, they will exhibit ergoregions
(by definition) despite being spherical in shape. In other words in the
presence of self-gravitating disks one can imagine energy extraction even from
“nonrotating” black holes.

Ansorg & Petroff Ansorg05  using their highly accurate multidomain
pseudospectral code Ansorg:2003br , and the Bardeen formalism
Bardeen73_Stars_Disks_Black_Holes , constructed a wide range of uniformly
rotating and constant density self-gravitating rings. They were able to
investigate various properties up to machine precision (for example the
authors claim that when $J_{h}=0$ the ratio $C_{p}/C_{e}$ is close but not
exactly one). They presented a BHD with $J_{h}/M_{h}^{2}>1$, while in a
follow-up work Ansorg2006  they calculate a model with $M_{h}<0$, lending more
evidence that the Komar mass for the black hole in the presence of matter is
not a good diagnostic.

More recently, the method of Nishida & Eriguchi was employed by Stergioulas
Stergioulas2011 ; Stergioulas2011c  in order to produce more accurate
solutions using a finite difference code. In particular he applied a
compactification in the numerical domain, as in RNS calculations Cook94b ,
using a redefinition of the radial coordinate. Self-gravitating sequences of
BHD were produced with vanishing horizon angular velocity and it was found
that self-gravitating heavy tori with constant specific angular momentum can
fill their Roche lobe only if $\ell<4M_{\rm bh}$, similar to massless disks
around a Schwarzschild black hole.

###  3.5 Black hole-toroid in the puncture framework

New ideas for the construction of self-gravitating BHDs have been introduced
by Shibata Shibata2007 , where methods from binary black hole initial data
calculations were employed. In particular the author uses the 3+1 formalism
and the puncture framework with inversion symmetric boundary conditions at the
black hole throat, together with assumptions on stationarity and axisymmetry
to derive a new set of equations (4 elliptic plus one first order) that solve
the Einstein system. His starting point was the formulation of Krivan & Price
Krivan:1998td  according to which in axisymmetric spacetimes around rotating
black holes, a nonconformally flat form of the 3-geometry can be chosen which
allows a simple superposition of Kerr black holes with arbitrary mass and
spin. For a spacetime element in quasi-isotropic form, Eq. (56), the
nontrivial components of the extrinsic curvature are

$K_{i\varphi}=\frac{\psi^{4}}{2\alpha}r^{2}\sin^{2}\theta\partial_{i}\beta,\qquad\mbox{where}\qquad
i\in\\{r,\theta\\},$ (68)

and the slicing is maximal, i.e. $K=0$. As a consequence there is only one
component for the momentum constraint (the $\varphi$ component) that needs to
be satisfied and results to

$\frac{1}{r^{2}}\partial_{r}(r^{2}\psi^{2}K_{r\varphi})+\frac{1}{r^{2}\sin\theta}\partial_{\theta}(\sin\theta\psi^{2}K_{\theta\varphi})=8\pi
T^{t}_{\ \varphi}\alpha\psi^{6}e^{2q}\,.$ (69)

The equation above is linear in the conformal frame, i.e., in
$\psi^{2}K_{i\varphi}$ Krivan:1998td . Thus one can make a decomposition that
separates the contribution of the Kerr black hole from the torus as

$\psi^{2}K_{i\varphi}:=K^{\scriptscriptstyle
K}_{i\varphi}+K^{\scriptscriptstyle
T}_{i\varphi},\qquad\mbox{and}\qquad\beta:=\beta_{\scriptscriptstyle
K}+\beta_{\scriptscriptstyle T}$ (70)

with

$K^{\scriptscriptstyle
K}_{i\varphi}:=\frac{H_{i}\sin\theta}{r},\qquad\mbox{and}\qquad
K^{\scriptscriptstyle
T}_{i\varphi}:=\frac{\psi^{6}}{2\alpha}r^{2}\sin^{2}\theta\partial_{i}\beta_{\scriptscriptstyle
T}.$ (71)

with $H_{i}=(H_{r},H_{\theta})=(H_{E}\sin\theta/r,H_{F})$ and $H_{E}$, $H_{F}$
having well-known expressions in Kerr spacetime Krivan:1998td . The “Kerr”
contribution to the shift is calculated from

$\partial_{r}\beta_{\scriptscriptstyle K}:=\frac{2\alpha
H_{E}}{\psi^{6}r^{4}},$ (72)

although this does not mean that $\beta_{\scriptscriptstyle K}$ will be the
Kerr shift, since the conformal factor $\psi$ in Eq. (72) will be a solution
of the Hamiltonian constraint and will have a contribution from the torus.
Equation (70) and the momentum constraint will then yield an elliptic equation
for the torus shift $\beta_{\scriptscriptstyle T}$. Elliptic equations for the
lapse $\alpha$, and the conformal factor $\psi$, are derived from the
$\partial_{t}K=0$ equation, and the Hamiltonian constraint respectively. On
the other hand a combination of the $\partial_{t}K_{ij}=0$ equations,

$\gamma^{rr}\partial_{t}K_{rr}+\gamma^{\theta\theta}\partial_{t}K_{\theta\theta}-3\gamma^{\phi\phi}\partial_{t}K_{\phi\phi}-\partial_{t}K=0,$
(73)

will yield an elliptic equation for $q$. Together with Eq. (72), this new
system in 3+1 and axisymmetry has one more elliptic equation (in the absence
of a torus one has 4 elliptic equations for $\alpha,\beta,\psi,q$) than the
Bardeen-Wagoner system Bardeen71 ; Bardeen73_Stars_Disks_Black_Holes ;
Butterworth1975  and completely determines the Einstein equations.

Shibata introduces the puncture by a transformation for the conformal factor
and the lapse to a new set of variables $s,\,B$

$\psi:=\left(1+\frac{r_{h}}{r}\right)e^{s},\qquad\alpha\psi:=\left(1-\frac{r_{h}}{r}\right)e^{-s}B,$
(74)

where $r_{h}=\sqrt{m^{2}-a^{2}}/2$ is the black hole horizon in quasi-
isotropic coordinates, and $m$, $a$ the mass and spin parameter of the initial
Kerr black hole. This transformation leads to a new elliptic system in terms
of $\\{s,B,\beta_{\scriptscriptstyle T},q\\}$. Because the vacuum equations
are inversion symmetric with respect to the $r=r_{h}$ surface, Shibata extends
that symmetry in the presence of the fluid and imposes

$\partial_{r}s=\partial_{r}B=\partial_{r}\beta_{\scriptscriptstyle
T}=\partial_{r}q=0,\qquad\mbox{at}\qquad r=r_{h}.$ (75)

With the above boundary conditions the 2-surface $r=r_{h}$ becomes a
marginally outer trapped surface, or an apparent horizon, which in stationary
spacetimes agrees with the event horizon Hawking73a .

Another important contribution of this work is the careful study of the
diagnostics for the BHD that corrected many of the problems that plagued
previous works, and in particular the diagnostic for the black hole mass.
Contrary to the angular momentum which can unambiguously be separated into the
black hole and the torus components, Eq. (64), the case of mass needs more
care. A similar separation for the mass, Eq. (65), does not lead to a good
measure of the black hole mass $M_{h}$ in the presence of a massive disk. This
point is already being discussed by Bardeen Bardeen73_Stars_Disks_Black_Holes
who notes that $M_{h}$ defined as such is smaller than the true value of the
black hole mass, by an amount first order in the mass of the torus, since it
includes the binding energy of the system. Numerically this was already
reported in the self-gravitating BHD computed by Lanza Lanza1992a  who found
$M_{h}$ to be decreasing as the rest mass of the torus increased. Assigning
$M_{h}$ to the black hole mass is problematic and Shibata showed explicitly
the following erroneous conclusions that such a diagnostic entails: i) For
heavy tori $M_{h}$ becomes smaller than the irreducible mass of the black hole
(which is impossible by definition of $M_{\rm irr}$) and ii) The sum of the
rest mass of the torus and $M_{h}$ is smaller than the ADM mass of the system
(which is is impossible for a bound system). In order to remedy this problem
the author proposed two new diagnostics for the black hole mass both of which
are inspired from the isolated Kerr spacetime. In the first one the black hole
mass is estimated from $M_{\scriptscriptstyle C}:=C_{e}/(4\pi)$, where $C_{e}$
the proper equatorial circuference of the black hole horizon. In the second
one, since the angular momentum of the black hole is $J_{h}$ (Eq. (64)) and
the irreducible mass can be found from the apparent horizon area $M_{\rm
irr}:=\sqrt{A_{h}/(16\pi)}$, Shibata uses the Christodoulou formula
Christodoulou70 , $M_{\rm bh}:=M_{\rm irr}\sqrt{1+J_{h}^{2}/(4M_{\rm
irr}^{4})}$, to get another measure for the black hole mass. For disk masses
up to twice the black hole mass the two diagnostics, $M_{\scriptscriptstyle
C}$ and $M_{\rm bh}$ agree to $\sim 3\%$.

Table 1: Dependence of various black hole diagnostics as a function of the disk rest mass, $M_{0}$, along a sequence of $j=\rm const$ self-gravitating BHDs. $M_{\rm bh}\uparrow$ | $M_{\scriptscriptstyle C}\uparrow$ | $M_{h}\downarrow$ | $M_{\rm irr}\uparrow$ | $\frac{|J_{h}|}{M_{\rm bh}^{2}}\downarrow$ | $\frac{C_{p}}{C_{e}}\uparrow$ | $|M_{\rm bh}\Omega_{h}|\downarrow$ | $M_{\rm bh}\kappa_{h}\downarrow$
---|---|---|---|---|---|---|---

Using a $\Gamma=4/3$ polytropic EOS and a $j=hu_{\phi}={\rm const}$
differential rotating law the author calculates self-gravitating sequences of
disks around highly spinning black holes. The author fixes the mass $m$, and
the spin parameter $a$ of the Kerr black hole, therefore the radius of the
throat is also fixed. As the mass of the toroid increases, the mass and spin
of the black hole diverge from the fixed Kerr quantities $m$ and $a$. Also
because of the boundary condition on the shift $\beta$ at the throat,
$J_{h}=ma$ always, irrespective of the mass and the spin of the torus. In
Table 1 we summarize the behavior of various diagnostics, as the mass of the
disk increases. In particular:

  1. 1.

The angular velocity of the black hole horizon, $\Omega_{h}$, decreases with
increasing disk mass, due to the fact that the black hole frame-dragging is
reduced by the slower rotation of the disk (the disk plays the role of an
“anchor” Lanza1992a ; NishidaEriguchi94:ToroidAroundBH ).

  2. 2.

The strength of gravity on the event horizon and the magnitude of its surface
gravity $\kappa_{h}$ is weakened by the tidal force of the torus (less angular
momentum is needed to stay in a Keplerian orbit) Lanza1992a ;
NishidaEriguchi94:ToroidAroundBH ).

  3. 3.

The area of the black hole horizon and therefore its irreducible mass is
increasing. This results to an actual increase in the black hole mass as
measured by the diagnostic $M_{\rm bh}$, since $J_{h}$ is constant. At the
same time the dimensioless spin decreases with the disk mass, consistent with
finding (1).

  4. 4.

Also consistent with the behavior of the dimensioless spin is the fact that
the ratio $C_{p}/C_{e}$ increases as the torus becomes heavier, and the black
hole becomes more spherical.

More recently the method of Shibata has been employed by the Kraków group in a
series of works that investigate the properties of rotation around self-
gravitating axisymmetric disks as well as to construct magnetized equilibria
using a second order finite difference code Karkowski2017 ; Karkowski2018 ;
Mach2019 ; Kulczycki2019 ; Dyba2019 ; Kulczycki2021 . In particular Karkowski
et al.  Karkowski2017 ; Karkowski2018  proposed a new Keplerian rotation law
that holds for massless as well as massive disks according to

$u^{t}u_{\phi}=-\frac{3}{2\lambda}\frac{d}{d\Omega}\ln\left[1-\frac{\lambda}{3}[a^{2}\Omega^{2}+3w^{\frac{4}{3}}\Omega^{\frac{2}{3}}(1-a\Omega)^{\frac{4}{3}}]\right].$
(76)

Here $\lambda$ is a parameter that takes values around $\lambda_{\rm Kerr}=3$,
which corresponds to the exact formula that characterizes the motion of
circular geodesics at the equatorial plane of a Kerr black hole with spin
parameter $a$ and mass $m=w^{2}$. For self-gravitating tori $w^{2}\neq m$ in
general. In the Newtonian limit, Eq. (76) yields the Keplerian angular
velocity $\Omega=w/(r\sin\theta)^{3/2}$. Using the rotation law Eq. (76), and
the Shibata formulation, the authors constructed sequences of stationary and
axisymmetric equilibria around nonrotating as well as spinning black holes. In
Newtonian gravity von Zeipel’s theorem vonZeipel1924  states that for a
barotropic fluid the surfaces of constant $\Omega$ coincide with the surfaces
of constant $\ell_{N}=\Omega\varpi^{2}$, where $\varpi$ the distance from the
axis of rotation. This implies that $\ell_{N}=\ell_{N}(\varpi)$ or
$\Omega=\Omega(\varpi)$, i.e. the angular velocity depends only on the
distance from the axis of rotation (Poincaré-Wavre Tassoul-1978:theory-of-
rotating-stars ). In general relativity Abramowicz Abramowicz1974  showed that
for massless disks the surfaces of constant $\Omega$ have cylindrical
topology, therefore they depend not only on the distance from the rotation
axis but also on the distance from the equatorial plane of symmetry. As
expected this is also true for self-gravitating tori Karkowski2018 .

Many properties of self-gravitating BHDs have been investigated by Dyba et al.
Dyba2019  using the infrastructure developed in previous works. In particular
the authors find that by fixing the black hole parameters ($a$ and $m$), the
polytropic exponent $\Gamma$, and the inner and outer coordinate equatorial
radii of the torus, as well as its maximal rest-mass density, there exist two
solutions differing in the ADM mass. In other words, it is possible to obtain
a sequence of tori with a decreasing maximum rest-mass density and increasing
mass, simply because their size is also growing. Applying Seguin’s Seguin1975
criterion for linear stability against axially symmetric perturbations the
authors found that the massive branch must be dynamically unstable. The
location of the ISCO for a variety of self-gravitating equilibria showed non-
negligible differences from the corresponding Kerr value, even for light
toroids. A typical behavior is the increase of the circumferential radius of
the ISCO relative to that of the Kerr black hole. On the other hand Dyba et
al.  found that for sufficiently massive disks the effective potential
$\mathcal{V}_{\rm eff}$ due to its nonmonotonic behavior further from the
black hole, can exhibit a region outside the ISCO in which circular geodesics
can be unstable ($\mathcal{V}^{{}^{\prime\prime}}_{\rm eff}(r)<0$). This
provides yet another reason why very massive disks are dynamically unstable.

###  3.6 Black hole-toroid with magnetic fields

Mach et al.  Mach2019  constructed self-gravitating BHDs with toroidal
magnetic fields in the general relativistic ideal magnetohydrodynamics (IMHD)
regime Komissarov2006a ; Shibata2007  under stationarity and axisymmetry. The
main difference here is the integral of the Euler equation, since for the
gravitational field equations the analysis follows previous works Shibata2007
; Karkowski2017 ; Karkowski2018  (apart from the modification of the sources
in the Einstein equations). Mach et al.  assume that
$T^{\alpha\beta}=T^{\alpha\beta}_{\rm\scalebox{0.6}{M}}+T^{\alpha\beta}_{\rm\scalebox{0.6}{EM}}$
with

$T^{\alpha\beta}_{\rm\scalebox{0.6}{EM}}=b^{2}(u^{\alpha}u^{\beta}+\frac{1}{2}g^{\alpha\beta})-b^{\alpha}b^{\beta},$
(77)

the IMHD stress-energy tensor, and $T^{\alpha\beta}_{\rm\scalebox{0.6}{M}}$
the perfect fluid stress-energy tensor Eq. (9). Here
$b^{\alpha}=B^{\alpha}/\sqrt{4\pi}$ is the magnetic field measured by an
observer with 4-velocity $u^{\alpha}$ (at rest with respect to the fluid), and
$b^{2}=b_{\mu}b^{\mu}$ the magnetic pressure $p_{\rm mag}=b^{2}/2$. Since

$b_{\mu}u^{\mu}=0\,,$ (78)

assuming a toroidal magnetic field with $b^{r}=b^{\theta}=0$ yields the
following realation:

$b_{t}=-\Omega b_{\phi}.$ (79)

In this set up, the continuity equation, the induction equation, as well as
the $t$ and $\phi$ components of $\nabla_{\alpha}T^{\alpha\beta}=0$ are
trivially satisfied, while the $r$ and $\theta$ components satisfy

$\nabla_{\mu}\ln\frac{h}{u^{t}}+u^{t}u_{\phi}\nabla_{\mu}\Omega+\frac{1}{2\rho
h\mathcal{A}}\nabla_{\mu}(b^{2}\mathcal{A})=0,$ (80)

where $\mathcal{A}=g_{\phi\phi}g_{tt}-g_{t\phi}^{2}$. The above equation can
be integrated if and only if $u^{t}u_{\phi}$ is a function of $\Omega$, and
$b^{2}\mathcal{A}$ a function of $\rho h\mathcal{A}$. When the magnetic field
is zero one recovers the integral of the Euler Eq. (62). The authors assumed
$b^{2}|\mathcal{A}|=f(\rho h|\mathcal{A}|)$ with
$f^{\prime}(x)=2nC_{1}x/(1+C_{1}x)$, so that

$\int\frac{d(b^{2}\mathcal{A})}{2\rho h\mathcal{A}}=\ln(1-C_{1}\rho
h\mathcal{A})^{n}\,.$ (81)

Constants $n,\,C_{1}$ characterize the topology of the magnetic field and the
authors construct a series of solutions with variable $C_{1}$ while fixing
$n=1$. They found that the larger the $C_{1}$, the smaller the thermal
pressure, even in cases in which the maximum of the baryonic density is fixed.
In general they also observe an increase of the value of the magnetic
pressure. For the rotation law ($u^{t}u_{\phi}$ term) they employed Eq. (76)
with $\lambda=3$. Each solution is specified by the black hole parameters ($m$
and $a$), the inner $r_{\rm in}$ and outer $r_{\rm out}$ radii of the disk,
the polytropic exponent of the EOS, the maximum rest-mass density, and the
magnetic field constants $C_{1}$ and $n$. On the other hand, the constant $w$
that appears in the differential rotation law, Eq. (76), the constant that
appears in the integral of the Euler Eq. (80), and the two angular velocities
at the inner and outer points are determined during the iteration scheme. As
in Dyba et al.  Dyba2019  fixing $a,\,m,\,r_{\rm in},\,r_{\rm
out},\,\Gamma,\,\rho_{\rm max}$ does not lead to unique solutions. Two
solutions exist with different ADM masses, with one of them being larger than
the mass of the central black hole. The authors measure the magnetization of
their solutions by the parameter $\beta_{\rm mag}=2p/b^{2}$ and they construct
a large set of solutions with $\beta_{\rm mag}$ as small as $\sim 6\times
10^{-4}$. The main characteristic of the toroidal magnetic field is to shift
the location of the maximum rest-mass density towards the black hole.

###  3.7 Arbitrary spinning black hole-toroid

In an effort to go beyond axisymmetry or stationarity Tsokaros et al.
Tsokaros2018a  presented a new formalism for 3-dimensional self-gravitating
BHD solutions. This method goes beyond the minimal construction of initial
data (solution of constraints in binary black hole calculations) York71 ;
York72 ; Pfeiffer:2002iy ; Brandt97b ; York99 ; Isenberg08 ; Wilson95 ;
Marronetti99  and solves the full Einstein system. In particular the 5
equations related to the conformal geometry which are associated with the true
dynamical degrees of freedom of the gravitational field York72  are resolved.
The following suite of benchmarks have been used to assess the new formalism:
1) In the absence of matter it can reproduce the exact Kerr-Schild solution,
even for extreme spins, which makes the method suitable for a broad range of
nonaxisymmetric problems, such as tilted disks or binary systems. 2) The
domain of the solution extends inside the apparent horizon, which is well-
suited for evolution simulations. 3) In the presence of massless disks around
the black hole our method reproduces well-known solutions (e.g.
Chakrabarti1985 ; DeVilliers03b ) even with black hole tilt. 4) The first
self-consistent, self-gravitating, and tilted BHD solutions are presented that
satisfy not only the constraint equations but the whole Einstein system.

The starting point of Tsokaros et al.  Tsokaros2018a  is a line element in its
general 3+1 form (Eq. (2)) together with a conformal metric decomposition as
in Eq. (50). Instead of a maximal slicing and the Dirac gauge Eq. (51), the
coordinates now used were Kerr-Schild with

$K=K_{\rm\scriptscriptstyle KS}=\frac{2H\alpha^{3}_{\rm\scriptscriptstyle
KS}}{r}\left(1+H+\frac{2H^{2}r}{m}\right),$ (82)

and

${\raise 4.30554pt\hbox{${}^{\
\circ}$}}\\!\\!\\!\\!\\!D_{i}\tilde{\gamma}^{ij}={\raise 4.30554pt\hbox{${}^{\
\circ}$}}\\!\\!\\!\\!\\!D_{i}h^{ij}_{\rm\scriptscriptstyle KS}.$ (83)

Here $\alpha_{\rm\scriptscriptstyle KS},\,h^{ij}_{\rm\scriptscriptstyle KS}$
are the lapse and non flat part of the conformal metric in Kerr-Schild
coordinates $ds^{2}=(\eta_{\mu\nu}+2H\ell_{\mu}\ell_{\nu})dx^{\mu}dx^{\nu}$.
Assuming Kerr-Schild boundary conditions for the potentials
$\psi,\,\alpha,\,\tilde{\beta}_{i},\,h_{ij}$ the authors found that the
corresponding system Shibata04a ; Bonazzola:2003dm  could not converge in the
neighborhood outside the horizon. The failure was due to the fact that the
equations for $h_{ij}$ were losing their elliptic character close to the
horizon. The region of nonconvergence was larger for higher spins but still
relatively small (less than twice the horizon radius). In order to evercome
this difficulty, and even more to obtain horizon penetrating solutions the
authors introduced a new decomposition for the traceless extrinsic curvature
as

$\tilde{A}_{ij}=\tilde{A}^{\rm\scriptscriptstyle
KS}_{ij}+\tilde{\sigma}(\tilde{\mathbb{L}}\tilde{W})_{ij},$ (84)

where $\tilde{A}^{\rm\scriptscriptstyle KS}_{ij}$ is the Kerr-Schild part,
$\tilde{W}_{i}$ an unknown spatial vector, and $\tilde{\sigma}$ a scalar. Here
$\tilde{\mathbb{L}}$ is the conformal Killing operator:
$(\tilde{\mathbb{L}}\tilde{W})_{ij}=\tilde{D}_{i}\tilde{W}_{j}+\tilde{D}_{j}\tilde{W}_{i}-\frac{2}{3}\tilde{\gamma}_{ij}\tilde{D}_{k}\tilde{W}^{k}$.
Equation (84) implies that the momentum constraint is now solved twice once
for the shift vector $\tilde{\beta}_{i}$ and once for $\tilde{W}_{i}$. Solving
for $\beta^{i}$ is necessary since it will be used in the computation of
$\mathcal{L}_{\alpha\mathbf{n}}K$. On the other hand decomposition (84) with
zero boundary conditions for $\tilde{W}_{i}$ on the horizon results into
convergence for the $h_{ij}$ potentials even inside the horizon. The Kerr-
Schild gauge conditions Eq. (83) are satisfied through a transformation of the
form Eq. (53) which introduces 3 more potentials $\xi^{i}$. The augmented
system of the 17 elliptic equations with zero boundary conditions for the
gauge potentials and the vector $\tilde{W}_{i}$ converges smoothly in vaccum
or in the presence of matter (like a massive disk), even for near maximally-
spinning black holes, thus not only solving for the constraint equations but
providing also a way to control the gravitational wave content of the initial
data in a self-consistent way.

For the Euler equations the authors assumed a stationary and axisymmetric
fluid flow which in the presence of a tilted black hole (with respect to the
angular momentum of the disk) is a valid approximation if the inner point of
the disk is further away from the black hole horizon. Using the KEH scheme for
black holes Tsokaros2007 , as implemented within the cocal code Uryu2012  the
authors computed self-gravitating tilted BHD solutions with a disk mass larger
than the black hole mass and with the black hole having almost extremal spin.

##  4 Magnetized rotating neutron stars

While the typical surface magnetic field of pulsars is $\sim
10^{12}-10^{13}\;\rm G$ lyne_graham-smith_2012  there exist neutron stars, the
so-called magnetars, with extremely large magnetic fields $\gtrsim
10^{14}\;\rm G$ Duncan1992 ; ThompsonDuncan1995 ; ThompsonDuncan1996 . The
magnetar model has been invoked to explain soft $\gamma$-ray repeaters (object
that emits large bursts of $\gamma$-rays and X-rays at irregular intervals)
and the related anomalous X-ray pulsars Mazets1979 ; Paczynski92 ; Duncan1992
; ThompsonDuncan1995 ; ThompsonDuncan1996 ; Kouveliotou1994 ; Kouveliotou1998
; Kouveliotou1999 ; Mereghetti2008 . Soft $\gamma$-ray repeater bursts are
caused by the cracking of the neutron star crust due to magnetic stresses,
which leads to injection of Alfvén waves into the magnetosphere, particle
acceleration, and formation of an optically thick pair plasma. The decay of
the magnetic field heats the neutron-star interior, and gives rise to
persistent thermal soft X-ray emission from the surface. In addition, the
strong magnetic field causes the spindown of the neutron star which may end up
with rotation periods of $\sim 10\;\rm s$. From the microscopic point of view
when the cyclotron energy equals the electron rest mass, one reaches the
quantum critical magnetic field strength of $\sim 10^{13}\;\rm G$ beyond which
the magnetic field affects physical processes and is responsible for many
exotic phenomena, such as vacuum birefringence, photon splitting, and the
distortion of atoms (see Harding2006  for an extensive review).

Magnetic fields beyond magnetar strength that can reach values of $\sim
10^{17}\;\rm G$ can be developed in the merger of two neutron stars due to a
number of different mechanisms: (i) The Kelvin-Helmholtz instability
Chandrasekhar81 ; Bodo1994  which occurs in the shear layer that forms between
the two neutron stars and can grow on a timescale of a couple of milliseconds
Price06 ; Anderson2008 ; Kiuchi2015a . (ii) The magnetorotational instability
Velikhov1959 ; Chandrasekhar1960 ; Balbus1991  during the merger as well as in
the postmerger compact object Shibata:2005mz ; Duez2006a ; Siegel2013 ;
Kiuchi2015 . (iii) Magnetic winding which is due to differential rotation
Baumgarte00bb ; Shapiro00 ; Kiuchi2015a . Differential rotation generates
toroidal Alfvén waves which convert rotational kinetic energy into magnetic
field energy. In fact all of the mechanisms above are responsible for
converting poloidal to toroidal magnetic fields leading to a remnant with
comparable poloidal and toroidal components. Strong magnetic fields affect the
neutron star in at least two ways. First, they result in an anisotropy through
a modification of the energy-momentum tensor. Second, they affect the EOS due