# Numerical-relativity simulation for tidal disruption of white dwarfs by a
supermassive black hole

Alan Tsz-Lok Lam1 Masaru Shibata1,2 Kenta Kiuchi1,2 1Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam-
Golm 14476, Germany  2Center for Gravitational and Quantum-Information
Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto,
606-8502, Japan

###### Abstract

We study tidal disruption of white dwarfs in elliptic orbits with the
eccentricity of $\sim 1/3$–$2/3$ by a nonspinning supermassive black hole of
mass $M_{\rm BH}=10^{5}M_{\odot}$ in fully general relativistic simulations
targeting the extreme mass-ratio inspiral leading eventually to tidal
disruption. Numerical-relativity simulations are performed by employing a
suitable formulation in which the weak self-gravity of white dwarfs is
accurately solved. We reconfirm that tidal disruption occurs for white dwarfs
of the typical mass of $\sim 0.6M_{\odot}$ and radius $\approx 1.2\times
10^{4}$ km near the marginally bound orbit around a nonspinning black hole
with $M_{\rm BH}\lesssim 4\times 10^{5}M_{\odot}$.

###### pacs:

04.25.D-, 04.30.-w, 04.40.Dg

##  I Introduction

Tidal disruption of ordinary stars and/or white dwarfs by supermassive black
holes has been revealed to be one of the major sources of bright
electromagnetic transients (see, e.g., Refs. [1, 2, 3]), which have been
actively observed in the last decade. In addition, gravitational waves emitted
by tidal disruption of white dwarfs closely orbiting supermassive black holes
could be observable by Laser Interferometer Space Antenna (LISA) [4].
Electromagnetic signals associated with tidal excitation (e.g., Ref. [5]) or
mass stripping (e.g., Refs. [6, 7, 8, 9] for related works) or tidal
disruption (e.g., Refs. [10, 11]) of white dwarfs can be an important
electromagnetic counterpart of gravitational waves. Because the expected event
rate is not so high [12] that the signal-to-noise ratio of gravitational waves
for the LISA sensitivity is unlikely to be very high, the discovery of the
possible electromagnetic counterparts will help extract gravitational waves
from the noisy data in the LISA mission.

The condition for mass shedding and tidal disruption during the cross
encounter of stars with supermassive black holes is often described by the so-
called $\beta$-parameter defined by

$\beta:=\frac{r_{t}}{r_{p}},$ (1)

where $r_{p}$ is the periastron radius for the orbit and $r_{t}$ is the Hill’s
radius [13] defined by

$r_{t}:=R_{*}\left(\frac{M_{\rm BH}}{M_{*}}\right)^{1/3},$ (2)

with $R_{*}$ the stellar radius, $M_{*}$ the stellar mass, and $M_{\rm BH}$
the mass of the supermassive black hole, respectively. Since the early 1990s
(see, e.g., Refs. [14, 15, 16]), a large number of numerical simulations have
been performed in the last three decades (see, e.g., Refs. [17, 18] for
reviews of the latest works and Refs. [19, 20, 21, 22, 23, 24] for some of the
most advanced works). They have shown that mass stripping can take place at
the close encounter if $\beta$ is larger than about 0.5, and tidal disruption
can take place if $\beta\gtrsim 1$ for stars in parabolic orbits (see, e.g.,
Refs. [17, 25, 26] for Newtonian simulation works, and also early
semianalytical work [27, 28]). It is also shown that for close orbits around a
black hole, the general relativistic effect can significantly reduce the
critical value of $\beta$ for the tidal disruption [23]. Indeed, general
relativistic works show that for circular orbits near the innermost stable
circular orbit of black holes, the mass shedding can occur even for $\beta\sim
0.4$ [29, 30].

However, the previous analyses have been carried out in Newtonian gravity or
in relativistic gravity of a black hole with Newtonian (or no) gravity for the
companion star or in a tidal approximation with a relativistic tidal potential
[29, 31, 30]. To date, no fully general relativistic (the so-called numerical-
relativity) simulation, i.e., a simulation with no approximation except for
the finite differencing, has been done for the tidal disruption problem with
$\beta\lesssim 1$ (but see Refs. [32, 33, 34] for a head-on and an off-axis
collision).

Numerical-relativity simulation is suitable for the tidal disruption problem
for the case that the orbit at the tidal disruption is highly general-
relativistic. This is particularly the case for tidal disruption of white
dwarfs by supermassive black holes because it can occur only for orbits very
close to the black-hole horizon. Advantages of the numerical-relativity
simulation are: (i) the redistribution of the energy and angular momentum of
the star can be followed in a straightforward manner and (ii) we can directly
follow the matter motion after the tidal disruption including the subsequent
disk formation.

In this paper, we present a result of numerical-relativity simulations for
tidal disruption of white dwarfs of typical mass ($0.6$–$0.8M_{\odot}$) by a
supermassive black hole with relatively low mass ($M_{\rm
BH}=10^{5}M_{\odot}$) for the first time. For simplicity, the white dwarfs are
modeled by the $\Gamma=5/3$ polytropic equation of state. As a first step
toward more detailed and systematic studies, we focus on tidal disruption of
white dwarfs in mildly elliptic orbits aiming at confirming that our
numerical-relativity approach is suitable for reproducing the criteria of
tidal disruption, which has been already investigated in many previous works
referred to above.

The paper is organized as follows. In Sec. II, we describe our formulation for
evolving gravitational fields, matter fields, and for providing initial data
of a star in elliptic orbits around supermassive black holes. In Sec. III,
numerical results are presented paying particular attention to the criterion
for tidal disruption. Section IV is devoted to a summary. Throughout this
paper, we use the geometrical units of $c=1=G$ where $c$ and $G$ are the speed
of light and gravitational constant, respectively. The Latin and Greek indices
denote the space and spacetime components, respectively.

##  II Basic equations for the time evolution

###  II.1 Gravitational field

First, we reformulate the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism
[35, 36] in numerical relativity to a form suitable for the simulation of
high-mass ratio binaries, in particular for accurately computing a weak self-
gravity of white dwarfs. Throughout this paper, high-mass ratio binaries imply
those composed of a very massive black hole of mass $M_{\rm BH}\gtrsim
10^{5}M_{\odot}$ and a white dwarf (or an ordinary star) of mass of
$M_{*}=O(M_{\odot})$ with the radius $R_{*}\gtrsim 10^{3}$ km, for which the
compactness defined by $M_{*}/R_{*}$ is smaller than $10^{-3}$.

We consider the two-body problem with a compact orbit of the orbital
separation $r\lesssim 30M_{\rm BH}$. With such setting, the magnitude of the
gravitational field generated by the black hole, which is defined by
$g_{\mu\nu}-\eta_{\mu\nu}$ is, of order $M_{\rm BH}/r>10^{-2}$. Here
$g_{\mu\nu}$ and $\eta_{\mu\nu}$ are the spacetime metric and Minkowski
metric, respectively. On the other hand, the magnitude of the gravitational
field generated by white dwarfs and ordinary stars is of order
$M_{*}/R_{*}<10^{-3}$, which is much smaller than that by the black hole. To
accurately preserve the nearly equilibrium state of such stars during their
orbits, an accurate computation of the gravitational field by them is
required. However, if we simply solve Einstein’s equation, a numerical error
for the computation of the black-hole gravitational field can significantly
affect the gravitational field for the white dwarfs/ordinary stars. To avoid
this numerical problem, we separate out the gravitational field into the black
hole part and other part, although we still solve fully nonlinear equations.
The idea employed here is similar to that of Ref. [33], but we develop a
formalism based on the BSSN formalism.

In a version of the BSSN formalism [36], the basic equations are written in
the form:

