# Black hole-neutron star mergers with massive neutron stars in numerical
relativity

Shichuan Chen  Department of Astronomy, School of Physics, Peking University,
Beijing 100871, China  Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Am Mühlenberg 1, Postdam-Golm 14476,
Germany  Luohan Wang  School of Physics, Peking University, Beijing 100871,
China  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Am Mühlenberg 1, Postdam-Golm 14476, Germany  Kota Hayashi  Max
Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, Postdam-Golm 14476, Germany  Center for Gravitational Physics
and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502, Japan  Kyohei Kawaguchi  Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Postdam-
Golm 14476, Germany  Center for Gravitational Physics and Quantum Information,
Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502,
Japan  Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Kenta Kiuchi  Max Planck
Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg
1, Postdam-Golm 14476, Germany  Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto
606-8502, Japan  Masaru Shibata  Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Am Mühlenberg 1, Postdam-Golm 14476,
Germany  Center for Gravitational Physics and Quantum Information, Yukawa
Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

###### Abstract

We study the merger of black hole-neutron star (BH-NS) binaries in numerical
relativity, focusing on the properties of the remnant disk and the ejecta,
varying the mass of compactness of the NS and the mass and spin of the BH. We
find that within the precision of our numerical simulations, the remnant disk
mass and ejecta mass normalized by the NS baryon mass ($\hat{M}_{\rm{rem}}$
and $\hat{M}_{\rm{eje}}$, respectively), and the cutoff frequency
$f_{\rm{cut}}$ normalized by the initial total gravitational mass of the
system at infinite separation approximately agree among the models with the
same NS compactness $C_{\rm{NS}}=M_{\rm{NS}}/R_{\rm{NS}}$, mass ratio
$Q=M_{\rm{BH}}/M_{\rm{NS}}$, and dimensionless BH spin $\chi_{\rm{BH}}$
irrespective of the NS mass $M_{\rm{NS}}$ in the range of
$1.092$–$1.691\,M_{\odot}$. This result shows that the merger outcome depends
sensitively on $Q$, $\chi_{\rm BH}$, and $C_{\rm{NS}}$ but only weekly on
$M_{\rm{NS}}$. This justifies the approach of studying the dependence of NS
tidal disruptions on the NS compactness by fixing the NS mass but changing the
EOS. We further perform simulations with massive NSs of
$M_{\rm{NS}}=1.8M_{\odot}$, and compare our results of $\hat{M}_{\rm{rem}}$
and $\hat{M}_{\rm{eje}}$ with those given by existing fitting formulas to test
their robustness for more compact NSs. We find that the fitting formulas
obtained in the previous studies are accurate within the numerical errors
assumed, while our results also suggest that further improvement is possible
by systematically performing more precise numerical simulations.

##  I Introduction

The first gravitational wave (GW) detection for GW150914 from a binary black-
hole (BBH) merger heralded the opening of the era of GW astronomy [1].
Subsequently, at the event of GW170817, the first binary neutron-star (BNS)
merger was observed not only through GWs but also through signals of
electromagnetic (EM) counterparts by diverse instruments all over the world
[2, 3, 4, 5, 6]. GW data analysis for this event gave a constraint on the
tidal deformability of the NS which excludes very stiff equations of state
(EOS) [2, 7, 8, 9, 10]. The observation of short gamma-ray burst (GRB), GRB
170817A, and the kilonova, AT2017gfo, indicated that the NS is involved in the
merger, and it is suggested that a remnant formed after a BNS merger is likely
to be the central engine o the GRB [4, 5, 11]. The observation of the kilonova
indicated that an $r$-process nucleosynthesis could occur in the BNS merger.
All these facts suggest that future observations of GWs and EM counterparts
will provide us valuable opportunities to deepen the knowledge of the
mechanism of short GRBs, kilonovae, and $r$-process nucleosyntheses of heavy
elements.

As in the BNS mergers, the BH-NS mergers can also, in principle, generate EM
counterparts of GW such as kilonovae when a significant amount of matter is
ejected in the tidal disruption of the NS. Such EM counterparts can give
information about the NS EOS [12, 8] , the origin of heavy elements produced
through $r$-process nucleosynthesis [13, 14, 15, 16, 17], and cast light on
physics beyond the nuclear saturation density. However, although two events
GW200105 and GW200115 [18] are reported to be the BH-NS merger in 2021, and in
addition, GW190425 (with masses of 1.12-1.68 $M_{\odot}$ and 1.61-2.52
$M_{\odot}$) could possibly be a BH-NS merger [19, 20], no multi-messenger
detection of BH-NS mergers has been confirmed yet. No EM counterpart
observation for the previous events is consistent with theoretical predictions
because the mass ratios of these binaries are as high as 4–5 and the BH spins
are likely to be zero or retrograde, and thus, tidal disruption is unlikely to
occur [21]. However, for interpreting future events, in which tidal disruption
of NSs may occur, it is crucial to prepare theoretical models that quantify
the relation between the properties of the BH-NS merger and the observables.

Numerical simulations of BH-NS mergers play an important role in understanding
the tidal disruption and mass ejection processes. Since 2006 [22, 23, 24, 25,
26], many numerical simulations have been performed for this to study the
dependence of the merger behavior on binary parameters, focusing on the
characteristic results of tidal disruptions [27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42, 43]. Additionally, it is revealed that the
cutoff frequency, one of the notable features of the GW signal in tidal
disruptions, encodes information on the NS EOS, especially when it is stiff
[29, 44, 31, 45]. More importantly, previous studies [32, 38, 46, 33, 47] show
that the restmass of the remnant matter located outside the apparent horizon
and ejecta mass normalized by the NS baryon mass ($\hat{M}_{\rm{rem}}$ and
$\hat{M}_{\rm{eje}}$, respectively) as well as the cutoff frequency of GWs
depend sensitively only on the mass ratio $Q$, the BH dimensionless spin
$\chi_{\rm BH}$, the compactness of NSs $C_{\rm{NS}}$ but not on the NS mass
$M_{\rm{NS}}$. Based on this approximate parameter dependence, previous
studies [38, 47, 33, 46] of BH-NS mergers give some fitting formulas for the
remnant mass $M_{\rm{rem}}$ and the ejecta mass $M_{\rm{eje}}$. These semi-
analytical models for the merger results are valuable because numerical
simulations of BH-NS mergers are resource-intensive. These fitting formulas
for the remnant disk mass and ejecta mass are used to assess whether the EM
counterparts would be present [48, 38], and constrain the binary parameters
after the observation of EM signals [49, 50, 51, 52, 53, 54].

However, the typical range of NS mass in these previous simulations is limited
to $M_{\rm{NS}}\approx 1.2$–$1.5M_{\odot}$ [55, 56], with NS compactness
between 0.12 and 0.19, although more massive NSs at least up to $\sim
2M_{\odot}$ exist in nature [57, 58]. It is not clear whether the fitting
formulas obtained by the previous studies also hold quantitatively for more
massive and compact NSs because they are not tested for more massive NSs.
Given possible GW detections of such massive NSs in BH-NS binaries [19, 59],
quantitative investigation of the BH-NS mergers in such a parameter space is
essential.

In this paper, we study BH-NS mergers in numerical relativity varying the NS
mass and NS compactness for a range wider than in the previous studies. We
first focus on BH-NS binaries with the several fixed combinations of mass
ratio $Q$, BH dimensionless spin $\chi_{\rm BH}$, and NS compactness
$C_{\rm{NS}}$ but with various NS masses $M_{\rm NS}$ in order to clarify how
$\hat{M}_{\rm{rem}}$, $\hat{M}_{\rm{eje}}$, and $f_{\rm{cut}}m_{0}$ (where
$m_{0}=M_{\mathrm{BH}}+M_{\mathrm{NS}}$) are sensitive to the NS mass $M_{\rm
NS}$. We then perform simulations for a large value of the NS mass ($M_{\rm
NS}=1.8\,M_{\odot}$), and examine the accuracy of the previous fitting
formulas of $\hat{M}_{\rm{rem}}$, $\hat{M}_{\rm{eje}}$, and ejecta velocity
$v_{\rm{eje}}$ for the system with a large value of the NS compactness
($C_{\rm NS}\geq 0.19$). We will show that within the range of the numerical
accuracy, $M_{\rm{rem}}$, $M_{\rm{eje}}$, and $f_{\rm{cut}}m_{0}$ can be
driven from the fitting formulas if the values of $\chi_{\rm{BH}}$, $Q$, and
$C_{\rm{NS}}$ are identical: The fitting formulas previously obtained in Refs.
[38, 47, 46] still work well for larger values of $C_{\rm{NS}}$ above $0.19$.

The paper is organized as follows. In Sec. II we briefly summarize numerical
methods used in our study. Numerical results of simulations are presented in
Sec. III, focusing on the remnant disk mass, the ejecta mass, and the
gravitational waveforms. In Secs. III.2 and III.3, we present the results for
the simulations with the fixed NS compactness and those for massive and
compact NSs, respectively. Sec. IV is devoted to a summary and discussion.
Throughout this paper, all the quantities in tables and figures are shown in
units of $c=G=M_{\odot}=1$ unless otherwise stated. Our convention of notation
is summarized in Table 1.

Table 1: Our convention of notation for physically important quantities, geometric variables, and hydrodynamic variables. Symbol | Description
---|---
$M_{\rm{BH}}$ | Gravitational mass of the black hole in isolation
$\chi_{\rm{BH}}$ | Dimensionless spin parameter of the black hole
$M_{\rm{NS}}$ | Gravitational mass of the neutron star in isolation
$R_{\rm{NS}}$ | Circumferential radius of the neutron star in isolation
$M_{\rm{b}}$ | Rest baryon mass of the neutron star in isolation
$C_{\rm{NS}}$ | Compactness of the neutron star: $C_{\rm{NS}}=M_{\rm{NS}}/R_{\rm{NS}}$
$Q$ | Mass ratio of BH-NS: $Q=M_{\rm{BH}}/M_{\rm{NS}}$
$M_{\rm{rem}}$ | Remnant disk mass after the merger
$M_{\rm{eje}}$ | Ejecta mass after the merger
$\hat{M}$ | Mass normalized by $M_{\rm{b}}$
$m_{0}$ | Initial total gravitational mass of the system
| at infinite separation
$\gamma_{ij}$ | Induced metric on the $t$ = const hypersurface
$\alpha$ | Lapse function
$\beta^{i}$ | Shift vector
$\gamma$ | Determinant of $\gamma_{ij}$
$\rho$ | Baryon rest-mass density
$u^{\mu}$ | Four velocity of the fluid
$P$ | Pressure
$\epsilon$ | Specific internal energy
$\rho_{*}$ | Conserved baryon rest mass density: $\rho\alpha\sqrt{\gamma}u^{t}$
$\Gamma$ | Adiabatic index
$h$ | Specific enthalpy: $1+\epsilon+P/\rho$
$\hat{e}$ | Specific energy: $h\alpha u^{t}-P/\rho\alpha u^{t}$

##  II Numerical methods

This section describes the numerical methods used in our study. Sec. II.1
summarizes our numerical simulation methods, Sec. II.2 describes the NS EOS,
Sec. II.3 presents the diagnostics of simulations, Sec. II.4 presents the
fitting formulas given from previous studies, and Sec. II.5 presents the
models employed in this study.

###  II.1 Simulation methods

We carry out numerical simulations using the SACRA-MPI code [60], which uses
an adaptive-mesh-refinement (AMR) algorithm [27] and MPI/OpenMP hybrid
parallelization to speed up the computation [60]. SACRA solves the Einstein
equation in a moving puncture version [61, 62, 63] of the Baumgarte-Shapiro-
Shibata-Nakamura (BSSN) formulation [64, 65], incorporating a Z4c constraint-
propagation prescription locally [66]. Together with the Einstein equation, we
solve hydrodynamics equations using the Harten-Lax-van Leer contact (HLLC)
solver [67, 68, 69] in this paper. More details of the formulation, the gauge
conditions, and the numerical scheme are described in Refs. [31, 45, 32, 42].
For actual computation, we employ the public spectral code FUKA (Frankfurt
University Kadath) [70] to generate initial data. We do not take into account
magnetohydrodynamics and neutrino-transfer effects, since we focus on the
evolution of the system up to 30 ms after the merger for which these effects
are not very important.

The grid structure of SACRA is summarized as follows. Computational domains
are composed of nested equidistant Cartesian grids, and each grid has
$(2N,2N,N)$ points in $(x,y,z)$ directions. We employ a cell-centered grid
structure and the $x$-coordinate at $j$th grid point is given by
$x_{j}=(j+1/2)\Delta x$, where $\Delta x$ is the grid resolution. The
equatorial-plane symmetry is imposed on the $z=0$ plane. We adopt $N=82$ as a
fiducial value, with which the NS radius is covered by $\approx 66$ points in
the finest grid. We also perform simulations with $N=62$ and $N=102$ for
selected models to evaluate the numerical error and check the convergence of
the numerical results (see Appendix A). The physical quantities in tables and
figures, unless specified, are taken from results with $N=82$. In our study,
the inspiral motion before the merger is followed at least for 5 orbits.

We prepare 10 refinement levels for the AMR computational domains. 6 coarser
domains ($l_{c}$) cover both BH and NS with the origins fixed at the center of
the mass of the binary system. Two sets of 4 finer domains ($l_{f}$) comove
with either BH or NS, covering the region of their vicinity. Namely, we
prepare $2l_{f}+l_{c}$ computational grid domains spanning $l_{c}+l_{f}$
refinement levels. Starting from the coarsest level as $l=0$, the $l$th level
has a grid spacing, $\Delta x_{l}=L/(2^{l}N)$, and at the finest level,
$\Delta x=L/(2^{l_{f}+l_{c}-1}N)$, where $L$ is the size of the computational
domain, which covers $[-L:L]\times[-L:L]\times[0:L]$ for $x$-$y$-$z$.

###  II.2 Equation of state

Since the lifetime of NS binaries is typically much longer than the cooling
time of NSs, NSs in the late inspiral stage are believed to be cold enough to
be modeled by zero-temperature EOS [71]. Thus, for modeling the NSs prior to
the merger, we use the zero-temperature EOS in our simulations. In this paper,
we model such EOSs by piecewise polytropes [72, 73, 74], which are written in
the form:

$P(\rho)=\kappa_{i}\rho^{\Gamma_{i}}\ (\rho_{i}\leq\rho\leq\rho_{i+1}),$ (1)