$\displaystyle(\partial_{t}-\beta^{k}\partial_{k})\tilde{\gamma}_{ij}=-2\alpha\tilde{A}_{ij}+\tilde{\gamma}_{ik}\partial_{j}\beta^{k}+\tilde{\gamma}_{jk}\partial_{i}\beta^{k}-\frac{2}{3}\tilde{\gamma}_{ij}\partial_{k}\beta^{k},~{}~{}~{}~{}~{}$
(3)
$\displaystyle(\partial_{t}-\beta^{l}\partial_{l})\tilde{A}_{ij}=W^{2}\biggl{[}\alpha\left(R_{ij}-\frac{\gamma_{ij}}{3}R_{k}^{~{}k}\right)$
(5) $\displaystyle\hskip
56.9055pt-\left(D_{i}D_{j}\alpha-\frac{\gamma_{ij}}{3}D_{k}D^{k}\alpha\right)\biggr{]}$
$\displaystyle\hskip
56.9055pt+\alpha\left(K\tilde{A}_{ij}-2\tilde{A}_{ik}\tilde{A}_{j}^{~{}k}\right)$
$\displaystyle\hskip
56.9055pt+\tilde{A}_{kj}\partial_{i}\beta^{k}+\tilde{A}_{ki}\partial_{j}\beta^{k}-\frac{2}{3}\tilde{A}_{ij}\partial_{k}\beta^{k}$
$\displaystyle\hskip 56.9055pt-8\pi\frac{G}{c^{4}}\alpha
W^{2}\left(S_{ij}-\frac{1}{3}\gamma_{ij}S_{k}^{~{}k}\right),$
$\displaystyle(\partial_{t}-\beta^{l}\partial_{l})W=\frac{W}{3}\left(\alpha
K-\partial_{k}\beta^{k}\right),$
$\displaystyle(\partial_{t}-\beta^{l}\partial_{l})K=\alpha\left[\tilde{A}_{ij}\tilde{A}^{ij}+\frac{1}{3}K^{2}\right]$
$\displaystyle\hskip
62.59596pt-W^{2}\left(\tilde{D}_{k}\tilde{D}^{k}\alpha-\frac{\partial_{i}W}{W}\tilde{\gamma}^{ij}\partial_{j}\alpha\right)$
$\displaystyle\hskip 62.59596pt+4\pi\frac{G}{c^{4}}\alpha\left(\rho_{\rm
h}+S_{k}^{~{}k}\right),$ (6)
$\displaystyle(\partial_{t}-\beta^{l}\partial_{l})\tilde{\Gamma}^{i}=-2\tilde{A}^{ij}\partial_{j}\alpha+2\alpha\left[\tilde{\Gamma}^{i}_{~{}jk}\tilde{A}^{jk}\right.$
$\displaystyle\hskip
42.67912pt\left.-\frac{2}{3}\tilde{\gamma}^{ij}\partial_{j}K-8\pi\frac{G}{c^{4}}\tilde{\gamma}^{ik}J_{k}-3\frac{\partial_{j}W}{W}\tilde{A}^{ij}\right]$
$\displaystyle\hskip
34.14322pt-\tilde{\Gamma}^{j}\partial_{j}\beta^{i}+\frac{2}{3}\tilde{\Gamma}^{i}\partial_{j}\beta^{j}+\frac{1}{3}\tilde{\gamma}^{ik}\partial_{k}\partial_{j}\beta^{j}+\tilde{\gamma}^{jk}\partial_{j}\partial_{k}\beta^{i},$
(7)

where $\alpha$ is the lapse function, $\beta^{j}$ is the shift vector,
$\tilde{\gamma}_{ij}$ is the conformal three metric defined from the three
metric $\gamma_{ij}$ by $\tilde{\gamma}_{ij}:=\gamma^{-1/3}\gamma_{ij}$ with
$\gamma={\rm det}(\gamma_{ij})$, $W:=\gamma^{-1/6}$, $\tilde{A}_{ij}$ is the
conformal trace-free extrinsic curvature defined from the extrinsic curvature
$K_{ij}$ by $\tilde{A}_{ij}=W^{2}(K_{ij}-\gamma_{ij}K/3)$ with
$K:=K_{k}^{~{}k}$, and $\tilde{\Gamma}^{i}:=-\partial_{j}\tilde{\gamma}^{ij}$.
$\rho_{\rm h}$, $J_{i}$, and $S_{ij}$ are quantities defined from the energy-
momentum tensor, $T_{\mu\nu}$, by $\rho_{\rm h}=T^{\mu\nu}n_{\mu}n_{\nu}$,
$J_{i}=-T^{\mu\nu}n_{\mu}\gamma_{\nu i}$, and $S_{ij}=T^{\mu\nu}\gamma_{\mu
i}\gamma_{\nu j}$ with $n^{\mu}$ the timelike unit vector normal to spatial
hypersurfaces.

In this problem, we employ the so-called puncture gauge [37], in which the
evolution equations for $\alpha$ and $\beta^{i}$ are written as

$\displaystyle\partial_{t}\alpha$ $\displaystyle=$ $\displaystyle-2\alpha K,$
(8) $\displaystyle\partial_{t}\beta^{i}$ $\displaystyle=$
$\displaystyle\frac{3}{4}B^{i},$ (9) $\displaystyle\partial_{t}B^{i}$
$\displaystyle=$ $\displaystyle\partial_{t}\tilde{\Gamma}^{i}-\eta_{B}B^{i},$
(10)

where $B^{i}$ is an auxiliary three-component variable and $\eta_{B}$ is a
constant of order $M_{\rm BH}^{-1}$.

By introducing a static black-hole solution for the geometric variables,
$\alpha_{0}$, $\beta_{0}^{i}$, $\tilde{\gamma}^{0}_{ij}$, $W_{0}$,
$\tilde{A}^{0}_{ij}$, and $K_{0}$ and by writing all the variables by

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle\alpha_{0}+\alpha_{\rm
s},$ (11) $\displaystyle\beta^{i}$ $\displaystyle=$
$\displaystyle\beta_{0}^{i}+\beta_{\rm s}^{i},$ (12)
$\displaystyle\tilde{\gamma}_{ij}$ $\displaystyle=$
$\displaystyle\tilde{\gamma}^{0}_{ij}+\tilde{\gamma}^{\rm s}_{ij},$ (13)
$\displaystyle W$ $\displaystyle=$ $\displaystyle W_{0}+W_{\rm s},$ (14)
$\displaystyle\tilde{A}_{ij}$ $\displaystyle=$
$\displaystyle\tilde{A}^{0}_{ij}+\tilde{A}^{\rm s}_{ij},$ (15) $\displaystyle
K$ $\displaystyle=$ $\displaystyle K_{0}+K_{\rm s},$ (16)
$\displaystyle\tilde{\Gamma}^{i}$ $\displaystyle=$
$\displaystyle\tilde{\Gamma}_{0}^{i}+\tilde{\Gamma}_{{\rm s}}^{i},$ (17)

we then write down the equations for $\alpha_{\rm s}$, $\beta_{\rm s}^{i}$,
$\tilde{\gamma}^{\rm s}_{ij}$, $W_{\rm s}$, $\tilde{A}^{\rm s}_{ij}$, $K_{\rm
s}$, and $\tilde{\Gamma}_{\rm s}^{i}$ (these are denoted by a representative
variable $Q_{s}$ as follows). Specifically, the evolution equations (5)–(7)
and (8)–(10) of the geometrical variables (denoted by a representative
variable $Q$) are schematically written in the form

$\partial_{t}Q=F(Q).$ (18)

Then, for the decomposition of $Q=Q_{0}+Q_{\rm s}$ with $F(Q_{0})=0$ (under
the conditions of $\partial_{t}Q_{0}=0$), we write the equation for $Q_{\rm
s}$ as

$\partial_{t}Q_{\rm s}=F(Q_{0}+Q_{\rm s})-F(Q_{0}).$ (19)

In numerical simulation, $F(Q_{0})$ obtained from finite difference is nonzero
which contains the truncation error of evolving the stationary background
metric numerically. Here, we added the second-term in the right-hand side to
explicitly subtract the leading error of evolving the background metric so
that the right-hand side of the evolution equation of $Q_{\rm s}$ does not
have the zeroth order terms in $Q_{\rm s}$

Any static black-hole solutions can be used for $\alpha_{0}$,
$\beta_{0}^{i},\cdots$, but in the BSSN formalism with the puncture gauge, the
metric relaxes to a solution in the limit hypersurface with $K_{0}=0$. Using
such a trumpet-puncture black hole also allows us to construct the initial
data in the conformal-thin-sandwich (CTS) formalism [38] (see Sec. II.3).
Thus, in the present formalism, it is appropriate to employ such a solution.
In the nonspinning black hole, the analytic solution is known and is written
as [39]

$\displaystyle\alpha_{0}$ $\displaystyle=$ $\displaystyle\sqrt{1-\frac{2M_{\rm
BH}}{R}+\frac{27M_{\rm BH}^{4}}{16R^{4}}},$ (20) $\displaystyle\beta_{0}^{i}$
$\displaystyle=$ $\displaystyle\frac{3\sqrt{3}M_{\rm BH}^{2}}{4R^{3}}x^{i},$
(21) $\displaystyle W_{0}$ $\displaystyle=$ $\displaystyle\frac{r}{R},$ (22)
$\displaystyle\tilde{\gamma}^{0}_{ij}$ $\displaystyle=$
$\displaystyle\delta_{ij},~{}~{}{\rm i.e.},~{}~{}\tilde{\Gamma}_{0}^{i}=0,$
(23) $\displaystyle\tilde{A}^{0}_{ij}$ $\displaystyle=$
$\displaystyle\frac{3\sqrt{3}M_{\rm
BH}^{2}}{4R^{3}}\left(\delta_{ij}-3\frac{x^{i}x^{j}}{r^{2}}\right),$ (24)

and $K_{0}=0$ where $R$ is a function of $r$ determined by [40]

$\displaystyle r$ $\displaystyle=$ $\displaystyle\left(\frac{2R+M_{\rm
BH}+\sqrt{4R^{2}+4M_{\rm BH}R+3M_{\rm BH}^{2}}}{4}\right)$ (25)
$\displaystyle\times\left[\frac{(4+3\sqrt{2})(2R-3M_{\rm BH})}{8R+6M_{\rm
BH}+3\sqrt{8R^{2}+8M_{\rm BH}R+6M_{\rm BH}^{2}}}\right]^{1/\sqrt{2}}.$

We note that $r=0$ corresponds to $R=3M_{\rm BH}/2$ and the event horizon is
located at $R=2M_{\rm BH}$ (i.e., $r\approx 0.78M_{\rm BH}$) in this solution.

###  II.2 Hydrodynamics

In this paper we model white dwarfs simply by the polytropic equation of
state,

$P=\kappa\rho^{\Gamma},$ (26)

where $P$ and $\rho$ are the pressure and rest-mass density, respectively,
$\kappa$ the polytropic constant, and $\Gamma$ adiabatic index for which we
set to be $5/3$. For the hydrodynamics, we solve the continuity and Euler
equations,

$\displaystyle\nabla_{\mu}(\rho u^{\mu})=0,$ (27)
$\displaystyle\nabla_{\mu}T^{\mu}_{~{}k}=0,$ (28)

with $\nabla_{\mu}$ the covariant derivative with respect to $g_{\mu\nu}$ and

$T_{\mu\nu}=(\rho+\rho\varepsilon+P)u_{\mu}u_{\nu}+Pg_{\mu\nu},$ (29)

where $\varepsilon$ and $u^{\mu}$ are the specific internal energy and four
velocity, respectively. In this work we do not solve the energy equation, and
determine $\varepsilon$ simply by
$\varepsilon=\kappa\rho^{\Gamma-1}/(\Gamma-1)$ which is derived from the
condition that the specific entropy is conserved for the fluid elements. The
continuity and Euler equations are solved in the same scheme as that used in
Refs. [41, 42].

The motivation for using the polytropic equation of state comes from the fact
that our primary purpose of this paper is to explore the tidal disruption
condition for a relatively low value of $\beta<1$ and the formation of shocks
by the tidal compression does not play any role. We here focus only on the
process of tidal disruption and subsequent short-term evolution of the tidally
disrupted material. After the tidal disruption, the fluid is highly elongated
and during the long-term evolution of the fluid elements with different
specific energy and angular momentum, they collide and shocks are likely to be
formed. For such a phase, the shock heating will play an important role. Our
plan is to follow this phase by solving the energy equation with a more
general equation of state.

###  II.3 Initial condition

First, we describe the formulation employed in this paper for computing the
initial data in which white dwarfs are approximately in an equilibrium state
in their comoving frame. From Eq. (28), we have

$\displaystyle\rho u^{\mu}\nabla_{\mu}(hu_{i})+\nabla_{i}P=0,$ (30)

where $h$ is the specific enthalpy defined by $h:=1+\varepsilon+P/\rho$. To
derive Eq. (30), we used Eq. (27).

For the isentropic fluid, the first law of thermodynamics is written as

$\rho dh=dP,$ (31)

where $dQ$ denotes the variation of a quantity $Q$ in the fluid rest frame. In
the polytropic equations of state employed in this work, we obtain the
relation

$\displaystyle h=\int\frac{dP}{\rho}~{}~{}{\rm and}~{}~{}\ln
h=\int\frac{dP}{\rho h}.$ (32)

In this situation, Eq. (30) is rewritten to

$\displaystyle u^{\mu}\nabla_{\mu}(hu_{i})+\partial_{i}h=0.$ (33)

Then, we define $k^{\mu}:=u^{\mu}/u^{t}$. Using this quantity, Eq. (33) is
written to

$u^{t}\mathscr{L}_{k}(hu_{i})-u^{t}hu_{\mu}\nabla_{i}k^{\mu}+\partial_{i}h=0,$
(34)

where $\mathscr{L}_{k}$ denotes the Lie derivative with respect to $k^{\mu}$.
The second term of Eq. (34) is written as

$\displaystyle
u^{t}hu_{\mu}\nabla_{i}k^{\mu}=u^{t}hu_{\mu}\nabla_{i}(u^{\mu}/u^{t})=h\partial_{i}\ln
u^{t},$ (35)

where we used $u^{\mu}u_{\mu}=-1$. Thus, Eq. (34) is written to

$\mathscr{L}_{k}(hu_{i})+\partial_{i}(h/u^{t})=0.$ (36)

We consider an initial condition for a system composed of a star of mass
$M_{*}$ and radius $R_{*}$, for which the center is located on the $x$-axis,
around a massive black hole of mass $M_{\rm BH}\gg M_{*}$ and $M_{\rm BH}\gg
R_{*}$ which is located at a coordinate origin. We assume that the star
predominantly moves toward the $y$-direction with the identical specific
momentum. Thus we set $v^{i}:=u^{i}/u^{t}=-\beta^{i}_{0}+V^{i}$ where
$V^{i}=V\delta^{i}_{~{}y}$ with $V$ being a constant to be determined. Here
the term of $\beta^{i}_{0}$ is added to simplify the iteration process for
computing quasiequilibrium states. Then, $u^{t}$ is calculated from

$u^{t}=\left[\alpha^{2}-\gamma_{ij}(v^{i}+\beta^{i})(v^{j}+\beta^{j})\right]^{-1/2}.$
(37)

In the present context, $\mathscr{L}_{k}(hu_{i})$ can be assumed to be zero
for $i=y$ and $z$, because the star has momentarily translation invariance for
the motion toward the $y$\- and $z$-directions. By contrast, with respect to
the $x$-direction, the star receives the force from the massive black hole.
Since the radius of the star, $R_{*}$, is much smaller than the orbital
separation, $x_{0}$, and $x_{0}$ is larger than the black-hole radius of $\sim
M_{\rm BH}$, $\mathscr{L}_{k}(hu_{i})$ for the $x$–direction can be
approximated by $\partial_{i}[A(x-x_{0})]$ where we take $A$ to be a constant,
which should be approximately written as $\sim-M_{\rm BH}/x_{0}^{2}$ for
$x_{0}>0$. Then, Eq. (33) is integrated to give

$A(x-x_{0})+\frac{h}{u^{t}}=C,$ (38)

where $C$ is an integration constant. We note that Eq. (38) is not an exact
first integral of the Euler equation but can be considered as an approximate
one for obtaining an initial condition in which the star is in an approximate
equilibrium state.

For computing initial conditions, we assume the line elements of the form

$ds^{2}=-(\alpha^{2}-\beta_{k}\beta^{k})dt^{2}+2\beta_{k}dtdx^{k}+\psi^{4}\delta_{ij}dx^{i}dx^{j},$
(39)

where $\psi$ is the conformal factor. Using the Isenberg-Wilson-Mathews
formalism [43, 44], the basic equations are written as

$\displaystyle\Delta\psi$ $\displaystyle=$ $\displaystyle-2\pi\rho_{\rm
H}\psi^{5}-\frac{\psi^{-7}}{8}\hat{A}_{ij}\hat{A}^{ij},$ (40)
$\displaystyle\Delta(\alpha\psi)$ $\displaystyle=$ $\displaystyle
2\pi\alpha\psi^{5}(\rho_{\rm
H}+2S)+\frac{7}{8}\alpha\psi^{-7}\hat{A}_{ij}\hat{A}^{ij},~{}~{}~{}$ (41)
$\displaystyle\partial_{i}\hat{A}^{i}_{~{}j}$ $\displaystyle=$ $\displaystyle
8\pi J_{j}\psi^{6},$ (42)