where $\kappa_{i}$ are constant and other quantities are described in Table 1.
We perform simulations for models with both two pieces and four pieces, with
$i\in\left\\{0,1\right\\}$ and $i\in\left\\{0,1,2,3\right\\}$ respectively.
$\rho_{0}=0$, and for $\rho\leq\rho_{1}$, we adopt $\Gamma_{0}=1.35692395$ and
$\kappa_{0}/c^{2}=3.99873692\times 10^{-8}\
\rm{g^{1-\Gamma_{0}}cm^{3\Gamma_{0}-3}}$ [72]. The boundary condition at
$\rho_{i}$ is
$\kappa_{i}\rho_{i}^{\Gamma_{i}}=\kappa_{i+1}\rho_{i}^{\Gamma_{i+1}}$, which
requires the continuity of pressure at the interface between $i^{\rm th}$ and
$(i+1)^{\rm th}$ region.

For different EOS models, we take different fiducial pressures $P_{\rm{fid}}$
at the fiducial density $\rho_{\rm{fid}}=10^{14.7}\ \rm{g\,cm^{-3}}$, and
thus, $\kappa_{2}$ is defined by
$\kappa_{2}=P_{\rm{fid}}/\rho_{\rm{fid}}^{\Gamma_{2}}$. For piecewise
polytropes with two pieces, we take
$\rho_{2}=+\infty,\Gamma_{1}=\Gamma_{2}=3.0$, and for piecewise polytropes
with four pieces, we take $\rho_{2}=\rho_{\rm{fid}}$, $\rho_{3}=10^{15}~{}{\rm
g\,cm^{-3}}$, and $\rho_{4}=\infty$. Other parameters for different EOS models
and NS initial properties related to the EOS are listed in Table 2.

In numerical simulations, we add the thermal part of the EOS to the zero-
temperature part described above. Our implementation for it is the same as
that described in Refs. [75] with an adiabatic index of the thermal part
$\Gamma_{\rm th}=1.8$.

Table 2: Piecewise polytropic EOS Models employed in the present simulations. Model EOS name, $P_{\rm fid}$: the pressure at the fiducial density $\rho_{\rm fid}=10^{14.7}\rm{g\,cm^{-3}}$, $\Gamma_{i}$ ($i=1$–3), and the maximum gravitational mass of spherical NSs for the given EOS. The first four EOSs are composed only of two peices. EOS | $\log_{10}P_{\rm fid}[\rm{dyne/cm^{2}}]$ | $\Gamma_{1}$ | $\Gamma_{2}$ | $\Gamma_{3}$ | $M_{\rm{max}}[M_{\odot}]$
---|---|---|---|---|---
15H | 34.700 | 3.000 | 3.000 | 3.000 | 2.53
H | 34.500 | 3.000 | 3.000 | 3.000 | 2.25
HB | 34.400 | 3.000 | 3.000 | 3.000 | 2.12
B | 34.300 | 3.000 | 3.000 | 3.000 | 2.00
H4 | 34.669 | 2.909 | 2.246 | 2.144 | 2.03
APR4 | 34.269 | 2.830 | 3.445 | 3.348 | 2.20

###  II.3 Diagnostics

####  II.3.1 Remnant disk and ejecta

After the merger, the fate of the matter originating from the NS can be
divided into three types. The matter directly falls into the BH, forms a
remnant disk, and becomes unbound from the system, i.e. ejecta. Evaluating the
properties of the disk and the ejecta is essential for discussing resultant
models of EM counterparts. Here we describe the method to evaluate these
quantities. In the following, we will define the time of merger $t_{\rm
merger}$ as the time when $0.01M_{\odot}$ of the NS matter falls into the
apparent horizon.

At each time slice, we evaluate the remnant disk by the rest mass outside the
apparent horizon, i.e. the remnant disk mass with the integral

$M_{\rm>AH}:=\int_{r>r_{\rm AH}}\rho_{*}\,d^{3}x$ (2)

where $r_{\rm{AH}}=r_{\rm{AH}}(\theta,\phi)$ denotes the coordinate radius of
the apparent horizon and $\rho_{*}$ is denfined in Table 1.

The unbound matter of the system, ejecta, is defined as the matter satisfying
$-u_{t}>1$, because in this paper we focus only on the dynamical ejecta, for
which the thermal effect is minor. The mass of ejecta at each time slice is
defined by the integral

$M_{\rm eje}:=\int_{-u_{t}>1,r>r_{\rm AH}}\rho_{*}\,d^{3}x.$ (3)

Assuming that the thermal energy of the ejecta is much smaller than the
kinetic energy, the kinetic energy of the ejecta may be defined by

$T_{\rm eje}:=\int_{-u_{t}>1,r>r_{\rm AH}}\rho_{*}(-u_{t})\,d^{3}x.$ (4)

Then the average velocity of the ejecta is estimated from the kinetic energy
$T_{\rm{eje}}$ and the ejecta mass $M_{\rm eje}$ as

$v_{\rm eje}:=\sqrt{\frac{2T_{\rm eje}}{M_{\rm eje}}}.$ (5)

Following the previous studies [45, 32, 33, 38, 47, 46], we evaluate all these
ejecta quantities at 10 ms after the onset of the merger.

####  II.3.2 Black hole parameters

During the inspiral or post-merger stage, the parameters of BHs are estimated
from the quantities of the apparent horizons. If we assume that the spacetime
near the apparent horizon is stationary, the equatorial circumferential radius
$C_{\rm{e}}$ and the area of the apparent horizon $A_{\rm{AH}}$ can be
approximated as [29]

$C_{\rm e}=4\pi M_{\rm BH},$ (6)

and

$A_{\rm AH}=8\pi M_{\rm BH}^{2}\left(1+\sqrt{1-\chi_{\rm BH}^{2}}\right).$ (7)

Then BH mass $M_{\rm BH}$ and dimensionless spin parameter
$\chi_{\mathrm{BH}}$ can be evaluated from Eqs. (6) and (7).

From the comparison of different estimations of the spin parameter, it has
been found that the systematic error with this method to estimate the spin is
much less than 0.01 [31, 45, 29]. Therefore, the spin parameter used in the
paper for developing fitting formulas is accurate enough for our study.

####  II.3.3 Gravitational waves

To derive gravitational waveforms, we extract the Weyl scalar $\Psi_{4}$ at
the coordinate radius of $D=480M_{\odot}$ from the coordinate origin and
extrapolate them to the null infinity with a method from the BH perturbation
theory [76]. The gravitational waveforms $h_{\rm GW}$ are obtained from the
time integration of the $l=|m|=2$ spherical harmonics modes of $\Psi_{4}$. All
the waveforms in this paper are shown for the observer along the $z$-axis. To
reduce the nonphysical low-frequency components, we do not perform the time
integration directly but adopt the method proposed by Reisswig and Pollney
[77]. The retard time $t_{\rm ret}$ is approximately obtained by

$t_{\rm ret}=t-D-2m_{0}\ln(D/m_{0}).$ (8)

The GW spectrum is obtained as the sum of the Fourier components of the two
polarizations of $l=|m|=2$ modes:

$\tilde{h}(f)=\sqrt{\frac{|\tilde{h}_{+}(f)|^{2}+|\tilde{h}_{\times}(f)|^{2}}{2}},$
(9)

where

$\tilde{h}_{+,\times}(f)=\int e^{2\pi ift}h_{+,\times}(t)\ dt.$ (10)

We use the normalized amplitude $f|\tilde{h}(f)|D/m_{0}$ as a function of a
dimensionless frequency $fm_{0}$ to show the GW spectrum in the following
plots.

###  II.4 Fitting formulas

Former studies give fitting formulas for $\hat{M}_{\rm{rem}}$ [38, 47] and
$\hat{M}_{\rm{eje}}$ [33, 46]. In these fitting formulas, $\hat{M}_{\rm{rem}}$
and $\hat{M}_{\rm{eje}}$ are determined by the mass ratio $Q$, the
dimensionless BH spin $\chi_{\rm{BH}}$, and the NS compactness $C_{\rm{NS}}$
except for the one derived in Ref. [33], which also depends on the value of
$M_{\rm NS}/M_{\rm b}$. In Sec. III.2 we will demonstrate that this
parametrization is sufficient for estimating $\hat{M}_{\rm{rem}}$ and
$\hat{M}_{\rm{eje}}$.

1

In Ref. [38], the fitting formula for the remnant mass (referred to as
rem_2012) is determined as

$\displaystyle\frac{M^{\rm{rem}}_{\rm{fit},2012}}{M_{\rm{b}}}$
$\displaystyle=\alpha(3Q)^{1/3}(1-2C_{\rm{NS}})-\beta\frac{R_{\rm{ISCO}}}{R_{\rm{NS}}},$
(11) $\displaystyle=\alpha(3Q)^{1/3}(1-2C_{\rm{NS}})-\beta
r_{\rm{ISCO}}(\chi_{\rm{BH}})QC_{\rm{NS}},$ (12)

where $\alpha$ and $\beta$ are determined by minimizing the $\chi$-square (as
defined in Eq. (23)) of the fitting formula, $R_{\rm{ISCO}}$ is the radius of
the innermost stable circular orbit (ISCO) of the remnant BH with
$r_{\mathrm{ISCO}}=R_{\mathrm{ISCO}}/M_{\mathrm{BH}}$, and other quantities
are defined in Table 1. The ranges of the initial NS and BH parameters covered
by the numerical simulations are $Q=3$–$7$, $\chi_{\rm{BH}}=0$–$0.9$,
$C_{\rm{NS}}=0.13$–$0.18$, and the best-fit model gives $\alpha\approx 0.288$
and $\beta\approx 0.148$. Here $R_{\rm{ISCO}}$ is given by

$\displaystyle Z_{1}=1+(1-\chi_{\rm{BH}})^{1/3}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times[(1+\chi_{\rm{BH}})^{1/3}+(1-\chi_{\rm{BH}})^{1/3}]$
(13) $\displaystyle Z_{2}=\sqrt{3\chi_{\rm{BH}}^{2}+Z_{1}^{2}}$ (14)
$\displaystyle\frac{R_{\rm{ISCO}}}{M_{\rm{BH}}}=3+Z_{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\rm{sign}(\chi_{\rm{BH}})\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})},$
(15)

as in Ref. [78].

2

Ref. [47] gives a modified fitting formula for the remnant mass (referred to
as rem_2018)

$\frac{M^{\rm{rem}}_{\rm{fit},2018}}{M_{\rm{b}}}=\\\
\left[\rm{Max}\left(\alpha\frac{1-2C_{\rm{NS}}}{\eta^{1/3}}-\beta
r_{\rm{ISCO}}\frac{C_{\rm{NS}}}{\eta}+\gamma,0\right)\right]^{\delta},$ (16)

where $\eta$ denotes the symmetric mass ratio $\eta=Q/(1+Q)^{2}$. The best
fitting gives $\alpha=0.406,\beta=0.139,\gamma=0.255,\delta=1.761$ and the
parameter space for the simulations was $Q=1$–$7$,
$\chi_{\rm{BH}}=-0.5$–$0.9$, $C_{\rm{NS}}=0.13$–$0.182$ and $M_{\rm{rem}}\leq
0.3M_{\rm{b}}$.

3

Ref. [33] gives the fitting formula for the ejecta mass (referred to as
eje_2016) as

$\displaystyle\frac{M^{\rm{eje}}_{\rm{fit},2016}}{M_{\rm{b}}}$
$\displaystyle=\text{Max}\Big{[}a_{1}Q^{n_{1}}\frac{1-2C_{\rm{NS}}}{C_{\rm{NS}}}-a_{2}Q^{n_{2}}r_{\rm{ISCO}}(\chi_{\rm{eff}})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+a_{3}(1-\frac{M_{\rm{NS}}}{M_{\rm{b}}})+a_{4},\
0\Big{]}$ (17) $\displaystyle\chi_{\rm{eff}}$
$\displaystyle=\chi_{\rm{BH}}\cos{i_{\rm{tilt}}}$ (18)

In our simulations, the angle between the BH spin and the orbital angular
momentum is $i_{\rm{tilt}}=0$. The best-fit gives $a_{1}=4.464\times 10^{-2},\
a_{2}=2.269\times 10^{-2},\ a_{3}=2.431,\ a_{4}=-0.4159,\ n_{1}=0.2497,\
n_{2}=1.352$, and the parameter space for simulations was $Q=3$–$7$,
$\chi_{\rm{BH}}=0$–$0.75$, and $C_{\rm{NS}}=0.138$–$0.18$ by fixing
$M_{\rm{NS}}=1.35M_{\odot}$ and changing EOS.

4

Ref. [46] gives a modified fitting formula for the ejecta mass (referred to as
eje_2020)

$\displaystyle\frac{M^{\rm{eje}}_{\rm{fit},2020}}{M_{\rm{b}}}=a_{1}Q^{n_{1}}\frac{1-2C_{\rm{NS}}}{C_{\rm{NS}}}-a_{2}Q^{n_{2}}r_{\rm{ISCO}}+a_{4}.$
(19)

The fitting results are $a_{1}=0.007116,\ a_{2}=0.001436,\ a_{4}=-0.02762,\
n_{1}=0.8636,\ n_{2}=1.6840$, and the parameter space for simulations was the
same as that of Ref. [33].

5

Ref. [33] also gives a fitting formula for the average ejecta velocity (cf,
Eq. (5)) as a simple linear function of $Q$:

$v^{\rm{eje}}_{\rm{fit,2016}}=(0.01533Q+0.1907)c.$ (20)

The definition of average ejecta velocity used to fit
$v^{\rm{eje}}_{\rm{fit,2016}}$ is given by Eq. (31) in Appendix B, which is
not the same as $v_{\rm{eje}}$ given in Eq. (5).

###  II.5 Models

Table 3: Parameters of initial data and grid structure for all the models
studied in this paper. The model name contains the EOS employed, the mass
ratio ($Q$), and the dimensionless spin of the black hole $\chi_{\rm BH}$.
$M_{\rm{NS}},M_{\rm{NS}}^{\rm{b}},R_{\rm{NS}},C_{\rm{NS}}$, and $M_{\rm BH,0}$
are the initial NS mass, NS restmass, NS radius, NS compactness, and BH mass,
respectively. $m_{0}\Omega_{0}$ is the dimensionless initial orbital angular
velocity of the system. $L$ is the box size of the simulation, $\Delta x$ is
the grid spacing at the finest level, and $R_{\rm diam}/\Delta x$ is the grid
number within the semimajor diameter of the NS. $N$ denotes the grid
resolution of the simulation. The table contains three parts. The first part
shows the parameters for the fixed value of the NS mass as
$M_{\rm{NS}}=1.8M_{\odot}$, the second and third parts show the parameters for
the fixed values of $C_{\rm{NS}}=0.182$ and $0.147$, respectively. $\Delta x$
and $R_{\rm{diam}}/\Delta x$ are calculated with the grid resolution of
$N=82$.