where

$\displaystyle\rho_{\rm H}$ $\displaystyle=$ $\displaystyle\rho h(\alpha
u^{t})^{2}-P,$ (43) $\displaystyle S$ $\displaystyle=$ $\displaystyle\rho
h[(\alpha u^{t})^{2}-1]+3P,$ (44) $\displaystyle J_{i}$ $\displaystyle=$
$\displaystyle\rho h\alpha u^{t}u_{i},$ (45)

and $\Delta$ is the flat Laplacian. $\hat{A}_{ij}$ is defined from the
extrinsic curvature, $\tilde{A}_{ij}$, by
$\hat{A}_{ij}=\psi^{6}\tilde{A}_{ij}$ and $K$ is set to be zero. Using the CTS
decomposition [45, 46] with trumpet-puncture

$\displaystyle\hat{A}^{ij}=\hat{A}^{ij}_{0}+\frac{\psi^{6}}{2\alpha}\left(\mathbf{L}\beta_{\rm
s}\right)^{ij},$ (46)

Equation (42) is rewritten as

$\displaystyle\partial_{j}\partial^{j}\beta^{i}_{\rm
s}+\frac{1}{3}\delta^{ij}\partial_{j}\partial_{k}\beta^{k}_{\rm s}$
$\displaystyle=$ $\displaystyle 16\pi\alpha\psi^{4}J^{i}$ (47)
$\displaystyle+$ $\displaystyle\left(\mathbf{L}\beta_{\rm
s}\right)^{ij}\partial_{j}\ln\left(\alpha\psi^{-6}\right),$

where $\left(\mathbf{L}\beta_{\rm
s}\right)^{ij}=\left(\delta^{ik}\partial_{k}\beta^{j}_{\rm
s}+\delta^{jk}\partial_{k}\beta^{i}_{\rm
s}-\frac{2}{3}\delta^{ij}\partial_{k}\beta^{k}_{\rm s}\right)$. Note that
although there are some works in constructing binary black holes initial data
with trumpet-puncture [47, 48, 49], this is, to our knowledge, the first
attempt combining the CTS decomposition and puncture method with the limit
(trumpet) hypersurface in constructing quasi-equilibrium initial data in
nonvacuum spacetime. We assume that the contribution to the extrinsic
curvature from the black hole is negligible because the orbital momentum of
the black hole is negligible in this problem, and thus, we set the black hole
at rest (however, it is straightforward to take into account the small black-
hole motion [50] in our formalism.)

For a solution of the initial data, we have to determine the free parameters,
$A$, $C$, and $V$. In the polytropic equation of state, we can consider
$\kappa$ as well as the central density $\rho_{c}$ as free parameters. In the
following, we first consider that $V$ and rest mass of the star are input
parameters and $A$, $C$, and $\kappa$ are parameters to be determined during
the iteration process in numerical computation. Our method to adjust $\kappa$
to a desired value will be described later.

To determine these three parameters we need three conditions, for which we
choose the following relations. First, we fix the location of the surface of
white dwarfs along the $x$-axis as $x=x_{1}$ (referred to as point 1) and
$x=x_{2}$ (point 2). Typically, we choose $x_{1}+x_{2}=2x_{0}$. At the
surface, $h=1$, and thus, Eq. (38) gives

$A(x_{1}-x_{0})+\frac{1}{u^{t}_{1}}=A(x_{2}-x_{0})+\frac{1}{u^{t}_{2}}=C,$
(48)

where $u^{t}_{1}$ and $u^{t}_{2}$ are the values of $u^{t}$ at points 1 and 2.
In addition, we fix the rest mass of the star which is defined by

$m_{*}=\int d^{3}x\rho\alpha\psi^{6}u^{t},$ (49)

where $m_{*}$ is approximately equal to the gravitational mass $M_{*}$ because
the star is only weakly self-gravitating.

Using the condition (48), the values of $C$ and $A$ are determined, and
subsequently, $h$ is determined from Eq. (38). In the polytropic equation of
state, the rest-mass density is written as

$\rho=\left(\frac{(h-1)(\Gamma-1)}{\kappa\Gamma}\right)^{1/(\Gamma-1)},$ (50)

and thus, from Eq. (49), $\kappa$ is determined for given values of $m_{*}$
and $x_{2}-x_{1}$. Once these free parameters are determined, the rest-mass
density are obtained from Eq. (50).

For realistic setting, we have to obtain the desired values of the mass of the
star and the value of $\kappa$. The value of $\kappa$ is controlled by varying
the stellar diameter $x_{2}-x_{1}$ for a given value of $m_{*}$.

To take into account the effect of the black-hole gravity, we employ the
puncture formulation by setting

$\displaystyle\psi$ $\displaystyle=$ $\displaystyle\psi_{0}+\phi,$ (51)
$\displaystyle\alpha\psi$ $\displaystyle=$
$\displaystyle\alpha_{0}\psi_{0}+X,$ (52) $\displaystyle\beta^{k}$
$\displaystyle=$ $\displaystyle\beta^{k}_{0}+\beta^{k}_{\rm s},$ (53)
$\displaystyle\hat{A}_{ij}$ $\displaystyle=$
$\displaystyle\hat{A}_{ij}^{0}+\hat{A}^{\rm s}_{ij},$ (54)

where $\psi_{0}$, $\alpha_{0}$, $\beta^{k}_{0}$, and $\hat{A}_{ij}^{0}$ denote
the solutions of vacuum Einstein’s equation shown already in Sec. II.1. Then
we numerically solve the equations for $\phi$, $X$, $\beta^{k}_{\rm s}$, and
$\hat{A}^{\rm s}_{ij}$ from Eqs. (40), (41), (47), and (46). The initial data
is prepared using the octree-mg code [51], an open source multigrid library
with an octree adaptive-mesh refinement (AMR) grid, which we modified to
support a fourth-order finite-difference elliptic solver.

##  III Numerical simulation

###  III.1 Setup

The simulation is performed using an AMR algorithm with the equatorial
symmetry imposed on the $z=0$ (equatorial) plane using the SACRA-TD code (for
SACRA see Refs. [41, 42]). We prepare two sets of finer domains, one of which
comoves with a white dwarf and the other of which is located around the center
and covers the massive black hole. Because the radius of the white dwarf,
$R_{*}$, is smaller than the black-hole horizon radius $\sim M_{\rm BH}$, we
need to prepare more domains for resolving the white dwarf. In addition to
these domains, we prepare coarser domains that contain both the finer domains
in their inside. All the domains are covered by $(2N+1,2N+1,N+1)$ grid points
for $(x,y,z)$ with $N$ being an even number.

Specifically, each domain is labeled by $i$ which runs as $0,1,2,\cdots,i_{\rm
fix},\cdots,i_{\rm BH},\cdots,i_{\rm max}$. The grid resolution for the
domains with $i_{\rm fix}\leq i\leq i_{\rm BH}$ is identical with that with
$i_{\rm BH}+1\leq i\leq 2i_{\rm BH}-i_{\rm fix}+1(<i_{\rm max})$,
respectively. For $0\leq i\leq i_{\rm BH}$, the center of the domain is
located at the origin, at which a black hole is present. Strictly speaking,
the black hole moves due to the backreaction against the motion of the
companion star, but this motion is tiny because of the condition $M_{\rm
BH}\gg M_{*}$. For these domains, the $i$th level covers a half cubic region
of $[-L_{i}:L_{i}]\times[-L_{i}:L_{i}]\times[0:L_{i}]$ where $L_{i}=N\Delta
x_{i}$, $\Delta x_{i}$ is the grid spacing for the $i$th level, and the grid
spacing for each level is determined by $\Delta x_{i+1}=\Delta x_{i}/2$
($i=0,1,2,\cdots,i_{\mathrm{BH}}-1$ and $i=i_{\mathrm{BH}}+1,\cdots,i_{\rm
max}-1$) with $\Delta x_{i_{\mathrm{BH}}+1}=\Delta x_{i_{\rm fix}}$ and
$L_{i_{\mathrm{BH}}}\sim 0.8M_{\rm BH}$.

For the moving domains that cover the white dwarf, the center is chosen to
approximately agree with the location of the density maximum. In the present
context, the local density maximum is approximately located along a geodesic
around the supermassive black hole. The size of the finest domain with
$i=i_{\rm max}$, $L_{\rm max}$, is chosen so that it is $1.3$–$1.5R_{*}$. We
check the convergence of two different models with three grid resolutions as
illustrated in Fig. 1. Higher resolution is used for model M8V17 to measure
the spin up of the white dwarf more accurately (see Sec. III.2). We obtain
good convergence for both models, and thus, we employ $N=60$ as the standard
resolution in this paper.

Figure 1: Maximum density as a function of time for model M7V165 (top) with $N=40$, $60$, $70$, and model M8V17 (bottom) with $N=60$, $82$, $102$. We find that a fair convergence is obtained with $N=60$.  Table 1: Models considered in this paper and the fate (last column). M7V16a and M7V16b correspond to the models with $R_{*}=8.5\times 10^{3}$ and $7.0\times 10^{3}$ km, respectively. For other models, $R_{*}\approx 10^{4}(M_{*}/0.7M_{\odot})^{-1/3}$ km. $r_{p}$ and $r_{p,A}$ are periastron radius in the present coordinates and the Schwarzschild coordinates, respectively. TD and OC denote tidal disruption and appreciable oscillation of white dwarfs, and NN denotes that no appreciable tidal effect is found.  ID |  $V$ | $M_{*}(M_{\odot})$ | $r_{p}/M_{\rm BH}$ | $r_{p,A}/M_{\rm BH}$ | $J/M_{\rm BH}$ |  $\beta$ |  Fate
---|---|---|---|---|---|---|---
M6V16 | 0.160 | 0.6 | 4.401 | 5.456 | 3.775 | 0.72 | TD
M7V16 | 0.160 | 0.7 | 4.401 | 5.456 | 3.775 | 0.65 | TD
M7V16a | 0.160 | 0.7 | 4.401 | 5.456 | 3.775 | 0.55 | TD
M7V16b | 0.160 | 0.7 | 4.401 | 5.456 | 3.775 | 0.45 | TD/OC
M8V16 | 0.160 | 0.8 | 4.401 | 5.456 | 3.775 | 0.59 | TD
M7V165 | 0.165 | 0.7 | 5.770 | 6.813 | 3.897 | 0.52 | TD
M8V165 | 0.165 | 0.8 | 5.770 | 6.813 | 3.897 | 0.47 | TD/OC
M6V17 | 0.170 | 0.6 | 7.030 | 8.065 | 4.019 | 0.49 | TD/OC
M7V17 | 0.170 | 0.7 | 7.030 | 8.065 | 4.019 | 0.44 | OC
M8V17 | 0.170 | 0.8 | 7.030 | 8.065 | 4.019 | 0.40 | OC
M6V175 | 0.175 | 0.6 | 8.317 | 9.346 | 4.142 | 0.42 | OC
M7V18 | 0.180 | 0.7 | 9.681 | 10.707 | 4.265 | 0.33 | NN

###  III.2 Numerical results

Figure 2: Geodesics for $V=0.160$, 0.165, 0.170, 0.175 and 0.180 in the
coordinates of $g_{\mu\nu}^{0}$. Those only for one orbital period started
from $x=20M_{\rm BH}$ and $y=z=0$ are plotted. The filled circle at the center
represents the black hole with the coordinate radius of its event horizon
$r\approx 0.78M_{\rm BH}$.

In the present paper we focus on the case that the black-hole mass is $M_{\rm
BH}=10^{5}M_{\odot}$, the white-dwarf mass is $M_{*}=0.6$, $0.7$, and
$0.8M_{\odot}$. For the polytropic equation of state, the stellar radius,
$R_{*}$, is proportional to $M_{*}^{(\Gamma-2)/(3\Gamma-4)}$ for a fixed value
of $\kappa$. Thus, for $\Gamma=5/3$, the stellar radius depends only weakly on
the stellar mass. In the present case we basically choose the value of
$\kappa$ so that $R_{*}\approx 1.0\times 10^{4}(M_{*}/0.7M_{\odot})^{-1/3}$
km. For $M_{*}=0.7M_{\odot}$ and $V=0.160$, we also prepare two additional
cases where $\kappa$ is chosen such that $R_{*}=8.5\times 10^{3}$ km and
$R_{*}=7.0\times 10^{3}$ km.

The initial separation is set to be $x_{0}=20M_{\rm BH}$ (it is $\approx
21.01M_{\rm BH}$ in the Schwarzschild coordinates), and $V$ is chosen to be
$0.160$, $0.165$, $0.170$, $0.175$, and $0.180$ (see Table 1). The
corresponding specific angular momentum of the white dwarf is $J\approx
3.7748$, $3.8968$, $4.0192$, $4.142$, and $4.2653M_{\rm BH}$, and the
resulting periastron radius is $r_{p}/M_{\rm BH}(r_{p,A}/M_{\rm
BH})=4.401(5.456)$, 5.770 (6.813), 7.030 (8.065), 8.317 (9.346), and 9.681
(10.707) where in the parenthesis the values in the Schwarzschild coordinate,
i.e., areal radius (hereafter denoted by $r_{p,A}$), are described. In Fig. 2,
we plot the geodesics only for one orbital period for $V=0.160$, 0.165, 0.170,
and 0.180.

With these settings, the white dwarf has an elliptic orbit around the black
hole with the periastron at $r_{p}\approx(4.4$–$10)M_{\rm BH}$, and thus, the
eccentricity is approximately defined by $e=(x_{0}-r_{p})/(x_{0}+r_{p})$ is
$\approx 1/3$–$2/3$. Here, $x_{0}(=20M_{\rm BH})$ and $r_{p}$ are defined in
the radial coordinates of the metric of $g^{0}_{\mu\nu}$, and thus, the values
of $e$ slightly change if we define it in the areal coordinate (Schwarzschild
radial coordinate).

For the models mentioned above, the value of $\beta$ is in the range between
0.33 and 0.72 and estimated by

$\displaystyle\beta$ $\displaystyle\approx$ $\displaystyle
0.59\left(\frac{R_{*}}{10^{4}\,{\rm
km}}\right)\left(\frac{M_{*}}{0.7M_{\odot}}\right)^{-1/3}$ (55)
$\displaystyle~{}~{}~{}~{}\times\left(\frac{r_{p,A}}{6M_{\rm
BH}}\right)^{-1}\left(\frac{M_{\rm BH}}{10^{5}M_{\odot}}\right)^{-2/3},$

where the areal radius $r_{p,A}$ is used for the definition of $\beta$ in this
section. For $V=0.160$ and $0.165$ with $M_{*}=0.6$–$0.8M_{\odot}$, we find
$0.50\leq\beta\leq 0.7$, and thus, the white dwarf is expected to be strongly
perturbed by the black-hole tidal field for $M_{*}=0.6$–$0.8M_{\odot}$. By
contrast, for $V=0.180$, $\beta<0.35$ with $M_{*}=0.7M_{\odot}$, and thus, the
tidal force of the black hole is likely to be too weak to perturb the white
dwarf.

For $V=0.170$, $\beta\approx 0.49$, $0.44$, and $0.40$ with $M_{*}=0.6$,
$0.7$, and $0.8M_{\odot}$ respectively. In these cases, tidal disruption is
not very likely to take place but the tidal force from the black hole should
induce the stellar oscillation on the white dwarf. Because for $\Gamma=5/3$,
the stellar radius depends only weakly on the stellar mass, the presence or
absence of the tidal disruption is likely to depend primarily on the value of
$V$ (or the specific angular momentum of the white dwarfs) in the present
setting. In the following, we will show that our code can reproduce all these
expected phenomena.

Figure 3: The maximum density as a function of time for $V=0.160$, 0.165,
0.170, and 0.180 with $M_{*}=0.7M_{\odot}$. The maximum density is normalized
by the initial value denoted by $\rho_{0}$.

Figure 3 plots the evolution of the maximum density for $V=0.160$, 0.165,
0.170, and 0.180 with $M_{*}=0.7M_{\odot}$. We note that for $M_{\rm
BH}=10^{5}M_{\odot}$, the orbital period for these parameters are in the range
from $\approx 220$ s for $V=0.160$ to $\approx 250$ s for $V=0.180$. The
figure shows the results expected in the previous paragraphs: For M7V16 and
M7V165, the white dwarfs are tidally disrupted while approaching the black
hole irrespective of the white-dwarf mass. For M7V17 ($\beta=0.44$), the white
dwarf is perturbed by the black hole near the periastron but it is not tidally
disrupted. After the close encounter, the white dwarf is in an oscillating
state due to the instantaneous tidal force received from the black hole. By
contrast, for M7V18, the maximum density is approximately preserved to be
constant, suggesting no disruption occurs and the tidal effect is negligible.
Note that such tidal field may still perturb the white dwarf and produce
detectable electromagnetic or gravitational-wave signal if a sufficient
amplitude of oscillation is induced.