Model | EOS | $M_{\rm NS}[M_{\odot}]$ | $M_{\rm{NS}}^{\rm{b}}[M_{\odot}]$ |  $R_{\rm{NS}}$[km] | $C_{\rm NS}$ | $Q$ | $M_{\rm BH,0}[M_{\odot}]$ | $\chi_{\rm BH}$ | $m_{0}\Omega_{0}$ | $\Delta x[\rm m]$ | $R_{\rm diam}/\Delta x$ | $L[\rm km]$ | $N$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
15H-Q3a5M18 | 15H | 1.800 | 2.022 | 13.70 | 0.194 | 3.0 | 5.400 | 0.50 | 0.0377 | 204 | 134 | 8553 | 62, 82
H-Q3a5M18 | H | 1.800 | 2.056 | 12.14 | 0.219 | 3.0 | 5.400 | 0.50 | 0.0377 | 180 | 135 | 7565 | 62, 82
15H-Q3a75M18 | 15H | 1.800 | 2.022 | 13.70 | 0.194 | 3.0 | 5.400 | 0.75 | 0.0375 | 204 | 134 | 8553 | 62, 82, 102
H-Q3a75M18 | H | 1.800 | 2.056 | 12.14 | 0.219 | 3.0 | 5.400 | 0.75 | 0.0375 | 180 | 135 | 7565 | 62, 82
B-Q3a75M18 | B | 1.800 | 2.103 | 10.55 | 0.252 | 3.0 | 5.400 | 0.75 | 0.0375 | 163 | 129 | 6842 | 62, 82
15H-Q5a75M18 | 15H | 1.800 | 2.022 | 13.70 | 0.194 | 5.0 | 9.000 | 0.75 | 0.0445 | 203 | 134 | 8523 | 62, 82
H-Q5a75M18 | H | 1.800 | 2.056 | 12.14 | 0.219 | 5.0 | 9.000 | 0.75 | 0.0445 | 180 | 135 | 7565 | 62, 82
15H-Q3a5M1691 | 15H | 1.691 | 1.884 | 13.72 | 0.182 | 3.0 | 5.073 | 0.50 | 0.0369 | 204 | 135 | 8553 | 62, 82
HB-Q3a5M1428 | HB | 1.428 | 1.590 | 11.59 | 0.182 | 3.0 | 4.284 | 0.50 | 0.0353 | 171 | 136 | 7182 | 62, 82
15H-Q3a75M1691 | 15H | 1.691 | 1.884 | 13.72 | 0.182 | 3.0 | 5.073 | 0.75 | 0.0367 | 204 | 135 | 8553 | 62, 82
HB-Q3a75M1428 | HB | 1.428 | 1.590 | 11.59 | 0.182 | 3.0 | 4.284 | 0.75 | 0.0351 | 171 | 136 | 7182 | 62, 82, 102
APR4-Q3a75M1366 | APR4 | 1.366 | 1.522 | 11.08 | 0.182 | 3.0 | 4.098 | 0.75 | 0.0358 | 164 | 135 | 6901 | 62, 82
H4-Q3a75M1651 | H4 | 1.651 | 1.836 | 13.40 | 0.182 | 3.0 | 4.950 | 0.75 | 0.0356 | 199 | 135 | 8347 | 62, 82
15H-Q5a75M1691 | 15H | 1.691 | 1.884 | 13.72 | 0.182 | 5.0 | 8.455 | 0.75 | 0.0411 | 204 | 135 | 8553 | 62, 82, 102
HB-Q5a75M1428 | HB | 1.428 | 1.590 | 11.59 | 0.182 | 5.0 | 7.140 | 0.75 | 0.0404 | 171 | 136 | 7182 | 62, 82, 102
H-Q3a5M1220 | H | 1.220 | 1.327 | 12.26 | 0.147 | 3.0 | 3.660 | 0.50 | 0.0351 | 182 | 135 | 7639 | 62, 82
B-Q3a5M1092 | B | 1.092 | 1.187 | 10.97 | 0.147 | 3.0 | 3.276 | 0.50 | 0.0357 | 163 | 135 | 6842 | 62, 82
15H-Q3a75M1363 | 15H | 1.363 | 1.484 | 13.69 | 0.147 | 3.0 | 4.089 | 0.75 | 0.0349 | 203 | 135 | 8523 | 62, 82
H-Q3a75M1220 | H | 1.220 | 1.327 | 12.26 | 0.147 | 3.0 | 3.660 | 0.75 | 0.0349 | 182 | 135 | 7639 | 62, 82, 102
B-Q3a75M1092 | B | 1.092 | 1.187 | 10.97 | 0.147 | 3.0 | 3.276 | 0.75 | 0.0356 | 163 | 135 | 6842 | 62, 82
APR4-Q3a75M1094 | APR4 | 1.094 | 1.190 | 10.99 | 0.147 | 3.0 | 3.282 | 0.75 | 0.0356 | 163 | 135 | 6842 | 62, 82
H-Q5a75M1220 | H | 1.220 | 1.327 | 12.26 | 0.147 | 5.0 | 6.100 | 0.75 | 0.0449 | 182 | 135 | 7639 | 62, 82
B-Q5a75M1092 | B | 1.092 | 1.187 | 10.97 | 0.147 | 5.0 | 5.460 | 0.75 | 0.0451 | 163 | 135 | 6842 | 62, 82

Table 3 lists BH-NS binary models and initial parameters used in our study.
The models are labeled by EOS-$Q\chi_{\rm{BH}}M_{\rm{NS}}$, i.e. 15H-Q3a75M18
means the model with NS EOS “15H”, mass ratio $Q=3$, NS mass
$M_{\rm{NS}}=1.8M_{\odot}$, and dimensionless BH spin $\chi_{\rm{BH}}=0.75$
aligned with the orbital angular momentum.

In Sec. III.2, to verify our hypothesis that the normalized disk mass and
ejecta mass only depend on $Q$, $\chi_{\rm{BH}}$, and $C_{\rm{NS}}$, we
perform simulations by varying the NS mass with fixing these three parameters,
and explore whether the differences of results among different models are
within the range of numerical accuracy. To compare our results with previous
ones, we choose three configurations of $Q$ and $\chi_{\rm{BH}}$ as
$(Q,\chi_{\rm{BH}})=\ (3,0.5),\ (3,0.75)$ and $(5,0.75)$. We study the models
fixing the NS compactness with $C_{\rm{NS}}=0.182$, but employing EOS models
with two-pieces 15H, HB, and four-pieces EOS models H4 and APR4 [79, 80, 81]
to vary the NS mass. Since the ejecta mass is small in the
$C_{\mathrm{NS}}=0.182$ cases, we also perform simulations with fixed
compactness $C_{\rm{NS}}=0.147$, with two-pieces EOS models 15H, H, B, and
four-pieces EOS model APR4, for which we expect to get a larger disk and
ejecta masses. In the following, we will refer to the models with fixed
compactness $C_{\rm{NS}}=0.147$ and $C_{\rm{NS}}=0.182$ as C147 and C182,
respectively.

We note that the purpose of studying the models with four-pieces EOS is to
test whether the simplified two-pieces EOSs can be used to accurately describe
the NS tidal disruption in BH-NS mergers and to test how the difference in the
detailed internal structure of the NS can change the results of the
simulations. Two-pieces polytopes are simplified EOSs, which are typically
employed in addressing less massive NSs. It can be proved that for an NS with
mass less than $1.4M_{\odot}$, its central density is relatively low, which
means the EOS at higher density plays a minor role. Hence, two-pieces and
four-pieces EOS should provide approximately the same results for the C147
models. On the other hand, for the C182 models, the NS structure in high-
density regions can be sensitive to the values of $\Gamma_{2}$ and
$\Gamma_{3}$. Therefore, it is necessary to perform simulations with both two-
pieces and four-pieces EOS models to test whether our conclusions can be
applied despite the detailed structure of the NS.

In Sec. III.3, we further perform simulations with a high value of
$M_{\rm{NS}}=1.8M_{\odot}$ while fixing $Q$ and $\chi_{\mathrm{BH}}$ to check
the validity of the fitting formulas with large values of NS compactness. As
we demonstrate in Sec. III.2 that merger dynamics depends on $M_{\rm{NS}}$ and
$R_{\rm{NS}}$ approximately only through the NS compactness ($C_{\rm
NS}=M_{\rm{NS}}/R_{\rm{NS}}$), it is an interesting exploration to see the
dependence of the merger dynamics by varying the compactness of NSs while
fixing $M_{\rm{NS}}$, i.e. by changing the EOS of NSs. We choose EOSs, 15H, H,
and B, with the compactness varying from 0.194 to 0.252.

To evaluate the validity of fitting formulas, we calculate the $\chi$-square
of their fittings, with the estimation of numerical error from both Refs. [33,
38] and our study. Numerical simulations are performed with grid resolutions
$N=62,\ 82$ for all the models, and $N=102$ for some of them (see Appendix A)
to estimate the numerical errors due to the finite grid resolution. We take
$N=82$ as the fiducial resolution of our simulations and other results are
summarized in the Appendix A.

##  III Results

Figure 1: Time evolution of the rest-mass density profile for model
HB-Q3a75M1428 with $N=82$ run at $t-t_{\rm merger}\approx-0.94\,\rm{ms}$ (top
left), $0.07\,\rm{ms}$ (top right), $0.90\,\rm ms$ (bottom left), and $10.01\
\rm ms$ (bottom right). The white circle indicates the apparent horizon. The
white dashed lines show the region of unbound components in which the relation
$-u_{t}\geq 1$ is satisfied.

We present our numerical results in this section, focusing primarily on their
dependence on the NS compactness, the mass ratio, and the dimensionless BH
spin. After we describe an overview of the BH-NS merger process in Sec. III.1,
we summarize the result for the case of a fixed NS compactness in Sec. III.2,
and then, for the case of massive and compact NSs in Sec. III.3. The numerical
results are shown in Tables 4, 5, and 6, respectively.

###  III.1 Overview of the merger

First we briefly describe an overview of the BH-NS merger process in our
numerical simulations. Broadly speaking, we find the same features as for
$M_{\rm NS}\approx 1.4\,M_{\odot}$ studied in the previous works [31, 45, 32,
22, 23, 24, 29, 35, 82, 25, 26, 28, 83, 84, 30, 39, 38, 85, 86, 37, 40, 87,
41, 88, 42, 33, 89, 90, 36, 91, 92]: The BH-NS mergers can have two types of
fate. One is that the NS is tidally disrupted by the BH, and subsequently, an
accretion disk is formed around the remnant BH and a part of the NS matter is
ejected as the dynamical ejecta; the other is that the NS simply plunges into
the BH with no tidal disruption and no mass ejection. In the following, we pay
attention only to the former case because it can provide a lot of information
about the BH-NS systems through the cutoff frequency and EM counterparts
powered by the ejecta.

The condition of tidal disruption can be obtained by comparing the ISCO radius
$R_{\rm ISCO}$ and the mass-shedding radius $R_{\rm dis}$. If $R_{\rm
ISCO}\leq R_{\rm dis}$, the NS can be tidally disrupted before it is swallowed
into the BH. Assuming Newton’s gravity for simplicity, we have [21]

$\frac{R_{\rm{dis}}}{R_{\rm
ISCO}}\propto\frac{1}{C_{\rm{NS}}Q^{2/3}r_{\rm{ISCO}}(\chi_{\rm{BH}})}.$ (21)

This ratio approximately describes the significance of tidal disruption: The
smaller this ratio is, the weaker the tidal disruption is (if it is less than
unity, the tidal disruption is not expected). Thus, the significance of the
tidal disruption appears to depend on $Q$, $\chi_{\mathrm{BH}}$, and
$C_{\rm{NS}}$, and the tidal disruption is weaker for larger values of $Q$ and
$C_{\rm{NS}}$ and smaller value of $\chi_{\rm{BH}}$ since $r_{\rm{ISCO}}$
decreases as $\chi_{\rm{BH}}$ increases. However, the relation of this ratio
may be modified in general relativity, in particular for high-mass NSs for
which the general relativistic effect is enhanced. The dependence of ejecta
mass and the remnant mass on these three parameters as well as on the NS mass
are discussed in detail in Secs. III.2 and III.3.

Figure 1 shows the snapshot of the rest-mass density profile, unbound matter,
and the apparent horizon at selected time slices for model HB-Q5a75M1428. The
separation between the NS and BH decreases due to GW emission, and after it
reaches the tidal disruption radius, the NS is tidally disrupted (see the top
panels in Fig. 1). While most of the disrupted NS matter is swallowed by the
BH, a fraction of the matter will gain angular momentum by gravitational
torque to form the tidal tail which remains outside the BH. The
gravitationally-bounded matter in the tidal tail gradually falls back and
forms a remnant disk around the remnant BH (see the bottom panels in Fig. 1),
while the matter in the outer part of the tidal tail becomes dynamical ejecta.

Figure 2 shows the time evolution for the rest mass of the matter located
outside the BH for three models, H-Q3a75M1220, 15H-Q3a5M18, and H-Q5a75M18
models, which represent the cases with significant tidal disruption, weak
tidal disruption, and no tidal disruption, respectively. It shows that in
tidal disruptions, indeed most of the NS matter is swallowed by the BH soon
after the onset of the merger, i.e., within $\sim$1 ms, and the rest mass of
the matter remaining outside the BH depends strongly on the binary parameter.
For the case with significant tidal disruption (H-Q3a75M1220), 1%–20% of the
mass surrounds the BH to form a remnant disk, but only matter less than
$0.01\%$ remains outside the BH after the merger for the no tidal disruption
case (H-Q5a75M18).

At the time of the merger, a fraction of the matter outside the BH becomes
gravitationally unbound, i.e., the dynamical ejecta. The region in Fig. 1
surrounded by the dashed curves denotes the matter with gravitationally
unbound orbits. Figure 1 shows that the ejecta is formed around the outer edge
of the tidal tail. The time evolution of the ejecta mass is shown in Fig. 3.
The ejecta are formed within the dynamical time scale ($\sim 1$ ms) after the
onset of the merger, and, as is the case for the remnant matter located
outside the BH, the amount of the ejecta depends strongly on the binary
parameters.