In Fig. 3, the results of M7V16 ($\beta=0.65$), M7V16a ($\beta=0.55$) and
M7V16b ($\beta=0.45$) are also compared. As expected, for the first two
models, the white dwarfs are tidally disrupted, while for the most compact
white dwarf, the tidal disruption does not occur although it is perturbed
significantly by the black-hole tidal force. This illustrates that the $\beta$
parameter is a good indicator for assessing whether tidal disruption takes
place or not irrespective of the white-dwarf radius.

Figure 4: The maximum density as a function of time for the cases that stellar
oscillation occurs $(0.5\gtrsim\beta\gtrsim 0.4)$. The red, blue, and green
curves show the results with $M_{*}=0.6M_{\odot}$, $0.7M_{\odot}$ and
$0.8M_{\odot}$, respectively. The maximum density is normalized by the initial
value denoted by $\rho_{0}$.

Figure 4 shows the evolution of the maximum density when stellar oscillation
is induced. For M6V17 ($\beta=0.49$), the white dwarf is significantly
elongated by the tidal force from the black hole; the central density is
decreased to less than 50% of the original value after passing through the
periastron. Associated with the tidal effect, the mass is lost from the white
dwarf. However, with the increase of the orbital radius, the central density
increases again, resulting in a less massive white dwarf. This is also the
case for M8V165 ($\beta=0.47$) and M7V16b ($\beta=0.45$). These results
indicate that the critical value of $\beta$ for the tidal disruption is $\sim
0.50$ and the threshold value for exciting a high-amplitude oscillation is
$\beta\sim 0.45$. Figure 4 also shows that even for $0.40\lesssim\beta\lesssim
0.45$ an appreciable oscillation is excited by the tidal force.

Figure 5: A summary for the fate of the white dwarfs in the plane of $M_{*}$
and $\beta$. TD and OC denote that tidal disruption and appreciable
oscillation of the white dwarfs are observed after the close encounter of the
white dwarfs with the black hole. NN denotes that no appreciable tidal effect
is observed.

In Fig. 5 and Table 1, we summarize the fates of white dwarfs as a result of
the tidal interaction. It is found that for $\beta\gtrsim 0.5$ tidal
disruption takes place and for $\beta\gtrsim 0.4$, the white dwarfs are
perturbed appreciably by the black-hole tidal field. All these results agree
approximately with the expectation from the previous studies.

Figure 6: The density profiles of the tidally-disrupted white dwarf for the
model $V=0.160$ and $M_{*}=0.7M_{\odot}$ (M7V16). The units of the length
scale for the density plots are $GM_{\rm BH}/c^{2}\approx 1.48\times 10^{5}$
km. The solid and dashed curves show the time evolution for the location of
the maximum density and the elliptic orbit shown in Fig. 2 for $V=0.160$
(i.e., geodesic). The length scale of $x$ and $y$ axes is shown in units of
$M_{\rm BH}$.

For M7V16, tidal disruption takes place but only a small fraction of the white
dwarf matter falls into the black hole because the fluid elements have
specific angular momentum large enough to escape capturing by the black hole.
Most of the tidally disrupted matter approximately maintains the original
elliptic orbit (see Fig. 6) although the matter has an elongated profile. To
clarify the eventual matter distribution around the black hole, we will need
to follow the matter motion for more than 10 orbits. This topic is one of our
major research targets in the future.

For $0.4\lesssim\beta\lesssim 0.5$, the white dwarf will be continuously
perturbed by the black-hole tidal force whenever it passes through the
periastron. In addition the angular momentum is transported during the tidal
interaction, and it will lead to the transport of the orbital angular momentum
to the white dwarf resulting in a spin-up of it. According to a perturbation
study for the stellar encounter, the energy deposition during the tidal
interaction in one orbit is written approximately as [52]

$\displaystyle\Delta E_{\rm tid}$ $\displaystyle=$ $\displaystyle f_{\rm
tid}\left(\frac{M_{*}^{2}}{R_{*}}\right)\left(\frac{M_{\rm
BH}}{M_{*}}\right)^{2}\left(\frac{R_{*}}{r_{p,A}}\right)^{6}$ (56)
$\displaystyle=$ $\displaystyle f_{\rm
tid}\left(\frac{M_{*}^{2}}{R_{*}}\right)\beta^{6}.$

where $f_{\rm tid}$ is a factor of $O(0.1)$, which depends on $\beta$ and the
equation of state. Associated with the energy deposition near the periastron,
the angular momentum deposition is also deposited. In one orbit it is
approximately estimated by $\Delta J_{\rm spin}\approx\Delta E_{\rm
tid}/\Omega_{p}$ [53] where $\Omega_{p}=\sqrt{M_{\rm BH}/r_{p,A}^{3}}$, and
thus,

$\displaystyle\Delta J_{\rm spin}$ $\displaystyle=$ $\displaystyle f_{\rm
spin}M_{*}\sqrt{M_{*}R_{*}}\left(\frac{M_{\rm
BH}}{M_{*}}\right)^{3/2}\left(\frac{R_{*}}{r_{p,A}}\right)^{9/2}$ (57)
$\displaystyle=$ $\displaystyle f_{\rm
spin}M_{*}\sqrt{M_{*}R_{*}}\beta^{9/2},$

where $f_{\rm spin}$ is a coefficient of the same order of the magnitude of
$f_{\rm tid}$.

Figure 7: The rescaled change in angular momentum $\Delta
J_{spin}/\left(M_{*}\sqrt{M_{*}R_{*}}\beta^{9/2}\right)$ as a function of time
for the stellar oscillation scenario (M8V17). This agrees with the analytic
expression Eq. (57) if $f_{\rm spin}\sim 0.1$–0.3.

Because the maximum spin angular momentum of the star is approximately written
as $M_{*}\sqrt{M_{*}R_{*}}$, we find that $\Delta J_{\rm spin}$ can be more
than 0.1% of the maximum spin if a white dwarf passes through a close orbit
with $\beta\gtrsim 0.4$. We approximate the orbit angular momentum $J_{\rm
orbit}$ and the spin angular momentum $J_{\rm spin}$ of the white dwarf as

$\displaystyle J_{\rm orbit}=$ $\displaystyle
M_{h}\left(\left<x\right>\left<u_{y}\right>-\left<y\right>\left<u_{x}\right>\right),$
(58) $\displaystyle J_{\rm spin}=$ $\displaystyle\int d^{3}x\psi^{6}\rho
h\left[\left(x-\left<x\right>\right)\left(u_{y}-\left<u_{y}\right>\right)\right.$
(59)
$\displaystyle\left.-\left(y-\left<y\right>\right)\left(u_{x}-\left<u_{x}\right>\right)\right],$

where $M_{h}:=\int d^{3}x\psi^{6}\rho h$ and the volume average of quantity
$q$ is defined as $\left<q\right>:=\frac{1}{M_{h}}\int d^{3}x\psi^{6}\rho hq$.
In such decomposition, the sum of orbital and spin angular momentum equals to
the total angular momentum of the white dwarfs. We analyzed the spin angular
momentum gain of the white dwarfs for M8V17, and we indeed find $\Delta J_{\rm
spin}/(M_{*}\sqrt{M_{*}R_{*}}\beta^{9/2})\approx 0.1$–0.3 as shown in Fig. 7.
Note that the spin up of white dwarf $\Delta J_{\rm spin}$ is about $10^{-6}$
of the total angular momentum, and hence, it is not easy to determine $\Delta
J_{\rm spin}$ accurately. Although we cannot achieve a good convergence in
$\Delta J_{\rm spin}$, we are able to obtain a noticeable rise in $J_{\rm
spin}$ during the close encounter, which suggests $f_{\rm spin}\sim 0.1$–0.3,
which is consistent with the above analytic result.

For close orbits, the tidal angular-momentum transport can dominate over the
orbital angular momentum loss by gravitational-radiation reaction. Assuming
that gravitational waves are most efficiently emitted near the periastron at
which we may approximate the orbit to be circular, the angular momentum
dissipation by gravitational waves in one orbit can be written as [54]

$\displaystyle\Delta J_{\rm GW}\approx\frac{64\pi}{5}\frac{M_{\rm
BH}^{2}M_{*}^{2}}{r_{p,A}^{2}}\left(1+\frac{7e^{2}}{8}\right),$ (60)

where $e$ denotes the eccentricity. Thus, the ratio of $\Delta J_{\rm tid}$ to
$\Delta J_{\rm GW}$ is written as

$\displaystyle\frac{\Delta J_{\rm spin}}{\Delta J_{\rm GW}}$
$\displaystyle\approx$ $\displaystyle 23f_{\rm
spin}\left(\frac{r_{p,A}}{4M_{\rm BH}}\right)^{-5/2}\left(\frac{M_{\rm
BH}}{10^{5}M_{\odot}}\right)^{-3}\left(\frac{M_{*}}{0.7M_{\odot}}\right)^{-2}$
(61) $\displaystyle\times\left(\frac{R_{*}}{10^{4}\,{\rm
km}}\right)^{5}\left(1+\frac{7e^{2}}{8}\right)^{-1}.$

Thus it is larger than unity for $r_{p,A}\lesssim 7M_{\rm BH}/c^{2}$,
$R_{*}\approx 10^{4}$ km, $M_{\rm BH}=10^{5}M_{\odot}$, $M_{*}=0.7M_{\odot}$,
and $f_{\rm spin}=0.2$. This is also the case for the ratio of $\Delta E_{\rm
tid}/\Delta E_{\rm GW}$ where $\Delta E_{\rm GW}$ is the energy dissipated by
gravitational waves in one orbit. Thus, near the tidal disruption orbit, the
orbital evolution would be primarily determined not by the gravitational-wave
emission but by the tidal effect. To clarify the eventual fate of such a white
dwarf, we obviously need a long-term accurate simulation. Such a topic is one
of our future targets.

We note that both $\Delta J_{\rm spin}$ and $\Delta J_{\rm GW}$ are much
smaller than the orbital angular momentum of order $M_{*}\sqrt{M_{\rm
BH}r_{p,A}}$. Thus, the cumulative effect of the tidal angular momentum
transport plays an important role just prior to the tidal disruption. By
repeated tidal interaction, the spin angular velocity of the white dwarfs is
likely to be enhanced up to $\sim M_{\rm
BH}^{1/2}/r_{p}^{3/2}=\beta^{3/2}M_{*}^{1/2}/R_{*}^{3/2}$. In addition, the
stellar oscillation for which the oscillation energy is comparable to or
larger than the rotational kinetic energy should be excited. As a result, mass
loss could be induced, resulting in the increase of the stellar radius and
enhancing the importance of tidal interaction. In this type of the system, the
tidal disruption is unlikely to take place by one strong impact by the black-
hole tidal force but likely to do as a result of a secular increase of the
stellar radius (see, e.g., Refs. [55, 9] for related studies).

##  IV Summary

We reported a new numerical-relativity code which enables us to explore tidal
disruption of white dwarfs by a relatively low-mass supermassive black hole.
As a first step toward more detailed future studies, we paid attention to the
condition for tidal disruption of white dwarfs with typical mass range in
elliptic orbits by a nonspinning supermassive black hole. We showed that our
code is capable of determining the condition for the tidal disruption. As
expected from previous general relativistic works (e.g., Refs. [30, 23]), the
tidal disruption takes place for $\beta\gtrsim 0.5$ and an appreciable
oscillation of the white dwarfs are induced by the black-hole tidal effect for
$\beta\gtrsim 0.4$ for orbits close to the black hole in the $\Gamma=5/3$
polytropic equation of state. The critical value for the onset of the tidal
disruption is smaller than that obtained by Newtonian analysis. For white
dwarfs with $M_{*}=0.6M_{\odot}$ and $R_{*}=1.2\times 10^{4}$ km, $\beta$ can
be larger than 0.4 even for $M_{\rm BH}\approx 4\times 10^{5}M_{\odot}$ if the
periastron radius is $r_{p,A}=4M_{\rm BH}$. Our result indicates that in such
systems with a relatively low-mass (but not intermediate-mass) supermassive
black hole for which gravitational waves in the late inspiral phase can be
detected by LISA [4], tidal disruption can occur for typical white dwarfs. For
spinning black holes with the dimensionless spin parameter of $\gtrsim 0.9$,
$r_{p,A}$ can be smaller than $\sim 1.7M_{\rm BH}$ [56]. For such black holes,
tidal disruption of typical-mass white dwarfs may occur even for $M_{\rm
BH}\approx 10^{6}M_{\odot}$. Investigation of this possibility is a future
issue.

There are several issues to be explored. The first one is to extend our
implementation for spinning black holes. Since no analytic solution is known
for the spacetime of spinning black holes on the limit hypersurface, we need
to develop a method to provide $g_{\mu\nu}^{0}$ for employing the formulation
introduced in this paper. One straightforward way to prepare such data is just
to numerically perform a simulation for a spinning black hole (in vacuum)
until the hypersurface reaches the limit hypersurface as a first step, and
then, the obtained data are saved and used in the subsequent simulations with
white dwarfs. A more subtle issue along this line is to prepare the initial
condition. For nonspinning black holes, we can assume that the conformal
flatness of the three metric, and as a result, the initial-value equations are
composed only of elliptic-type equations with the flat Laplacian. For the
spinning black holes, the basic equations are composed of elliptic-type
equations of complicated Laplacian, and hence, the numerical computation could
be more demanding, although in principle it would be still possible to obtain
an initial condition. We plan to explore this strategy in the subsequent work.

For modeling realistic white dwarfs it is necessary to implement a realistic
equation of state. If we assume that the temperature of the white dwarfs is
sufficiently low and the pressure is dominated by that of degenerate
electrons, it is straightforward to implement this.

More challenging issue is to follow the hydrodynamics of tidally disrupted
white-dwarf matter for a long term. After the tidal disruption, the matter of
the white dwarf is likely to move around the black hole for many orbits.
During such orbits, the matter collides each other, and eventually, a compact
disk will be formed after the circularization. Such disks are likely to be hot
due to the shock heating, and thus, it can be a source of electromagnetic
counterparts of the tidal disruption. In the presence of magnetic fields,
magnetorotational instability [57] occurs in the disk, and the magnetic fields
will be amplified. If the amplified magnetic field eventually penetrates the
black hole and if the black hole is appreciably spinning, a jet may be
launched through the Blandford-Znajek effect [58]. After the amplification of
the magnetic fields, a turbulent state will be developed in the disk and mass
ejection could occur by the effective viscosity or magneto-centrifugal force
[59]. The ejecta may be a source of electromagnetic signals. One long-term
issue is to investigate such scenarios by general relativistic
magnetohydrodynamics.

###### Acknowledgements.

We thank Kenta Hotokezaka for helpful discussion. This work was in part
supported by Grant-in-Aid for Scientific Research (Grant No. JP20H00158) of
Japanese MEXT/JSPS. Numerical computations were in part performed on Sakura
clusters at Max Planck Computing and Data Facility.