Figure 2: Time evolution of $\hat{M}_{\rm{rem}}$ for representative models,
H-Q3a75M1220, 15H-Q3a5M18, and H-Q5a75M18 (solid, dot-dashed, and dashed
lines, respectively). Each model represents the cases with significant tidal
disruption, weak tidal disruption, and no tidal disruption, respectively.
Figure 3: Time evolution of $\hat{M}_{\rm{eje}}$ for the same models as in
Fig. 2. The $y$-axis shows the normalized ejecta mass. The tiny amount of
ejecta before and right after the merger is due to the mass conservation error
around the NS atomsphere.

###  III.2 Fixed compactness case

####  III.2.1 Properties of remnant disks and ejecta

Table 4: Characteristic physical quantities for the model with $C_{\rm{NS}}=0.182$. $\hat{M}_{\rm{rem}}$ is the normalized remnant mass for the matter located outside the apparent horizon, $\hat{M}_{\rm{eje}}$ is the ejecta mass during the merger, normalized by the baryon mass of the NS $M_{\rm{b}}$. $f_{\rm{cut}}m_{0}$ is the normalized cutoff frequency. $T_{\rm{eje}}$ and $v_{\rm{eje}}$ are the kinetic energy and velocity of the ejecta. $\chi_{\rm{BH}}^{f}$ and $M_{\rm{BH}}^{f}$ are the dimensionless spin and mass of BHs after the merger. Models with EOS ”B” are taken from Ref. [45], with $C_{\rm{NS}}=0.1819$, $M_{\rm{NS}}=1.35M_{\odot}$, and $N=50$. The significant figures of some quantities from Ref. [45] are different from that in our paper. All physical quantities are evaluated at 10 ms after the onset of merger. The grid resolution is $N=82$. Model ID |  | EOS | $\hat{M}_{\rm{rem}}$ | $\hat{M}_{\rm{eje}}$ | $f_{\rm{cut}}m_{0}$ | $T_{\rm{eje}}[\rm{erg}]$ | $v_{\rm{eje}}[c]$ | $\chi_{\rm{BH}}^{f}$ | $M_{\rm{BH}}^{f}[M_{\odot}]$
---|---|---|---|---|---|---|---|---|---
1 | Q3a5 | 15H | 0.052 | 0.0015 | 0.079 | 5.8$\times 10^{49}$ | 0.15 | 0.768 | 6.48
2 | HB | 0.047 | 0.0015 | 0.079 | 5.3$\times 10^{49}$ | 0.16 | 0.769 | 5.48
| Ref. [45] | B | 0.033 | N/A | N/A | N/A | N/A | 0.77 | N/A
3 | Q3a75 | 15H | 0.12 | 0.0046 | 0.071 | 2.3$\times 10^{50}$ | 0.18 | 0.871 | 6.37
4 | HB | 0.12 | 0.0059 | 0.069 | 2.8$\times 10^{50}$ | 0.18 | 0.870 | 5.38
5 | APR4 | 0.12 | 0.0065 | 0.070 | 3.1$\times 10^{50}$ | 0.19 | 0.870 | 5.15
6 | H4 | 0.13 | 0.0042 | 0.069 | 2.3$\times 10^{50}$ | 0.18 | 0.871 | 6.22
| Ref. [45] | B | 0.10 | N/A | N/A | N/A | N/A | 0.86 | N/A
7 | Q5a75 | 15H | 0.037 | 0.0047 | 0.101 | 2.8$\times 10^{50}$ | 0.19 | 0.845 | 9.81
8 | HB | 0.036 | 0.0054 | 0.101 | 2.8$\times 10^{50}$ | 0.19 | 0.845 | 8.29
| Ref. [45] | B | 0.021 | N/A | N/A | N/A | N/A | 0.85 | N/A
Table 5: Characteristic physical quantities for the models with $C_{\rm{NS}}=0.147$. Models with EOS ”H4” are taken from Ref. [32], with $C_{\rm{NS}}=0.147$, $M_{\rm{NS}}=1.35M_{\odot}$, and $N=60$. The significant figures of some quantities from Ref. [32] are different from that in our paper. $T_{\rm{eje}}$ and $v_{\rm{eje}}$ from Ref. [32] are not listed in the table since we use different definitions, which will be discussed in Appendix. B. Other quantities have the same definitions as in Table 4. Model ID |  | EOS | $\hat{M}_{\rm{rem}}$ | $\hat{M}_{\rm{eje}}$ | $f_{\rm{cut}}m_{0}$ | $T_{\rm{eje}}[\rm{erg}]$ | $v_{\rm{eje}}[c]$ | $\chi_{\rm{BH}}^{f}$ | $M_{\rm{BH}}^{f}[M_{\odot}]$
---|---|---|---|---|---|---|---|---|---
9 | Q3a5 | H | 0.16 | 0.024 | 0.040 | 9.4$\times 10^{50}$ | 0.18 | 0.755 | 4.59
10 | B | 0.16 | 0.022 | 0.041 | 8.4$\times 10^{50}$ | 0.19 | 0.758 | 4.11
| Ref. [32] | H4 | 0.16 | 0.02 | N/A | N/A | N/A | 0.76 | 5.07
11 | Q3a75 | 15H | 0.22 | 0.030 | 0.038 | 1.3$\times 10^{51}$ | 0.18 | 0.869 | 5.07
12 | H | 0.21 | 0.030 | 0.039 | 1.2$\times 10^{51}$ | 0.19 | 0.871 | 4.54
13 | B | 0.21 | 0.031 | 0.039 | 1.2$\times 10^{51}$ | 0.19 | 0.871 | 4.07
14 | APR4 | 0.21 | 0.031 | 0.039 | 1.2$\times 10^{51}$ | 0.19 | 0.870 | 4.07
| Ref. [32] | H4 | 0.22 | 0.03 | N/A | N/A | N/A | 0.88 | 4.99
15 | Q5a75 | H | 0.21 | 0.040 | 0.059 | 2.2$\times 10^{51}$ | 0.21 | 0.830 | 6.93
16 | B | 0.21 | 0.039 | 0.061 | 2.0$\times 10^{51}$ | 0.22 | 0.831 | 6.21
| Ref. [32] | H4 | 0.22 | 0.03 | N/A | N/A | N/A | 0.83 | 7.65

We list the key results of the simulations for C182 and C147 models in Tables
4 and 5, respectively. We also include the results for $\hat{M}_{\rm{rem}}$
and $\hat{M}_{\rm{eje}}$ with the same compactness from Refs. [32, 45]. Tables
4 and 5 show that $\hat{M}_{\rm{rem}}$ increases as $\chi_{\rm BH}$ increases,
and as the compactness $C_{\rm{NS}}$ and the mass ratio $Q$ decreases. This is
consistent with the predictions given in Sec. III.1. However, we note that the
dependence of $\hat{M}_{\rm{rem}}$ on $Q$ is rather weak in the case of the
lower compactness $C_{\rm{NS}}\approx 0.147$, which deviates from the naive
prediction of Newton approximation shown in Eq. (21). This weaker dependence
is also found and an interpretation for this is described in Ref. [34].

From Table 5 and 4, it is found that $\hat{M}_{\rm{eje}}$ increases as
$\chi_{\rm BH}$ increases and as $C_{\rm{NS}}$ decreases. This is the same
dependence as $\hat{M}_{\rm{rem}}$. However, being different from
$\hat{M}_{\rm{rem}}$, $\hat{M}_{\rm{eje}}$ increases as $Q$ increases. This
has also been found in the previous works [29, 42, 32, 35, 38, 41, 33],
indicating that the dynamics of the tidally disrupted matter in the vicinity
of the BH plays a key role in determining the ejecta mass.

Figure 4: The values of $\hat{M}_{\rm{rem}}$ for the models with
$C_{\rm{NS}}=0.182$ and $0.147$. The models with the same values of
$\chi_{\rm{BH}},Q$, and $C_{\rm{NS}}$ but with different NS mass $M_{\rm{NS}}$
are grouped with the same colors. The error bar of each model is determined by
the difference between $N=62$ and $N=82$ results. Model ID corresponds to the
column model ID in Tables 4 and 5.

Fig. 4 compares the values of $\hat{M}_{\rm{rem}}$ among the models with the
same values of $Q$, $\chi_{\rm{BH}}$, and $C_{\rm{NS}}$ but with different NS
mass $M_{\rm{NS}}$, i.e., with different NS EOSs. The error bars are taken to
be the discrepancy of the results between $N=62$ and $N=82$ runs of the same
model. Hydrodynamic quantities do not always show monotonic dependency on the
gird resolution and it is difficult to evaluate the error by checking the
convergence. Therefore we take the error estimation method shown above. We
note that the discrepancies between $N=82$ and $N=102$ results are mostly
smaller than those between $N=62$ and $N=82$ results (see Appendix. A). If the
range of error bars for different models with the same values of
$\chi_{\rm{BH}},Q$, and $C_{\rm{NS}}$ overlaps, we consider that the results
are in agreement within the range of the accuracy.

It is shown that, within the numerical accuracy, $\hat{M}_{\rm{rem}}$ with the
same values of $Q$, $\chi_{\rm{BH}}$, and $C_{\rm{NS}}$ agrees with each other
regardless of their different NS masses. We also find from Tables 4 and 5 that
the dimensionless spin of the remnant BH, $\chi_{\rm{BH}}^{f}$, is
approximately the same among the models with the same initial values of $Q$,
$\chi_{\rm{BH}}$, and $C_{\rm{NS}}$. This indicates that the BH swallows
approximately the same amount of the mass and angular momentum normalized by
$M_{\rm{NS}}$. These results show that $\hat{M}_{\rm{rem}}$ and the radius at
which the NS tidal disruption occurs depends sensitively on $Q$, $\chi_{\rm
BH}$, and $C_{\rm{NS}}$ while $M_{\rm{NS}}$ plays a role only in scaling the
system.

Figure 5: The values of $\hat{M}_{\rm{eje}}$ for models with
$C_{\rm{NS}}=0.182$ and $0.147$. The models with the same values of
$\chi_{\rm{BH}},Q$, and $C_{\rm{NS}}$ but with different NS mass $M_{\rm{NS}}$
are grouped with the same colors. The error bar of each model is determined by
the difference between $N=62$ and $N=82$ results. Model ID corresponds to the
column model ID in Tables 4 and 5.

Fig. 5 shows the results for the ejecta mass, with its error bar defined in
the same way as Fig. 4. $\hat{M}_{\rm{eje}}$ with the same values of $Q$,
$\chi_{\rm{BH}}$, and $C_{\rm{NS}}$ also agrees with each other. The results
show that $\hat{M}_{\rm{eje}}$ also depends sensitively on $Q$, $\chi_{\rm
BH}$, and $C_{\rm{NS}}$ but only minorly on $M_{\rm{NS}}$.

It is worth noting that the values of $\hat{M}_{\rm{rem}}$ and
$\hat{M}_{\rm{eje}}$ for the models which employ the four-pieces EOSs (model
ID $=5,6$, and 14) are also approximately in agreement with the results for
the models which employ the two-pieces EOSs with the same values of $Q$,
$\chi_{\rm BH}$, and $C_{\rm{NS}}$. This implies that the detailed structure,
i.e., the value of $\Gamma_{2}$ in Table 2, has only a minor effect on the
results of $\hat{M}_{\rm{rem}}$ and $\hat{M}_{\rm{eje}}$. This conclusion also
justifies us to study NSs with $M_{\mathrm{NS}}=1.8M_{\odot}$ using two-pieces
polytropes in Sec. III.3.

Figure 6: The value of $v_{\rm{eje}}$ at 10 ms after the merger for models
with $C_{\rm{NS}}=0.182$ and $0.147$. The models with the same values of
$\chi_{\rm{BH}},Q$, and $C_{\rm{NS}}$ but with different NS mass $M_{\rm{NS}}$
are grouped with the same colors. The error bar of each model is determined by
the difference between $N=62$ and $N=82$ results. Model ID corresponds to the
column model ID in Tables 4 and 5.

Fig. 6 shows the average ejecta velocity $v_{\rm eje}$ and their error bars.
Broadly speaking, $v_{\rm{eje}}$ increases with $Q$ and depends weakly on
other parameters. This is consistent with the results of the previous studies
for $M_{\mathrm{NS}}\approx 1.4\,M_{\odot}$ [33, 32, 35, 89], which is also
described in fitting formula of Eq. (20). As is the case for
$\hat{M}_{\rm{rem}}$ and $\hat{M}_{\rm{eje}}$, given the same values of $Q$,
$\chi_{\rm{BH}}$, and $C_{\rm{NS}}$, the values of $v_{\rm{eje}}$ for
different models approximately agree with each other within the numerical
accuracy.

The coefficients of the fitting formula of Eq. (20) given in Ref. [33] were
obtained by fitting the average ejecta velocity $v_{\rm eje,old}$ obtained by
various configurations of BH-NS binaries. They assumed a linear relation
between $v_{\rm eje,old}$ and $Q$: $v_{\rm eje,old}=(a_{1}Q+a_{2})c$, and used
the least squares fitting to derive the coefficients $a_{1}$ and $a_{2}$.
However, the definition of $v_{\rm eje,old}$ in Ref. [33], which is different
from the present definition in Eq. (5), did not take into account the
gravitational potential. Thus, the value of $v_{\rm eje,old}$ can be
overestimated.

Here, we correct the effect of the gravitational potential energy in the data
summarized in Ref. [33] by employing the method introduced in Ref. [34], and
recalibrate the ejecta-velocity fitting formula employing those corrected data
as well as the results obtained in this paper. Our method to correct the
ejecta velocity $v_{\rm eje,cor}$ is described in Appendix B. We note that,
since the estimated value of the ejecta velocity for the models with small
ejecta mass is not reliable, they employed only the models with $\hat{M}_{\rm
eje}>0.003$. In this paper, we restrict the models in the same manner.

By employing all the data, including $v_{\rm{eje,cor}}$ in Ref. [33] and
$v_{\rm{eje}}$ in this paper, we obtain the fitting formula as

$v_{\rm{eje,fit}}=(0.01108Q+0.1495)c.$ (22)

Fig. 7 shows that the average ejecta velocity in Ref. [33] is generally larger
than $v_{\rm eje}$ obtained in this paper, while after correction,
$v_{\rm{eje,cor}}$ is generally closer to $v_{\rm eje}$. This indeed shows
that the values of $v_{\rm eje,old}$ in Ref. [33] are overestimated by
neglecting the effect of the gravitational potential energy. The largest
relative fitting error of the fitting formula is around $10\%$ for Eq. (20)
and is around $18\%$ for Eq. (22). On the other hand, the relative numerical
error of $v_{\rm{eje}}$ in this paper estimated from the different resolution
runs is around $10\%$. This implies that there is a approximate linear
relation between $v_{\rm{eje}}$ and $Q$. It may indicate that the ejecta
velocity is not only determined by the mass ratio $Q$, but also influenced by
other parameters, such as the BH spin, NS compactness, and initial spin
misalignment angle.