## References

  * van Velzen  _et al._ [2020] S. van Velzen, T. W. S. Holoien, F. Onori, T. Hung, and I. Arcavi, Space Sci. Rev. 216, 124 (2020), arXiv:2008.05461 [astro-ph.HE] .
  * Alexander  _et al._ [2020] K. D. Alexander, S. van Velzen, A. Horesh, and B. A. Zauderer, Space Sci. Rev. 216, 81 (2020), arXiv:2006.01159 [astro-ph.HE] .
  * Saxton  _et al._ [2021] R. Saxton, S. Komossa, K. Auchettl, and P. G. Jonker, Correction to: X-Ray Properties of TDEs, Space Science Reviews, Volume 217, Issue 1, article id.18 (2021), arXiv:2103.15442 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), arXiv e-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Vick  _et al._ [2017] M. Vick, D. Lai, and J. Fuller, MNRAS 468, 2296 (2017).
  * Metzger  _et al._ [2022] B. D. Metzger, N. C. Stone, and S. Gilbaum, Astrophys. J. 926, 101 (2022), arXiv:2107.13015 [astro-ph.HE] .
  * King [2022] A. King, MNRAS 515, 4344 (2022), arXiv:2206.04698 [astro-ph.GA] .
  * Lu and Quataert [2022] W. Lu and E. Quataert, (2022), arXiv:2210.08023 [astro-ph.HE] .
  * Linial and Sari [2022] I. Linial and R. Sari, (2022), arXiv:2211.09851 [astro-ph.HE] .
  * Rosswog  _et al._ [2008] S. Rosswog, E. Ramirez-Ruiz, and W. R. Hix, ApJ 679, 1385 (2008), arXiv:0712.2513 [astro-ph] .
  * Rosswog  _et al._ [2009] S. Rosswog, E. Ramirez-Ruiz, and W. R. Hix, ApJ 695, 404 (2009), arXiv:0808.2143 [astro-ph] .
  * MacLeod  _et al._ [2014] M. MacLeod, J. Goldstein, E. Ramirez-Ruiz, J. Guillochon, and J. Samsing, ApJ 794, 9 (2014), arXiv:1405.1426 [astro-ph.HE] .
  * Hills [1975] J. G. Hills, Nature 254, 295 (1975).
  * Laguna  _et al._ [1993] P. Laguna, W. A. Miller, and W. H. Zurek, ApJ 404, 678 (1993).
  * Khokhlov  _et al._ [1993] A. Khokhlov, I. D. Novikov, and C. J. Pethick, ApJ 418, 163 (1993).
  * Diener  _et al._ [1997] P. Diener, V. P. Frolov, A. M. Khokhlov, I. D. Novikov, and C. J. Pethick, ApJ 479, 164 (1997).
  * Rossi  _et al._ [2021] E. M. Rossi, N. C. Stone, J. A. P. Law-Smith, M. Macleod, G. Lodato, J. L. Dai, and I. Mandel, Space Sci. Rev. 217, 40 (2021), arXiv:2005.12528 [astro-ph.HE] .
  * Maguire  _et al._ [2020] K. Maguire, M. Eracleous, P. G. Jonker, M. MacLeod, and S. Rosswog, Space Sci. Rev. 216, 39 (2020), arXiv:2004.00146 [astro-ph.HE] .
  * Gafton and Rosswog [2019] E. Gafton and S. Rosswog, MNRAS 487, 4790 (2019), arXiv:1903.09147 [astro-ph.HE] .
  * Ryu  _et al._ [2020a] T. Ryu, J. Krolik, T. Piran, and S. C. Noble, ApJ 904, 98 (2020a), arXiv:2001.03501 [astro-ph.HE] .
  * Ryu  _et al._ [2020b] T. Ryu, J. Krolik, T. Piran, and S. C. Noble, ApJ 904, 99 (2020b), arXiv:2001.03502 [astro-ph.HE] .
  * Ryu  _et al._ [2020c] T. Ryu, J. Krolik, T. Piran, and S. C. Noble, ApJ 904, 100 (2020c), arXiv:2001.03503 [astro-ph.HE] .
  * Ryu  _et al._ [2020d] T. Ryu, J. Krolik, T. Piran, and S. C. Noble, ApJ 904, 101 (2020d), arXiv:2001.03504 [astro-ph.HE] .
  * Ryu  _et al._ [2022] T. Ryu, J. Krolik, and T. Piran, arXiv e-prints , arXiv:2211.00059 (2022), arXiv:2211.00059 [astro-ph.HE] .
  * Guillochon and Ramirez-Ruiz [2013] J. Guillochon and E. Ramirez-Ruiz, ApJ 767, 25 (2013), arXiv:1206.2350 [astro-ph.HE] .
  * Mainetti  _et al._ [2017] D. Mainetti, A. Lupi, S. Campana, M. Colpi, E. R. Coughlin, J. Guillochon, and E. Ramirez-Ruiz, A&A 600, A124 (2017), arXiv:1702.07730 [astro-ph.HE] .
  * Luminet and Marck [1985] J. P. Luminet and J. A. Marck, MNRAS 212, 57 (1985).
  * Luminet and Carter [1986] J. P. Luminet and B. Carter, ApJS 61, 219 (1986).
  * Fishbone [1973] L. G. Fishbone, ApJ 185, 43 (1973).
  * Ishii  _et al._ [2005] M. Ishii, M. Shibata, and Y. Mino, Phys. Rev. D 71, 044017 (2005), arXiv:gr-qc/0501084 [astro-ph] .
  * Marck [1983] J. A. Marck, Proceedings of the Royal Society of London Series A 385, 431 (1983).
  * Haas  _et al._ [2012] R. Haas, R. V. Shcherbakov, T. Bode, and P. Laguna, ApJ 749, 117 (2012), arXiv:1201.4389 [astro-ph.HE] .
  * East and Pretorius [2013] W. E. East and F. Pretorius, Phys. Rev. D 87, 101502 (2013), arXiv:1303.1540 [gr-qc] .
  * East [2014] W. E. East, ApJ 795, 135 (2014), arXiv:1408.1695 [gr-qc] .
  * Shibata and Nakamura [1995] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995).
  * Baumgarte and Shapiro [1999] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 024007 (1999), gr-qc/9810065 .
  * Alcubierre  _et al._ [2003] M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 67, 084023 (2003), arXiv:gr-qc/0206072 .
  * Baumgarte [2011] T. W. Baumgarte, Classical and Quantum Gravity 28, 215003 (2011), arXiv:1108.3550 [gr-qc] .
  * Estabrook  _et al._ [1973] F. Estabrook, H. Wahlquist, S. Christensen, B. Dewitt, L. Smarr, and E. Tsiang, Phys. Rev. D 7, 2814 (1973).
  * Baumgarte and Naculich [2007] T. W. Baumgarte and S. G. Naculich, Phys. Rev. D 75, 067502 (2007), arXiv:gr-qc/0701037 [gr-qc] .
  * Yamamoto  _et al._ [2008] T. Yamamoto, M. Shibata, and K. Taniguchi, Phys. Rev. D 78, 064054 (2008), arXiv:0806.4007 [gr-qc] .
  * Kiuchi  _et al._ [2017] K. Kiuchi, K. Kawaguchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D 96, 084060 (2017), arXiv:1708.08926 [astro-ph.HE] .
  * Isenberg [2008] J. A. Isenberg, International Journal of Modern Physics D 17, 265 (2008), arXiv:gr-qc/0702113 [gr-qc] .
  * Wilson and Mathews [1995] J. R. Wilson and G. J. Mathews, Phys. Rev. Lett. 75, 4161 (1995).
  * York [1999] J. York, James W., Phys. Rev. Lett. 82, 1350 (1999), arXiv:gr-qc/9810051 [gr-qc] .
  * Baumgarte  _et al._ [2007] T. W. Baumgarte, N. Ó. Murchadha, and H. P. Pfeiffer, Phys. Rev. D 75, 044009 (2007), arXiv:gr-qc/0610120 [gr-qc] .
  * Immerman and Baumgarte [2009] J. D. Immerman and T. W. Baumgarte, Phys. Rev. D 80, 061501 (2009), arXiv:0908.0337 [gr-qc] .
  * Dietrich and Brügmann [2014] T. Dietrich and B. Brügmann, Phys. Rev. D 89, 024014 (2014), arXiv:1309.3087 [gr-qc] .
  * Clemente and de Oliveira [2017] P. C. M. Clemente and H. P. de Oliveira, Phys. Rev. D 96, 024035 (2017), arXiv:1703.09131 [gr-qc] .
  * Slinker  _et al._ [2018] K. Slinker, C. R. Evans, and M. Hannam, Phys. Rev. D 98, 044014 (2018), arXiv:1806.08364 [gr-qc] .
  * Teunissen and Keppens [2019] J. Teunissen and R. Keppens, Computer Physics Communications 245, 106866 (2019), arXiv:1901.11370 [physics.comp-ph] .
  * Press and Teukolsky [1977] W. H. Press and S. A. Teukolsky, ApJ 213, 183 (1977).
  * Kumar and Quataert [1998] P. Kumar and E. J. Quataert, ApJ 493, 412 (1998).
  * Peters [1964] P. C. Peters, Physical Review 136, 1224 (1964).
  * Dai and Blandford [2013] L. Dai and R. D. Blandford, Mon. Not. Roy. Astron. Soc. 434, 2948 (2013), arXiv:1110.2614 [astro-ph.HE] .
  * Bardeen  _et al._ [1972] J. M. Bardeen, W. H. Press, and S. A. Teukolsky, ApJ 178, 347 (1972).
  * Balbus and Hawley [1998] S. A. Balbus and J. F. Hawley, Reviews of Modern Physics 70, 1 (1998).
  * Blandford and Znajek [1977] R. D. Blandford and R. L. Znajek, MNRAS 179, 433 (1977).
  * Blandford and Payne [1982] R. D. Blandford and D. G. Payne, MNRAS 199, 883 (1982).