Figure 7: Original average ejecta velocity in Ref. [33], corrected velocity
$v_{\rm{eje,cor}}$, and $v_{\rm{eje}}$ obtained in this paper as functions of
mass ratio $Q$. “KK” denotes the data obtained from Ref. [33], “KK-cor”
denotes the same data as “KK” but with the gravitational potential energy
correction, $v_{\rm eje,cor}$, “paper” denotes the data newly obtained in this
paper. The line with a label “KK” denotes the fitting formula obtained in the
previous study [33] (Eq. (20)). The line with the label “all” denotes the
fitting formula obtained by employing both data of “KK-cor” and “paper” (Eq.
(22)). The data points are slightly shifted in $x$-direction to make different
sets of points separated.

Fig. 8 shows the velocity distributions of the ejecta normalized by its mass
for selected models. The distributions are evaluated at $10$ ms after the
onset of the merger. The models shown in the figure are 15H-Q3a75M1691 and
HB-Q3a75M1428 which have $(Q,\chi_{\mathrm{BH}},C_{\rm NS})=(3,0.75,0.182)$ in
common, and H-Q5a75M1220 and B-Q5a75M1092 which have
$(Q,\chi_{\mathrm{BH}},C_{\rm NS})=(5,0.75,0.147)$ in common. The
distributions we obtained from our simulations show similar features to those
obtained by Refs. [32, 89, 36]; they peak at $v_{\rm eje}\approx 0.1c$–$0.3c$
and have a steep cut-off for both the lower and higher velocity sides. The
models with the same values of $(Q,\chi_{\mathrm{BH}},C_{\rm NS})$ essentially
show the same distributions having a peak at $v_{\rm eje}\approx 0.15c$ and
$0.2c$ for the models with $(Q,\chi_{\mathrm{BH}},C_{\rm NS})=(3,0.75,0.182)$
and $(5,0.75,0.147)$, respectively. The NS mass $M_{\rm NS}$ does not play a
role in changing the result significantly.

Figure 8: Ejecta velocity distributions normalized by the ejecta mass for
selected models. The models 15H-Q3a75M1691 and HB-Q3a75M1428 have
$(Q,\chi_{\mathrm{BH}},C_{\rm NS})=(3,0.75,0.182)$ in common, and the models
H-Q5a75M1220 and B-Q5a75M1092 have $(Q,\chi_{\mathrm{BH}},C_{\rm
NS})=(5,0.75,0.147)$ in common. The models with the same values of
$(Q,\chi_{\mathrm{BH}},C_{\rm NS})$ approximately show the same distributions.

We also perform a careful resolution study in this paper for the ejecta
velocity distribution and find a result worthy to be remark. Fig. 9 shows the
ejecta velocity distribution for model H-Q3a75M1220 with resolutions N62, N82,
and N102. It shows that for relatively low values of the escape velocity
convergence is achieved. However, for the high-velocity side with $\gtrsim
0.3c$ the convergence is not fully achieved. Specifically, for higher grid
resolutions, more matter is present for the higher velocity side. Our
interpretation for this is that with the higher grid resolution, the surface
of the one-armed spiral, where the high-velocity ejecta exists, is better
resolved and the pressure gradient plays a role in additionally accelerating
the matter there.

Figure 9: Ejecta velocity distributions normalized by the ejecta mass for
model 15H-Q3a75M18 with resolutions N62, N82, and N102. For the higher
resolutions, the distribution is extended to the higher velocity side.

####  III.2.2 Properties of gravitational waves

Figure 10: The example that shows how we determine the value of $f_{\rm
cut}m_{0}$. The blue curve is the GW spectrum obtained from model
HB-Q5a75M1428, and the purple line is half of the GW spectrum obtained from
the quadrupole formula. We identify the intersection frequency as the value of
the cutoff frequency $f_{\rm cut}m_{0}$. Figure 11: The normalized cutoff
frequency $f_{\rm{cut}}m_{0}$ for models with the same values of $Q$,
$\chi_{\rm{BH}}$, $C_{\rm{NS}}$, and different EOSs. The error bar of each
model is taken to be the discrepancy of the results between $N=62$ and $N=82$
runs. Model ID corresponds to the column model ID in Tables 4 and 5.

The location of the cutoff frequency in the amplitude of the GW spectrum at
$f\gtrsim 1$ kHz contains rich information on the last stage of the inspiral
phase and mainly depends on the orbital frequency of the tidal disruption[93,
21, 94]. If the NS is not tidally disrupted, the GW spectrum is characterized
by the inspiral, merger, and ringdown waveforms, in which case the
gravitational waveforms are nearly identical to those for binary BH mergers
[95]. When tidal disruption happens, the wave amplitude quickly decreases
before the ringdown phase, and the wave amplitude by the inspiral motion shuts
down. This is because in the tidal disruptions, the NS matter becomes less
dense and diffused, and the excitation of the quasinormal modes is suppressed
by the phase cancellation [96, 97, 98, 29]. The GW spectra in tidal
disruptions are characterized by their damping above the cutoff frequency,
$f_{\rm{cut}}$. To quantify $f_{\rm{cut}}$, previous studies fit the GW
spectra by a function with seven parameters and define $f_{\rm{cut}}$ as one
of the parameters [29, 31, 45]. However, these two definitions are not robust
because the fitting formula is not flexible enough.

For $f\rightarrow 0$, $f|h(f)|$ is proportional to $f^{-1/6}$ and can be
approximated by the quadrupole formula $fh_{\rm quad}(f)$ (e.g., Ref. [99]).
Therefore, we determine the value of $f_{\rm{cut}}$ to be the value of the
intersection of the half of the GW spectrum obtained by the quadrupole formula
and the one obtained from simulations; see Fig. 10 as an example.

Fig. 11 shows the normalized cutoff frequency $f_{\rm cut}m_{0}$ and its error
bars for models with the same values of $Q$, $\chi_{\rm{BH}}$, and
$C_{\rm{NS}}$, and different EOSs. As in the cases of $M_{\rm rem}$ and
$M_{\rm eje}$ (see Sec. III.2.1), the error bars of $f_{\rm cut}m_{0}$ are
also taken to be the discrepancies of the results between $N=62$ and $N=82$
runs. By comparing the models with different values of $Q$ or $\chi_{\rm{BH}}$
or $C_{\rm{NS}}$, the figure shows that larger compactness $C_{\rm{NS}}$ and
higher mass ratio $Q$ induce higher cutoff frequency, while higher BH spin
$\chi_{\rm BH}$ induces lower cutoff frequency. These results are consistent
with Eq. (21), since larger values of $C_{\rm{NS}}$ make the tidal disruption
of the binary occur at a more inner orbit with a higher orbital frequency,
higher values of $Q$ induce weaker tidal force in the vicinity of $R_{\rm
ISCO}$, and higher values of $\chi_{\rm BH}$ gives larger spin-orbit
repulsion, thus reduces the orbital angular velocity to maintain a circular
orbit, decreases the radius of the ISCO, and therefore, makes NS more easily
to be disrupted. These parameter dependence of the cutoff frequency shown in
Fig. 11 are all consistent with the interpretation based on the comparison
between $R_{\rm ISCO}$ and $R_{\rm dis}$ in Eq. (21).

It is also seen that within the numerical accuracy, the normalized cutoff
frequency $f_{\rm cut}m_{0}$ for models with the same values of $Q$,
$\chi_{\rm{BH}}$, and $C_{\rm{NS}}$ is approximately identical even though
these models have different NS masses, i.e. different EOSs. The relative
deviation between $f_{\rm cut}m_{0}$ within the same group is small, and as a
quantity characterizing where tidal disruptions happen, the value of
$f_{\rm{cut}}m_{0}$ depends solely on the ratio of $R_{\rm{dis}}$ and
$R_{\rm{ISCO}}$. This dependence of $f_{\rm cut}m_{0}$ implies the relative
location of tidal disruptions depends solely on $Q$, $\chi_{\rm{BH}}$, and
$C_{\rm{NS}}$.

###  III.3 Massive NS case

Table 6 shows the simulation results with a fixed NS mass
$M_{\rm{NS}}=1.8M_{\odot}$ but with different EOSs. Massive NSs generally have
weak or no tidal disruptions due to their high compactness, as indicated in
our analysis in Sec. III.1.

Table 6: Characteristic physical quantities for the model with $M_{\rm{NS}}=1.8M_{\odot}$. Other quantities have the same definitions as in Table 4.  | EOS | $\hat{M}_{\rm rem}$ | $\hat{M}_{\rm fit}^{2012}$ | $\hat{M}_{\rm fit}^{2018}$ | $\hat{M}_{\rm eje}$ | $\hat{M}_{\rm fit}^{2016}$ | $\hat{M}_{\rm fit}^{2020}$ | $T_{\rm eje}[\rm erg]$ | $v_{\rm eje}[c]$ | $\chi_{\rm BH}^{f}$ | $M_{\rm BH}^{f}[M_{\odot}]$
---|---|---|---|---|---|---|---|---|---|---|---
Q3a0.5 | 15H | 0.014 | 1.88$\times 10^{-3}$ | 0.0117 | 3.3$\times 10^{-4}$ | 0 | 0 | 1.0$\times 10^{49}$ | 0.13 | 0.770 | 6.95
H | 5.8$\times 10^{-5}$ | 0 | 0 | 1.1$\times 10^{-6}$ | 0 | 0 | 3.8$\times 10^{46}$ | 0.13 | 0.761 | 6.93
Q3a0.75 | 15H | 0.081 | 0.0945 | 0.0779 | 1.5$\times 10^{-3}$ | 4.56$\times 10^{-3}$ | 1.49$\times 10^{-3}$ | 6.9$\times 10^{49}$ | 0.16 | 0.870 | 6.83
H | 8.9$\times 10^{-3}$ | 0.0283 | 0.0310 | 1.1$\times 10^{-4}$ | 5.48$\times 10^{-3}$ | 0 | 3.7$\times 10^{48}$ | 0.14 | 0.866 | 6.90
B | 1.1$\times 10^{-5}$ | 0 | 7.73-04 | 5.3$\times 10^{-7}$ | 0.0187 | 0 | 2.3$\times 10^{46}$ | 0.15 | 0.858 | 6.87
Q5a0.75 | 15H | 6.3$\times 10^{-3}$ | 0 | 0.0244 | 9.2$\times 10^{-4}$ | 0 | 0 | 5.6$\times 10^{49}$ | 0.18 | 0.846 | 10.48
H | 9.5$\times 10^{-6}$ | 0 | 6.58$\times 10^{-6}$ | 2.6$\times 10^{-6}$ | 0 | 0 | 1.4$\times 10^{47}$ | 0.17 | 0.844 | 10.47

To test the validity of fitting formulas derived in previous studies in more
compact regions quantitatively, we calculate the $\chi$-square defined as

$\chi^{2}=\frac{1}{N_{\text{model}}-N_{\text{para}}}\sum_{i=1}^{N_{\text{model}}}\Delta_{i,\rm{fit}}^{2},$
(23)

where $N_{\text{model}}$ is the number of models in simulations,
$\Delta_{i,\rm{fit}}$ is defined in Eq. (24), and $N_{\text{para}}$ is the
number of parameters used in the fitting formula.

$\Delta_{\rm{fit}}$ denotes the ratio of discrepancy of the results between
the fitting formulas and the simulation data to the estimated numerical error
of the simulation $\Delta_{\rm{NR}}$, which is defined by

$\Delta_{\rm{fit}}=\frac{\hat{M}_{\rm{NR}}-\hat{M}_{\rm{fit}}}{\Delta_{\rm{NR}}},$
(24)

where $\hat{M}_{\rm{NR}}$ and $\hat{M}_{\rm{fit}}$ denote the results from
numerical simulations and the prediction of the fitting formulas,
respectively. In Refs. [33, 46], $\Delta_{\rm{NR}}$ is assumed to be

$\Delta_{\rm{NR}}=\sqrt{(\hat{M}_{\rm{NR}}/10)^{2}+(1/50)^{2}},$ (25)

which holds for eje_2016 and eje_2020, while Refs. [38, 47] give

$\Delta_{\rm{NR}}=\sqrt{(\hat{M}_{\rm{NR}}/10)^{2}+(1/100)^{2}},$ (26)

which is the case for rem_2012 and rem_2018.

Table 7: The $\chi$-squares for different fitting formulas. M18, C182, and C147 denote models with $M_{\rm{NS}}=1.8M_{\odot},C_{\rm{NS}}=0.182$, and $C_{\rm{NS}}=0.147$, respectively. The first half is $\chi$-square calculated using $\Delta_{\rm{NR}}$ from Refs. [33, 38] to compare with the results in the previous simulations. The second half is $\chi$-square calculated using $\Delta_{\rm{NR}}$ obtained in our simulations. | $\chi_{\rm{rem}}^{2012}$ | $\chi_{\rm{rem}}^{2018}$ | $\chi_{\rm{eje}}^{2016}$ | $\chi_{\rm{eje}}^{2020}$
---|---|---|---|---
$\Delta_{\rm{NR}}$ given in Eqs. (25), (26)
Models in Ref. [32] | 4.75 | 5.56 | 0.0250 | 0.0601
M18 | 1.35 | 2.75 | 0.968 | 0.00122
C182 | 0.601 | 3.38 | 0.0320 | 0.0300
C147 | 1.77 | 5.47 | 0.0499 | 0.0162
$\Delta_{\rm{NR}}$ given in Eq. (27)
M18 | 11.7 | 28.6 | 43.0 | 0.0542
C182 | 2.53 | 10.9 | 1.39 | 1.31
C147 | 2.15 | 6.95 | 1.20 | 0.344

Since the fitting formulas were derived by minimizing the $\chi$-square in
previous simulations, to achieve a benchmark of their accuracy, we first
calculate the $\chi$-square of the fitting formulas with $\Delta_{\rm{NR}}$
from Eq. (25) for the ejecta mass and from Eq. (26) for the remnant mass for
numerical-relativity results given in Ref. [32], of which parameters of BH-NS
binaries are within the calibrated range. The results are shown in Table 7.
These results can be considered as typical values of the $\chi$-square. To
evaluate how well the fitting formulas work for our simulations, we compare
the $\chi$-square of our simulation data with that of Ref. [32].

Next, we calculate the $\chi$-square for our simulation data with
$M_{\rm{NS}}=1.8M_{\odot}$, $C_{\rm{NS}}$ ranging from 0.194 to 0.252, and
C182, C147 models with Eq. (23) and $\Delta_{\rm{NR}}$ given in Eqs. (25) and
(26). For the remnant mass, the $\chi$-squares of M18, C147, and C182 for both
fitting formulas are smaller than those of the models in Ref. [32]. For
rem_2012, its $\chi$-squares for our simulations are generally smaller than
previous simulations and close to the unity. $\chi_{\rm{rem}}^{2018}$
increases as the NS compactness decreases, and is generally larger than the
corresponding $\chi_{\rm{rem}}^{2012}$. As to the ejecta mass, for C182 and
C147 models, their $\chi$-squares are smaller than those of models in Ref.
[32]. For the M18 case, its $\chi_{\rm{eje}}^{2016}$ is significantly larger
than those of the other three cases, while its $\chi_{\rm{eje}}^{2020}$ is
much smaller. Nevertheless, we find that within the numerical accuracy
analyzed in previous studies, their fitting formulas for $\hat{M}_{\rm{rem}}$
and $\hat{M}_{\rm{eje}}$ work well even in the compact NS region since the
$\chi$-squares are close to unity.

We then calculate the $\chi$-square of fitting formulas using
$\Delta_{\rm{NR}}$ given by our simulations to test whether there could be
room for improving the fitting formulas employing higher resolution data of a
numerical-relativity simulation. In our research, the numerical error is
estimated to be

$\Delta_{\rm{NR}}=\sqrt{(\hat{M}_{\rm{NR}}/10)^{2}+(3/1000)^{2}}$ (27)

(see Appendix A). Eq. (27) shows that the numerical errors of our simulations
are smaller than those given by Eqs. (25) and (26), and hence, the
$\chi$-square calculated based on it gives more severe examination to the
fitting formula.

Results are shown in Table 7. The $\chi$-squares for M18 are significantly
larger than those of C147 and C182 except for eje_2020, which are much
smaller. The reason for this increase is that mergers with compact NSs result
in smaller values of $\hat{M}_{\rm{rem}}$ and $\hat{M}_{\rm{eje}}$ and
therefore sensitive to the lower limit of the estimated numerical error.
Assuming a smaller numerical error enhances their $\chi$-squares
significantly. $\chi^{2018}_{\rm{rem}}$ of C182, C147 also increases
significantly, while other $\chi-$squares for C147 and C182 models remain
close to the unity.

Therefore, if we assume a smaller numerical error of our simulations,
rem_2012, rem_2018, and eje_2016 have room for improvement, particularly for
compact NS cases. Their $\chi$-squares increase as the NS compactness
increases, and for the M18 case, their $\chi$-squares are significantly
larger. This is reasonable because these fitting formulas are calibrated in
the parameter space of $C_{\rm{NS}}<0.18$. Only $\chi_{\rm{eje}}^{2020}$
decreases as the compactness increases and has a better performance in the
compact NS region. Thus, if we want to estimate the remnant disk mass for
$C_{\rm{NS}}>0.194$ at higher accuracy, existing fitting formulas require
improvement. On the other hand, the fitting formula of the ejecta mass given
in Ref. [46] still works well in this region even if we assume a smaller
numerical error of our simulations since its $\chi$-square is less than the
unity. In addition, if we require a higher accuracy of the estimation of the
remnant mass, rem_2012 works better than rem_2018.

##  IV Conclusion & Discussion

Figure 12: Normalized remnant mass as a function of the compactness. For
models with the same values of $Q$, $\chi_{\rm{BH}}$, and $C_{\rm{NS}}$, since
the discrepancies between them are within the margin of error, the simulation
results are taken randomly from one of them, and the errorbar is taken to be
the largest error range of them.

We perform numerical-relativity simulations for various setups of BH-NS
mergers with $Q=3,5$, $\chi_{\rm{BH}}=0.5,0.75$, and
$C_{\rm{NS}}=0.147\mbox{--}0.252$ varying the NS mass. Our results with the
fixed NS compactness $C_{\rm{NS}}=0.182$ and $0.147$ show that given the
identical values of $Q$, $\chi_{\rm{BH}}$, and $C_{\rm{NS}}$,
$\hat{M}_{\rm{rem}}$ and $f_{\rm{cut}}m_{0}$ agree with each other within the
numerical accuracy of our study. We also found that this is the case for
$\hat{M}_{\rm{eje}}$, except for only one data set which is slightly out of
the numerical accuracy. Therefore, the hypothesis that $\hat{M}_{\rm{rem}}$,
$\hat{M}_{\rm{eje}}$, $v_{\rm{eje}}$ and $f_{\rm{cut}}m_{0}$ depend only on
$Q$, $\chi_{\mathrm{BH}}$, and $C_{\rm{NS}}$ is approximately true for the NS
mass in a wide range of $1.1\mbox{--}1.7\,M_{\odot}$. This justifies the
approach of studying the dependence of NS tidal disruptions on the NS
compactness by fixing the NS mass but changing the EOS, which is employed in
the previous studies.

We then performed numerical-relativity simulations of BH-NS mergers for a
large value of the NS mass ($M_{\rm NS}=1.8\,M_{\odot}$), and examined the
accuracy of the previous fitting formulas of $\hat{M}_{\rm{rem}}$,
$\hat{M}_{\rm{eje}}$, and ejecta velocity $v_{\rm{eje}}$ for the system with a
large value of the NS compactness ($C_{\rm NS}\geq 0.19$). As to the four
fitting formulas of $\hat{M}_{\rm{rem}}$ and $\hat{M}_{\rm{eje}}$, we find
that they still work well for compact NSs within the error range of their
studies. However, if we require a higher-precision prediction, the fitting
formulas of the remnant mass could give inaccurate predictions for
$C_{\rm{NS}}>0.194$ cases, leaving room for further improvements, and rem_2012
is generally more accurate than rem_2018. As to the ejecta mass, the fitting
formula eje_2020 from Ref. [46] still works well for more massive NSs with
$C_{\rm{NS}}$ between 0.194 and 0.252.

For the case of ejecta velocity, we corrected the data of average ejecta
velocity in Ref. [33] to take into account the gravitational potential effect,
and recalibrated the ejecta-velocity fitting formula employing these corrected
data as well as the results obtained in our present simulations. We obtained
the simple linear function as the fitting formula between $v_{\rm eje}$ and
$Q$ as in Ref. [33]. We found that the predicted ejecta velocity is
systematically smaller than those of the previous studies. The relative
fitting error of the fitting formula is around $10\%$ for the data in this
paper and in Ref. [33], and grows to $18\%$ if we use all the data to obtain a
new fitting formula. It is larger than $10\%$ of the relative numerical error
of $v_{\rm eje}$ estimated from different resolutions runs. This indicates
that the ejecta velocity has a dependence on the other parameters which cannot
be captured by the simple linear relations between $v_{\rm eje}$ and $Q$.

Our results give some guidance on how to improve the fitting formulas for the
remnant mass. From Fig. 12, which shows $\hat{M}_{\rm{rem}}$ as a function of
$C_{\rm{NS}}$, we can see that rem_2012 works best for Q3a5 models, even for
compact NSs. The inaccuracy of the fitting formula is caused mainly by the
large discrepancies between the fitting and the simulation results for Q3a75
models. The fitting formula generally overestimates the remnant mass for all
Q3a75 models, but this overestimation is only significant for mergers with
small remnant masses, i.e. with compact NSs. For example, H-Q3a75 has only
0.9% of its baryon mass remaining outside the BH after the merger, and the
rem_2012 gives a prediction of 2.8%. Rem_2018 overestimates the remnant mass
for compact NSs and is generally less accurate than rem_2012. Therefore, a
linear function of $C_{\rm{NS}}$ is still a good choice for the fitting
formula of the remnant mass and the second term in Eq. (11) requires
modification to get more accurate results.

The improvement in the fitting formulas in such a parameter region may be not
very important at this moment, since the overestimation is only significant
for the cases having weak tidal disruptions, for which only faint EM
counterparts can associate with the GW signals and are less likely to be
detected by observations. However, the improvement in such a parameter region
can play an important role in interpreting the observational data once faint
EM counterparts are detected or deep EM upper-limits are obtained in the
future BH-NS GW events (cf, GW190425 and GW190814 [100, 101]).

###### Acknowledgements.

Numerical computation was performed on Yukawa21 at Yukawa Institute for
Theoretical Physics, Kyoto University and the Yamazaki, Sakura, Cobra, and
Raven clusters at Max Planck Computing and Data Facility. This work was
supported by Grant-in-Aid for Scientific Research (20H00158, 23H04900, and
23H01172) of JSPS/MEXT. Shichuan Chen is supported by China Scholarship
Council.

##  Appendix A ESTIMATION OF ERROR CAUSED BY GRID RESOLUTION

Figure 13: Discrepancy of $M_{\rm{rem}}/M_{\rm{b}}$ and
$M_{\rm{eje}}/M_{\rm{b}}$ between runs with different grid resolutions. The
$x$-axis is taken to the results from $N=82$ runs. The scatters are taken to
be the minimum discrepancy among N62, N82, and N102 runs. The blue-dashed
curves are an upper estimation of the numerical error of the normalized
remnant and ejecta mass, and the orange-dashed lines are an estimation of the
relative numerical error.

The mass of the remnant disk and the ejecta are only a fraction of the NS
mass, and thus, can have large errors in the numerical results. Here we
estimate the errors due to finite grid resolutions, as shown in Fig. 13. Both
masses depend sensitively on the grid resolution.

As shown in the lower plots, for $\hat{M}>0.01M_{\odot}$, the relative errors
are generally less than 10%, while for very small values of $\hat{M}$, the
relative error can be up to 100%. We expect that it is possible to obtain a
function for the numerical errors, analogous to that given in Eq. (25) and Eq.
(26). Therefore, the coefficient in front of $\hat{M}$ in Eq. (27) can be
taken as 0.1, and the constant part is taken as 0.003. As shown in Fig. 13,
the error for $\hat{M}_{\rm{eje}}$ and $\hat{M}_{\rm{rem}}$ can be described
using a uniform estimation as it is in Eq. (27).

Table 8: Several numerical results for models 15H-Q3a75M18,H-Q3a75M1220, HB-Q3a75M1428, and HB-Q5a75M1428 with different grid resolutions $N=62,82$, and $102$. All the quantities are defined in the body text. N | $\hat{M}_{\rm rem}$ | $\hat{M}_{\rm eje}$ | $v_{\rm{eje}}$ | $f_{\rm{cut}}m_{0}$ | $\chi_{\rm{BH}}$
---|---|---|---|---|---
| 15H-Q3a75M18 |  |  |  |
62 | 0.0840 | 0.00144 | 0.155 | 0.0762 | 0.870
82 | 0.0807 | 0.00151 | 0.159 | 0.0849 | 0.870
102 | 0.0782 | 0.00203 | 0.172 | 0.0883 | 0.870
| H-Q3a75M1220 |  |  |  |
62 | 0.213 | 0.0290 | 0.181 | 0.0384 | 0.863
82 | 0.212 | 0.0296 | 0.185 | 0.0394 | 0.871
102 | 0.210 | 0.0303 | 0.188 | 0.0398 | 0.868
| HB-Q3a75M1428 |  |  |  |
62 | 0.128 | 0.00513 | 0.172 | 0.0666 | 0.873
82 | 0.124 | 0.00589 | 0.182 | 0.0687 | 0.870
102 | 0.122 | 0.005884 | 0.185 | 0.0701 | 0.870
| HB-Q5a75M1428 |  |  |  |
62 | 0.0393 | 0.000414 | 0.200 | 0.0993 | 0.845
82 | 0.0360 | 0.00542 | 0.191 | 0.1014 | 0.845
102 | 0.0327 | 0.00527 | 0.192 | 0.1016 | 0.845
| 15H-Q5a75M1691 |  |  |  |
62 | 0.0428 | 0.000463 | 0.177 | 0.0993 | 0.845
82 | 0.0369 | 0.00474 | 0.187 | 0.1011 | 0.845
102 | 0.0351 | 0.00491 | 0.185 | 0.1016 | 0.845

We list numerical results for some models with different grid resolutions
$N=62,\ 82$ and $102$ in Table 8. The relative error between N82 and N102 runs
ranges from 0.9% to 3% for $\hat{M}_{\rm{rem}}$, and ranges from 1% to 34% for
$\hat{M}_{\rm{eje}}$. The relative error of the ejecta velocity between $N=82$
and $N=62$ runs ranges from 2% to 5%. The discrepancies of $f_{\rm{cut}}m_{0}$
between $N=62$ and $N=82$ range from 0.4% to 5%. The discrepancies of the spin
are around 0.1% for the three resolutions. In general, due to the complexity
of hydrodynamics, hydrodynamic quantities have larger uncertainties.
$\hat{M}_{\rm{rem}}$ decreases as the grid resolution increases, while
systematic convergence properties of other quantities are not seen, and the
discrepancies of numerical results among different resolutions do not decrease
necessarily as the grid resolution increases. The reason for this unsystematic
behavior is likely that (i) the motion of the matters is affected by the
background atmosphere in the simulations [32] and (ii) in the presence of
shocks (e.g., in the tidal tail), numerical accuracy becomes the first-order
convergence and hence the error can be enhanced in an unpredicted manner.

##  Appendix B CORRECTION OF EJECTA VELOCITY

Figure 14: Time evolution of $v_{\rm eje}$ for model 15H-Q3a75M1691 with
different definitions. Here, the results are from $N=82$ runs. $v_{\rm eje}$
is defined by Eq. (20), $v_{\rm eje,old}$ by Eq. (31), and $v_{\rm eje,cor}$
by Eq. (32).

Different from Eq. (4), the kinetic energy of the ejecta defined in Ref. [33]
is

$T_{\rm eje,old}:=E_{\rm eje}-U_{\rm eje}-M_{\rm eje}.$ (28)

The total energy of the ejecta $E_{\rm eje}$ is defined by

$E_{\rm eje}:=\int_{-u_{t}>1,r>r_{\rm AH}}\rho_{*}\hat{e}\,d^{3}x.$ (29)

The internal energy of the ejecta $U_{\rm eje}$ is defined by

$U_{\rm eje}:=\int_{-u_{t}>1,r>r_{\rm AH}}\rho_{*}\epsilon\,d^{3}x.$ (30)

The mass of ejecta $M_{\rm eje}$ is defined as in Eq. (3). The definitions of
$\rho_{*}$, $\hat{e}$ and $\epsilon$ are listed in Table. 1.

Assuming the Newtonian dynamics, the averaged velocity of the ejecta can be
evaluated as

$v_{\rm eje,old}:=\sqrt{\frac{2T_{\rm eje,old}}{M_{\rm eje}}}.$ (31)

The gravitational potential energy is not excluded in $T_{\rm eje,old}$, and
$v_{\rm eje,old}$ is overestimated. To subtract the effect of this
gravitational potential energy, we can approximately correct the velocity as
in Ref. [34]

$v_{\rm eje,cor}:=\sqrt{v^{2}_{\rm eje,old}-\frac{2m_{0}}{v_{\rm
eje,old}(t-t_{\rm merger})}},$ (32)

where $v_{\rm eje,old}$ is evaluated at time $t$.

Fig. 14 shows the time evolution of $v_{\rm eje}$ with different definitions
for a representative model 15H-Q3a75M1691. If we choose $v_{\rm eje}$ defined
in Eq. (5) to be the fiducial value, we find that $v_{\rm eje,old}$ is still
overestimated by $\sim 50\%$ at 12 ms after the onset of the merger. The
corrected velocity of ejecta $v_{\rm eje,cor}$ is closer to $v_{\rm eje}$
compared to $v_{\rm eje,old}$, and asymptotically approaches $v_{\rm eje}$ at
10 ms after the onset of the merger. Thus, it is reasonable to use Eq. (32) to
correct previous data of the average velocity of ejecta.

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Observation of gravitational waves from a binary black hole merger, Phys. Rev. Lett. 116, 061102 (2016).
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Gw170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119, 161101 (2017a).
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ , Multi-messenger observations of a binary neutron star merger, Astrophys. J. 848, L12 (2017b).
  * Abbott  _et al._ [2017c] B. P. Abbott _et al._ , Gravitational waves and gamma-rays from a binary neutron star merger: GW170817 and GRB 170817a, Astrophys. J. 848, L13 (2017c).
  * Goldstein  _et al._ [2017] A. Goldstein, P. Veres, E. Burns, M. S. Briggs, R. Hamburg, D. Kocevski, C. A. Wilson-Hodge, R. D. Preece, S. Poolakkil, O. J. Roberts, C. M. Hui, V. Connaughton, J. Racusin, A. von Kienlin, T. D. Canton, N. Christensen, T. Littenberg, K. Siellez, L. Blackburn, J. Broida, E. Bissaldi, W. H. Cleveland, M. H. Gibby, M. M. Giles, R. M. Kippen, S. McBreen, J. McEnery, C. A. Meegan, W. S. Paciesas, and M. Stanbro, An ordinary short gamma-ray burst with extraordinary implications: Fermi-gbm detection of grb 170817a, The Astrophysical Journal Letters 848, L14 (2017).
  * LIGO Scientific Collaboration and Virgo Collaboration [2017] LIGO Scientific Collaboration and Virgo Collaboration, LIGO/Virgo G298048: Fermi GBM trigger 524666471/170817529: LIGO/Virgo Identification of a possible gravitational-wave counterpart, GRB Coordinates Network 21505, 1 (2017).
  * De  _et al._ [2018] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Tidal deformabilities and radii of neutron stars from the observation of gw170817, Phys. Rev. Lett. 121, 091102 (2018).
  * Abbott and et al [2018] B. P. Abbott and et al (The LIGO Scientific Collaboration and the Virgo Collaboration), Gw170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121, 161101 (2018).
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Properties of the binary neutron star merger gw170817, Phys. Rev. X 9, 011001 (2019).
  * Narikawa  _et al._ [2020] T. Narikawa, N. Uchikata, K. Kawaguchi, K. Kiuchi, K. Kyutoku, M. Shibata, and H. Tagoshi, Reanalysis of the binary neutron star mergers gw170817 and gw190425 using numerical-relativity calibrated waveform models, Phys. Rev. Res. 2, 043039 (2020).
  * Savchenko  _et al._ [2017] V. Savchenko, C. Ferrigno, E. Kuulkers, A. Bazzano, E. Bozzo, S. Brandt, J. Chenevez, T. J.-L. Courvoisier, R. Diehl, A. Domingo, L. Hanlon, E. Jourdain, A. von Kienlin, P. Laurent, F. Lebrun, A. Lutovinov, A. Martin-Carrillo, S. Mereghetti, L. Natalucci, J. Rodi, J.-P. Roques, R. Sunyaev, and P. Ubertini, Integral detection of the first prompt gamma-ray signal coincident with the gravitational-wave event gw170817, The Astrophysical Journal Letters 848, L15 (2017).
  * Tanaka  _et al._ [2013] M. Tanaka, K. Hotokezaka, K. Kyutoku, S. Wanajo, K. Kiuchi, Y. Sekiguchi, and M. Shibata, Radioactively powered emission from black hole-neutron star mergers, Astrophys. J. 780, 31 (2013).
  * Lattimer and Schramm [1974] J. M. Lattimer and D. N. Schramm, Black-Hole-Neutron-Star Collisions, Astrophys. J. Lett. 192, L145 (1974).
  * Eichler  _et al._ [1989] D. Eichler, M. Livio, T. Piran, and D. N. Schramm, Nucleosynthesis, neutrino bursts and $\gamma$-rays from coalescing neutron stars, Nature (London) 340, 126 (1989).
  * Freiburghaus  _et al._ [1999] C. Freiburghaus, S. Rosswog, and F.-K. Thielemann, r-process in neutron star mergers, The Astrophysical Journal 525, L121 (1999).
  * Drout  _et al._ [2017] M. R. Drout, A. L. Piro, B. J. Shappee, C. D. Kilpatrick, J. D. Simon, C. Contreras, D. A. Coulter, R. J. Foley, M. R. Siebert, N. Morrell, K. Boutsia, F. D. Mille, T. W.-S. Holoien, D. Kasen, J. A. Kollmeier, B. F. Madore, A. J. Monson, A. Murguia-Berthier, Y.-C. Pan, J. X. Prochaska, E. Ramirez-Ruiz, A. Rest, C. Adams, K. Alatalo, E. Bañados, J. Baughman, T. C. Beers, R. A. Bernstein, T. Bitsakis, A. Campillay, T. T. Hansen, C. R. Higgs, A. P. Ji, G. Maravelias, J. L. Marshall, C. M. Bidin, J. L. Prieto, K. C. Rasmussen, C. Rojas-Bravo, A. L. Strom, N. Ulloa, J. Vargas-González, Z. Wan, and D. D. Whitten, Light curves of the neutron star merger gw170817/sss17a: Implications for r-process nucleosynthesis, Science 358, 1570 (2017), https://www.science.org/doi/pdf/10.1126/science.aaq0049 .
  * Pian  _et al._ [2017] E. Pian, P. D’Avanzo, S. Benetti, M. Branchesi, E. Brocato, S. Campana, E. Cappellaro, V. D’Elia, J. Fynbo, F. Getman, G. Ghirlanda, G. Ghisellini, A. Grado, G. Greco, J. Hjorth, C. Kouveliotou, A. Levan, L. Limatola, and D. Vergani, Spectroscopic identification of r-process nucleosynthesis in a double neutron-star merger, Nature 551 (2017).
  * Abbott  _et al._ [2021] R. Abbott _et al._ , Observation of gravitational waves from two neutron star–black hole coalescences, Astrophys. J. 915, L5 (2021).
  * Abbott  _et al._ [2020a] B. P. Abbott _et al._ , GW190425: Observation of a compact binary coalescence with total mass $\sim$ 3.4 m $\odot$, The Astrophysical Journal Letters 892, L3 (2020a).
  * Kyutoku  _et al._ [2020a] K. Kyutoku, S. Fujibayashi, K. Hayashi, K. Kawaguchi, K. Kiuchi, M. Shibata, and M. Tanaka, On the possibility of GW190425 being a black hole–neutron star binary merger, Astrophys. J. 890, L4 (2020a).
  * Shibata and Taniguchi [2011] M. Shibata and K. Taniguchi, Coalescence of black hole-neutron star binaries, Living Rev. Relativity 14, 6 (2011).
  * Shibata and Uryū [2006] M. Shibata and K. Uryū, Merger of black hole-neutron star binaries: Nonspinning black hole case, Phys. Rev. D 74, 121503 (2006).
  * Shibata and Uryū [2007] M. Shibata and K. Uryū, Merger of black hole–neutron star binaries in full general relativity, Classical Quantum Gravity 24, S125 (2007).
  * Shibata and Taniguchi [2008] M. Shibata and K. Taniguchi, Merger of black hole and neutron star in general relativity: Tidal disruption, torus mass, and gravitational waves, Phys. Rev. D 77, 084015 (2008).
  * Etienne  _et al._ [2008] Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro, K. Taniguchi, and T. W. Baumgarte, Fully general relativistic simulations of black hole-neutron star mergers, Phys. Rev. D 77, 084002 (2008).
  * Duez  _et al._ [2008] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and S. A. Teukolsky, Evolving black hole-neutron star binaries in general relativity using pseudospectral and finite difference methods, Phys. Rev. D 78, 104015 (2008).
  * Yamamoto  _et al._ [2008] T. Yamamoto, M. Shibata, and K. Taniguchi, Simulating coalescing compact binaries by a new code (sacra), Phys. Rev. D 78, 064054 (2008).
  * Etienne  _et al._ [2009] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baumgarte, General relativistic simulations of black-hole–neutron-star mergers: Effects of black-hole spin, Phys. Rev. D 79, 044024 (2009).
  * Shibata  _et al._ [2009] M. Shibata, K. Kyutoku, T. Yamamoto, and K. Taniguchi, Gravitational waves from black hole-neutron star binaries: Classification of waveforms, Phys. Rev. D 79, 044030 (2009).
  * Foucart  _et al._ [2011] F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukolsky, Black hole-neutron star mergers: Effects of the orientation of the black hole spin, Phys. Rev. D 83, 024005 (2011).
  * Kyutoku  _et al._ [2010] K. Kyutoku, M. Shibata, and K. Taniguchi, Gravitational waves from nonspinning black hole-neutron star binaries: Dependence on equations of state, Phys. Rev. D 82, 044049 (2010).
  * Kyutoku  _et al._ [2015] K. Kyutoku, K. Ioka, H. Okawa, M. Shibata, and K. Taniguchi, Dynamical mass ejection from black hole-neutron star binaries, Phys. Rev. D 92, 044028 (2015).
  * Kawaguchi  _et al._ [2016] K. Kawaguchi, K. Kyutoku, M. Shibata, and M. Tanaka, Models of Kilonova/Macronova Emission from Black Hole-Neutron Star Mergers, Astrophys. J. 825, 52 (2016).
  * Hayashi  _et al._ [2021] K. Hayashi, K. Kawaguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Properties of the remnant disk and the dynamical ejecta produced in low-mass black hole-neutron star mergers, Phys. Rev. D 103, 043007 (2021).
  * Foucart  _et al._ [2019a] F. Foucart, M. D. Duez, L. E. Kidder, S. M. Nissanke, H. P. Pfeiffer, and M. A. Scheel, Numerical simulations of neutron star-black hole binaries in the near-equal-mass regime, Phys. Rev. D 99, 103025 (2019a).
  * Brege  _et al._ [2018] W. Brege, M. D. Duez, F. Foucart, M. B. Deaton, J. Caro, D. A. Hemberger, L. E. Kidder, E. O’Connor, H. P. Pfeiffer, and M. A. Scheel, Black hole-neutron star mergers using a survey of finite-temperature equations of state, Phys. Rev. D 98, 063009 (2018).
  * Foucart  _et al._ [2013a] F. Foucart, M. B. Deaton, M. D. Duez, L. E. Kidder, I. MacDonald, C. D. Ott, H. P. Pfeiffer, M. A. Scheel, B. Szilagyi, and S. A. Teukolsky, Black-hole–neutron-star mergers at realistic mass ratios: Equation of state and spin orientation effects, Phys. Rev. D 87, 084006 (2013a).
  * Foucart [2012] F. Foucart, Black-hole–neutron-star mergers: Disk mass predictions, Phys. Rev. D 86, 124007 (2012).
  * Foucart  _et al._ [2012] F. Foucart, M. D. Duez, L. E. Kidder, M. A. Scheel, B. Szilagyi, and S. A. Teukolsky, Black hole-neutron star mergers for $10\text{ }\text{ }{M}_{\bigodot}$ black holes, Phys. Rev. D 85, 044015 (2012).
  * Lovelace  _et al._ [2013] G. Lovelace, M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Massive disc formation in the tidal disruption of a neutron star by a nearly extremal black hole, Classical Quantum Gravity 30, 135004 (2013).
  * Foucart  _et al._ [2014] F. Foucart, M. B. Deaton, M. D. Duez, E. O’Connor, C. D. Ott, R. Haas, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi, Neutron star-black hole mergers with a nuclear equation of state and neutrino cooling: Dependence in the binary parameters, Phys. Rev. D 90, 024026 (2014).
  * Kawaguchi  _et al._ [2015] K. Kawaguchi, K. Kyutoku, H. Nakano, H. Okawa, M. Shibata, and K. Taniguchi, Black hole-neutron star binary merger: Dependence on black hole spin orientation and equation of state, Phys. Rev. D 92, 024014 (2015).
  * Kyutoku  _et al._ [2021] K. Kyutoku, M. Shibata, and K. Taniguchi, Coalescence of black hole–neutron star binaries, Living Rev. Rel. 24, 5 (2021).
  * Ferrari  _et al._ [2010] V. Ferrari, L. Gualtieri, and F. Pannarale, Neutron star tidal disruption in mixed binaries: The imprint of the equation of state, Phys. Rev. D 81, 064026 (2010).
  * Kyutoku  _et al._ [2011] K. Kyutoku, H. Okawa, M. Shibata, and K. Taniguchi, Gravitational waves from spinning black hole-neutron star binaries: dependence on black hole spins and on neutron star equations of state, Phys. Rev. D 84, 064018 (2011).
  * Krüger and Foucart [2020] C. J. Krüger and F. Foucart, Estimates for disk and ejecta masses produced in compact binary mergers, Phys. Rev. D 101, 103002 (2020).
  * Foucart  _et al._ [2018] F. Foucart, T. Hinderer, and S. Nissanke, Remnant baryon mass in neutron star-black hole mergers: Predictions for binary neutron star mimickers and rapidly spinning black holes, Phys. Rev. D 98, 081501 (2018).
  * Pannarale and Ohme [2014] F. Pannarale and F. Ohme, Prospects for joint gravitational-wave and electromagnetic observations of neutron-star–black-hole coalescing binaries, The Astrophysical Journal Letters 791, L7 (2014).
  * Radice  _et al._ [2018] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Gw170817: Joint constraint on the neutron star equation of state from multimessenger observations, The Astrophysical Journal Letters 852, L29 (2018).
  * Hinderer  _et al._ [2019a] T. Hinderer, S. Nissanke, F. Foucart, K. Hotokezaka, T. Vincent, M. Kasliwal, P. Schmidt, A. R. Williamson, D. A. Nichols, M. D. Duez, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Distinguishing the nature of comparable-mass neutron star binary systems with multimessenger observations: Gw170817 case study, Phys. Rev. D 100, 063021 (2019a).
  * Coughlin  _et al._ [2019] M. W. Coughlin, T. Dietrich, B. Margalit, and B. D. Metzger, Multimessenger Bayesian parameter inference of a binary neutron star merger, Monthly Notices of the Royal Astronomical Society: Letters 489, L91 (2019), https://academic.oup.com/mnrasl/article-pdf/489/1/L91/54701970/mnrasl_489_1_l91.pdf .
  * Barbieri  _et al._ [2019] C. Barbieri, O. S. Salafia, M. Colpi, G. Ghirlanda, A. Perego, and A. Colombo, Filling the mass gap: How kilonova observations can unveil the nature of the compact object merging with the neutron star, The Astrophysical Journal Letters 887, L35 (2019).
  * Andreoni  _et al._ [2020] I. Andreoni, D. A. Goldstein, M. M. Kasliwal, P. E. Nugent, R. Zhou, J. A. Newman, M. Bulla, F. Foucart, K. Hotokezaka, E. Nakar, S. Nissanke, G. Raaijmakers, J. S. Bloom, K. De, J. E. Jencson, C. Ward, T. Ahumada, S. Anand, D. A. H. Buckley, M. D. Caballero-García, A. J. Castro-Tirado, C. M. Copperwheat, M. W. Coughlin, S. B. Cenko, M. Gromadzki, Y. Hu, V. R. Karambelkar, D. A. Perley, Y. Sharma, A. F. Valeev, D. O. Cook, U. C. Fremling, H. Kumar, K. Taggart, A. Bagdasaryan, J. Cooke, A. Dahiwale, S. Dhawan, D. Dobie, P. Gatkine, V. Z. Golkhou, A. Goobar, A. G. Chaves, M. Hankins, D. L. Kaplan, A. K. H. Kong, E. C. Kool, S. Mohite, J. Sollerman, A. Tzanidakis, S. Webb, and K. Zhang, Growth on s190814bv: Deep synoptic limits on the optical/near-infrared counterpart to a neutron star–black hole merger, The Astrophysical Journal 890, 131 (2020).
  * Ascenzi  _et al._ [2019] S. Ascenzi, N. D. Lillo, C.-J. Haster, F. Ohme, and F. Pannarale, Constraining the neutron star radius with joint gravitational-wave and short gamma-ray burst observations of neutron star–black hole coalescing binaries, The Astrophysical Journal 877, 94 (2019).
  * Tauris  _et al._ [2017] T. M. Tauris, M. Kramer, P. C. C. Freire, N. Wex, H.-T. Janka, N. Langer, P. Podsiadlowski, E. Bozzo, S. Chaty, M. U. Kruckow, E. P. J. van den Heuvel, J. Antoniadis, R. P. Breton, and D. J. Champion, Formation of double neutron star systems, The Astrophysical Journal 846, 170 (2017).
  * Farrow  _et al._ [2019] N. Farrow, X.-J. Zhu, and E. Thrane, The mass distribution of galactic double neutron stars, The Astrophysical Journal 876, 18 (2019).
  * Demorest  _et al._ [2010] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, A two-solar-mass neutron star measured using shapiro delay, Nature 467, 1081–1083 (2010).
  * Antoniadis  _et al._ [2013] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, and D. G. Whelan, A massive pulsar in a compact relativistic binary, Science 340, 1233232 (2013), https://www.science.org/doi/pdf/10.1126/science.1233232 .
  * Abbott  _et al._ [2020b] R. Abbott _et al._ , GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object, Astrophys. J. 896, L44 (2020b).
  * Kiuchi  _et al._ [2017] K. Kiuchi, K. Kawaguchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Sub-radian-accuracy gravitational waveforms of coalescing binary neutron stars in numerical relativity, Phys. Rev. D 96, 084060 (2017).
  * Campanelli  _et al._ [2006] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Accurate evolutions of orbiting black-hole binaries without excision, Phys. Rev. Lett. 96, 111101 (2006).
  * Baker  _et al._ [2006] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, Gravitational-wave extraction from an inspiraling configuration of merging black holes, Phys. Rev. Lett. 96, 111102 (2006).
  * Marronetti  _et al._ [2008] P. Marronetti, W. Tichy, B. Brügmann, J. González, and U. Sperhake, High-spin binary black hole mergers, Phys. Rev. D 77, 064010 (2008).
  * Shibata and Nakamura [1995] M. Shibata and T. Nakamura, Evolution of three-dimensional gravitational waves: Harmonic slicing case, Phys. Rev. D 52, 5428 (1995).
  * Baumgarte and Shapiro [1998] T. W. Baumgarte and S. L. Shapiro, Numerical integration of einstein’s field equations, Phys. Rev. D 59, 024007 (1998).
  * Hilditch  _et al._ [2013] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy, and B. Brügmann, Compact binary evolutions with the z4c formulation, Phys. Rev. D 88, 084057 (2013).
  * Mignone and Bodo [2005] A. Mignone and G. Bodo, An HLLC Riemann solver for relativistic flows - I. Hydrodynamics, Monthly Notices of the Royal Astronomical Society 364, 126 (2005).
  * White  _et al._ [2016] C. J. White, J. M. Stone, and C. F. Gammie, An Extension of the Athena++ Code Framework for GRMHD Based on Advanced Riemann Solvers and Staggered-mesh Constrained Transport, Astrophys. J. Suppl. Ser. 225, 22 (2016).
  * Kiuchi  _et al._ [2022] K. Kiuchi, L. E. Held, Y. Sekiguchi, and M. Shibata, Implementation of advanced Riemann solvers in a neutrino-radiation magnetohydrodynamics code in numerical relativity and its application to a binary neutron star merger, Phys. Rev. D 106, 124041 (2022).
  * Papenfort  _et al._ [2021] L. J. Papenfort, S. D. Tootle, P. Grandclément, E. R. Most, and L. Rezzolla, New public code for initial data of unequal-mass, spinning compact-object binaries, Phys. Rev. D 104, 024057 (2021), arXiv:2103.09911 [gr-qc] .
  * Lattimer and Prakash [2004] J. M. Lattimer and M. Prakash, The physics of neutron stars, Science 304, 536 (2004).
  * Read  _et al._ [2009a] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Constraints on a phenomenologically parametrized neutron-star equation of state, Phys. Rev. D 79, 124032 (2009a).
  * Read  _et al._ [2009b] J. S. Read, C. Markakis, M. Shibata, K. Uryū, J. D. E. Creighton, and J. L. Friedman, Measuring the neutron star equation of state with gravitational wave observations, Phys. Rev. D 79, 124033 (2009b).
  * Lackey  _et al._ [2012] B. D. Lackey, K. Kyutoku, M. Shibata, P. R. Brady, and J. L. Friedman, Extracting equation of state parameters from black hole-neutron star mergers: Nonspinning black holes, Phys. Rev. D 85, 044061 (2012).
  * Hotokezaka  _et al._ [2013] K. Hotokezaka, K. Kiuchi, K. Kyutoku, H. Okawa, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Mass ejection from the merger of binary neutron stars, Phys. Rev. D 87, 024001 (2013).
  * Lousto  _et al._ [2010] C. O. Lousto, H. Nakano, Y. Zlochower, and M. Campanelli, Intermediate-mass-ratio black hole binaries: Intertwining numerical and perturbative techniques, Phys. Rev. D 82, 104057 (2010).
  * Reisswig and Pollney [2011] C. Reisswig and D. Pollney, Notes on the integration of numerical relativity waveforms, Classical Quantum Gravity 28, 195015 (2011).
  * Bardeen  _et al._ [1972] J. M. Bardeen, W. H. Press, and S. A. Teukolsky, Rotating black holes: Locally nonrotating frames, energy extraction, and scalar synchrotron radiation, Astrophys. J. 178, 347 (1972).
  * Akmal  _et al._ [1998] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Equation of state of nucleon matter and neutron star structure, Phys. Rev. C 58, 1804 (1998).
  * Glendenning and Moszkowski [1991] N. K. Glendenning and S. A. Moszkowski, Reconciliation of neutron-star masses and binding of the $\Lambda$ in hypernuclei, Phys. Rev. Lett. 67, 2414 (1991).
  * Lackey  _et al._ [2006] B. D. Lackey, M. Nayyar, and B. J. Owen, Observational constraints on hyperons in neutron stars, Phys. Rev. D 73, 024021 (2006).
  * Hinderer  _et al._ [2019b] T. Hinderer, S. Nissanke, F. Foucart, K. Hotokezaka, T. Vincent, M. Kasliwal, P. Schmidt, A. R. Williamson, D. A. Nichols, M. D. Duez, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Distinguishing the nature of comparable-mass neutron star binary systems with multimessenger observations: Gw170817 case study, Phys. Rev. D 100, 063021 (2019b).
  * Chawla  _et al._ [2010] S. Chawla, M. Anderson, M. Besselman, L. Lehner, S. L. Liebling, P. M. Motl, and D. Neilsen, Mergers of magnetized neutron stars with spinning black holes: Disruption, accretion, and fallback, Phys. Rev. Lett. 105, 111101 (2010).
  * Duez  _et al._ [2010] M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, and S. A. Teukolsky, Equation of state effects in black hole–neutron star mergers, Classical Quantum Gravity 27, 114106 (2010).
  * Etienne  _et al._ [2012a] Z. B. Etienne, Y. T. Liu, V. Paschalidis, and S. L. Shapiro, General relativistic simulations of black-hole–neutron-star mergers: Effects of magnetic fields, Phys. Rev. D 85, 064029 (2012a).
  * Etienne  _et al._ [2012b] Z. B. Etienne, V. Paschalidis, and S. L. Shapiro, General-relativistic simulations of black-hole–neutron-star mergers: Effects of tilted magnetic fields, Phys. Rev. D 86, 084026 (2012b).
  * Deaton  _et al._ [2013] M. B. Deaton, M. D. Duez, F. Foucart, E. O’Connor, C. D. Ott, L. E. Kidder, C. D. Muhlberger, M. A. Scheel, and B. Szilagyi, Black hole-neutron star mergers with a hot nuclear equation of state: Outflow and neutrino-cooled disk for a low-mass, high-spin case, Astrophys. J. 776, 47 (2013).
  * Paschalidis  _et al._ [2015] V. Paschalidis, M. Ruiz, and S. L. Shapiro, RELATIVISTIC SIMULATIONS OF BLACK HOLE–NEUTRON STAR COALESCENCE: THE JET EMERGES, Astrophys. J. 806, L14 (2015).
  * Foucart  _et al._ [2017] F. Foucart, D. Desai, W. Brege, M. D. Duez, D. Kasen, D. A. Hemberger, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Dynamical ejecta from precessing neutron star-black hole mergers with a hot, nuclear-theory based equation of state, Classical Quantum Gravity 34, 044002 (2017).
  * Kyutoku  _et al._ [2018] K. Kyutoku, K. Kiuchi, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Neutrino transport in black hole-neutron star binaries: Neutrino emission and dynamical mass ejection, Phys. Rev. D 97, 023009 (2018).
  * Ruiz  _et al._ [2018] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Jet launching from binary black hole-neutron star mergers: Dependence on black hole spin, binary mass ratio, and magnetic field orientation, Phys. Rev. D 98, 123017 (2018).
  * Foucart  _et al._ [2019b] F. Foucart, M. D. Duez, T. Hinderer, J. Caro, A. R. Williamson, M. Boyle, A. Buonanno, R. Haas, D. A. Hemberger, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Gravitational waveforms from spectral einstein code simulations: Neutron star-neutron star and low-mass black hole-neutron star binaries, Phys. Rev. D 99, 044008 (2019b).
  * Vallisneri [2000] M. Vallisneri, Prospects for gravitational-wave observations of neutron-star tidal disruption in neutron-star–black-hole binaries, Phys. Rev. Lett. 84, 3519 (2000).
  * Pannarale  _et al._ [2015] F. Pannarale, E. Berti, K. Kyutoku, B. D. Lackey, and M. Shibata, Gravitational-wave cutoff frequencies of tidally disruptive neutron star-black hole binary mergers, Phys. Rev. D 92, 081504 (2015).
  * Foucart  _et al._ [2013b] F. Foucart, L. Buchman, M. D. Duez, M. Grudich, L. E. Kidder, I. MacDonald, A. Mroue, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi, First direct comparison of nondisrupting neutron star-black hole and binary black hole merger simulations, Phys. Rev. D 88, 064017 (2013b).
  * Nakamura and Sasaki [1981] T. Nakamura and M. Sasaki, Is collapse of a deformed star always effectual for gravitational radiation?, Physics Letters B 106, 69 (1981).
  * Shapiro and Wasserman [1982] S. L. Shapiro and I. M. Wasserman, Gravitational radiation from nonspherical infall into black holes, The Astrophysical Journal 260, 838 (1982).
  * Nakamura and ichi Oohara [1983] T. Nakamura and K. ichi Oohara, Gravitational radiation emitted by n particles in circular orbits, Physics Letters A 98, 403 (1983).
  * Blanchet  _et al._ [1990] L. Blanchet, T. Damour, and G. Schäefer, Post-Newtonian hydrodynamics and post-Newtonian gravitational wave generation for numerical relativity, Monthly Notices of the Royal Astronomical Society 242, 289 (1990), https://academic.oup.com/mnras/article-pdf/242/3/289/3863050/mnras242-0289.pdf .
  * Kyutoku  _et al._ [2020b] K. Kyutoku, S. Fujibayashi, K. Hayashi, K. Kawaguchi, K. Kiuchi, M. Shibata, and M. Tanaka, On the Possibility of GW190425 Being a Black Hole–Neutron Star Binary Merger, Astrophys. J. Lett. 890, L4 (2020b), arXiv:2001.04474 [astro-ph.HE] .
  * Kawaguchi  _et al._ [2020] K. Kawaguchi, M. Shibata, and M. Tanaka, Constraint on the Ejecta Mass for Black Hole–Neutron Star Merger Event Candidate S190814bv, Astrophys. J. 893, 153 (2020), arXiv:2002.01662 [astro-ph.HE] .