# General-relativistic neutrino-radiation magnetohydrodynamics simulation of
seconds-long black hole-neutron star mergers: Dependence on initial magnetic
field strength, configuration, and neutron-star equation of state

Kota Hayashi1, Kenta Kiuchi2,1, Koutarou Kyutoku3,1,4, Yuichiro Sekiguchi5,1,
Masaru Shibata2,1 1Center for Gravitational Physics, Yukawa Institute for
Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
2Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, Postdam-Golm 14476, Germany
3Department of Physics, Kyoto University, Kyoto 606-8502, Japan
4Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
RIKEN, Wako, Saitama 351-0198, Japan
5Department of Physics, Toho University, Funabashi, Chiba 274-8510, Japan

###### Abstract

As a follow-up study of our previous work Hayashi  _et al._ (2022), numerical-
relativity simulations for seconds-long black hole-neutron star mergers are
performed for a variety of setups. Irrespective of the initial and symmetry
conditions, we find qualitatively universal evolution processes: The dynamical
mass ejection takes place together with a massive accretion disk formation
after the neutron star is tidally disrupted; Subsequently, the magnetic field
in the accretion disk is amplified by the magnetic winding, Kelvin-Helmholtz
instability, and magnetorotational instability, which establish a turbulent
state inducing the dynamo and angular momentum transport; The post-merger mass
ejection by the effective viscous processes stemming from the
magnetohydrodynamics turbulence sets in at $\sim 300$–$500$ ms after the
merger and continues for several hundred ms; A magnetosphere near the black-
hole spin axis is developed and the collimated strong Poynting flux is
generated with its lifetime of $\sim 0.5$–$2$ s.  We have newly found that the
model of no equatorial-plane symmetry shows the reverse of the magnetic-field
polarity in the magnetosphere, which is caused by the dynamo associated with
the magnetorotational instability in the accretion disk. The model with
initially toroidal fields shows the tilt of the disk and magnetosphere in the
late post-merger stage because of the anisotropic post-merger mass ejection.
These effects could terminate the strong Poynting-luminosity stage within the
timescale of $\sim 0.5$–$2$ s.

##  I Introduction

The first direct detection of gravitational waves from a binary black hole
merger, referred to as GW150914 Abbott  _et al._ (2016), opened the era of
gravitational-wave astronomy. To date, advanced LIGO and advanced Virgo have
observed $\sim 80$ binary black hole merger events Abbott  _et al._ (2021a,
b). Several neutron-star merger events have also been observed in addition to
the binary black hole mergers. A remarkable event is GW170817 Abbott  _et al._
(2017a), the first binary neutron star merger event. Associated with this
event, a wide variety of electromagnetic counterparts were successfully
observed Abbott  _et al._ (2017b, c), and provided us with invaluable
information for understanding the processes of the neutron-star merger and its
post-merger evolution. The electromagnetic observations also opened the era of
multi-messenger astronomy including gravitational-wave observation.

In the latest observational run (O3b), gravitational waves from black hole-
neutron star binaries, referred to as  (GW200105 and) GW200115 Abbott  _et
al._ (2021c), were observed. These events surely indicate that black hole-
neutron star binaries exist in nature. Although no electromagnetic counterpart
is observed for them, it is natural to expect that electromagnetic
counterparts will be observed in future events, if the binary parameters are
suitable for inducing tidal disruption of neutron stars. A number of
numerical-relativity simulations for black hole-neutron star binaries predict
that the neutron star could be disrupted by the tidal force of the companion
black hole if the black-hole mass is relatively small and/or the black-hole
spin is high (e.g., Refs. Shibata (2015); Kyutoku  _et al._ (2021)). The tidal
disruption is accompanied by disk formation and mass ejection, which will
result in the $r$-process nucleosynthesis for synthesizing heavy neutron-rich
elements Lattimer and Schramm (1974); Eichler  _et al._ (1989). If the remnant
black hole is rapidly spinning and surrounded by a massive magnetized disk or
torus, an ultra-relativistic jet could be launched and drive a short-hard
gamma-ray burst Eichler  _et al._ (1989); Nakar (2007); Berger (2014). Powered
by thermal energy generated by the radioactive decay of synthesized heavy
neutron-rich elements, the ejecta will shine with high luminosity as kilonovae
Li and Paczyński (1998); Metzger  _et al._ (2010).

The sensitivity of the gravitational-wave detectors is being improved for the
forthcoming observational runs (O4 and O5) Abbott  _et al._ (2020). Also,
large-scale telescopes such as JWST and Vera Rubin telescopes will be in
operation during such runs Gardner  _et al._ (2006); Ivezić  _et al._ (2019).
It is quite natural to expect simultaneous detection of gravitational waves
and electromagnetic counterparts from black hole-neutron star mergers if the
source is within a distance of several hundred Mpc from the earth. This
implies that black hole-neutron star mergers are among the most promising
sources for multi-messenger astronomy in the near future. In view of this
situation, it is urgent to theoretically develop the entire evolution scenario
from the merger to the post-merger stages, in order to predict observable
signals and to make a reliable model for the interpretation of the forthcoming
observational data.

In the last two decades, a variety of numerical-relativity simulations for
black hole-neutron star mergers have been performed Shibata and Uryū (2006,
2007); Shibata and Taniguchi (2008); Etienne  _et al._ (2008); Duez  _et al._
(2008); Shibata  _et al._ (2009); Etienne  _et al._ (2009); Chawla  _et al._
(2010); Duez  _et al._ (2010); Kyutoku  _et al._ (2010, 2011); Foucart  _et
al._ (2011, 2012); Etienne  _et al._ (2012a, b); Kyutoku  _et al._ (2013,
2015); Foucart  _et al._ (2013); Lovelace  _et al._ (2013); Deaton  _et al._
(2013); Foucart  _et al._ (2014); Paschalidis  _et al._ (2015); Kawaguchi  _et
al._ (2015); Kiuchi  _et al._ (2015); Foucart  _et al._ (2017); Kyutoku  _et
al._ (2018); Brege  _et al._ (2018); Ruiz  _et al._ (2018); Foucart  _et al._
(2019a, b); Hinderer  _et al._ (2019); Hayashi  _et al._ (2021); Foucart  _et
al._ (2021); Most  _et al._ (2021); Chaurasia  _et al._ (2021); Most  _et al._
(2021); Hayashi  _et al._ (2022). By improving the input physics and grid
resolution, the previous studies have extensively explored the process of the
tidal disruption, accretion disk formation, dynamical mass ejection,
gravitational-wave emission, and neutrino emission. However, most of the
previous works have focused only on the evolution from the inspiral to early
post-merger stages; the evolution was followed at longest for a few hundred ms
after the merger. Hence, the long-term seconds-long evolution of the system
has not been explored deeply. In order to compensate for this deficiency,
i.e., to explore the entire post-merger evolution processes, many long-term
numerical simulations for black hole-accretion disk systems have also been
performed, including viscous hydrodynamics or magnetohydrodynamics effects
Fernández and Metzger (2013); Metzger and Fernández (2014); Just  _et al._
(2015); Fernández  _et al._ (2015, 2017); Siegel and Metzger (2018); Fernández
_et al._ (2019); Janiuk (2019); Christie  _et al._ (2019); Miller  _et al._
(2019); Fujibayashi  _et al._ (2020a, b); Li and Siegel (2021); Fernández  _et
al._ (2020); Just  _et al._ (2022); Shibata  _et al._ (2021). These
simulations have qualitatively clarified the evolution processes in the post-
merger stage such as the post-merger mass ejection and jet launch. For
example, it is now widely accepted that the post-merger mass ejection is
likely to be driven by the effective viscous processes induced by the
magnetohydrodynamics turbulence. The jet outflow is also likely to be powered
by the Blandford-Znajek mechanism Blandford and Znajek (1977) by magnetic
fields penetrating a rapidly spinning black hole.  However, it is not clear
whether the initial conditions given in such simulations are appropriate for
exploring the post-merger evolution of the black hole-neutron star merger.
During the merger stage, the neutron-star matter is spread non-
axisymmetrically around the remnant black hole, and the non-axisymmetric fall-
back matter suppresses the system from being axisymmetric. Such non-
axisymmetric structure of the accretion disk could also significantly affect
the magnetic-field amplification process and field profile. Therefore, the
quantitative details of the post-merger process are not fully understood.

In order to overcome these deficiencies and acquire the self-consistent
evolution scenario of the black hole-neutron star merger starting from the
inspiral stage to the late post-merger stage, in a previous paper, we
performed seconds-long merger simulations including dynamical general-
relativity effect, neutrino-radiation effect, and magnetohydrodynamics effect
altogether for the first time Hayashi  _et al._ (2022). We confirmed that the
post-merger mass ejection is indeed driven by the magnetically-induced viscous
effect at several hundred milliseconds after the onset of the merger, and this
post-merger mass ejection is triggered by the decrease of the temperature and
neutrino luminosity in the accretion disk, as previous viscous hydrodynamics
studies have clarified (e.g., Refs. Fernández and Metzger (2013); Fujibayashi
_et al._ (2020a); Just  _et al._ (2022)). We also confirmed that the electron
fraction of the post-merger ejecta is not very low, typically with
$Y_{\mathrm{e}}\sim 0.2$–0.3. In addition, we found the development of the
magnetosphere near the spin axis of the remnant black hole as a result of the
infall of the amplified magnetic flux from the disk into the black hole. The
intensity of the outgoing Poynting flux in the magnetosphere is high and is
consistent with that of short-hard gamma-ray bursts.

Although the results of our previous work showed a self-consistent evolution
picture of black hole-neutron star binaries, some questions remain to be
answered. First, our previous work assumed initially strong magnetic fields
with its maximum strength $\geq 3\times 10^{16}$ G. The high magnetic-field
strength was given in order to get a high field strength at the formation of
the accretion disk and to enable numerical computation to resolve the
magnetorotational instability (MRI) Balbus and Hawley (1991, 1998) in the disk
right after the disk formation. We should clarify whether, at least
qualitatively, the same evolution is obtained even if we assume a much lower
strength of the initial magnetic field. Second, our previous work assumed a
poloidal magnetic field confined in the neutron star as an initial condition.
It is quite natural to ask whether the result is qualitatively changed or not
if we employ different magnetic-field configuration in the neutron star as the
initial condition. Third, in our previous work, we imposed the equatorial-
plane symmetry on all the simulations. It is necessary to understand what
happens in the absence of such a symmetry. Finally, our previous work used
only the DD2 equation of state (EOS) Banik  _et al._ (2014) to model the
neutron star. We should investigate how the evolution process  and properties
of the ejecta depend quantitatively on the EOSs employed.

The paper is organized as follows. In Sec. II, we briefly describe the method
and initial setup for the numerical simulation. In Sec. III, we present the
numerical results focusing on the entire evolution process, mass ejection
mechanisms, and collimated electromagnetic outflow developed near the spin
axis of the black hole. We pay particular attention to the quantitative
difference among the numerical results with different initial  magnetic-field
profile, different EOSs, and different computational setups. Finally, we
conclude this work in Sec. IV. Throughout this paper, we use the geometrical
units in which $G=c=1$, where $G$ and $c$ are the gravitational constant and
the speed of light, respectively.

##  II Numerical methods

Figure 1:  The white lines show the magnetic-field line for the two types of
the initial magnetic field configurations adopted in this work. The magnetic
field is confined in the neutron star initially. The top and bottom panels
show the cases for the poloidal and toroidal fields, respectively.

The numerical scheme for the present simulation is the same as that summarized
in Ref. Hayashi  _et al._ (2022) for which the readers may refer to.

For modeling the neutron-star matter, we employ two nuclear-theory-based
finite-temperature EOSs referred to as DD2 Banik  _et al._ (2014) and SFHo
Steiner  _et al._ (2013) for a high-density range and Helmholtz EOS Timmes and
Swesty (2000) for a low-density range. As initial data, we prepare black hole-
neutron star binaries in a quasi-equilibrium state assuming the neutrinoless
beta-equilibrium cold state Kyutoku  _et al._ (2018). The initial
gravitational mass of the neutron star is set to be $M_{\rm NS}=1.35M_{\odot}$
following Ref. Hayashi  _et al._ (2022). For the DD2 and SFHo EOSs, the
circumferential radius of the isolated spherical neutron star of mass
$1.3$–$1.4M_{\odot}$ is $\approx 13.2\,\mathrm{km}$ and $\approx
11.9\,\mathrm{km}$, respectively. These EOSs satisfy constraints imposed by
the observation of gravitational waves for GW170817 Abbott  _et al._ (2017a)
and by the X-ray observation by NICER Miller  _et al._ (2019).

For the initial black-hole mass, we choose $M_{\rm BH,0}=5.4M_{\odot}$; the
mass ratio of the black hole to the neutron star is $Q:=M_{\rm BH,0}/M_{\rm
NS}=4$. The initial dimensionless spin parameter of the black hole is set to
be $\chi=0.75$. This setting is the same as one of our previous settings
Hayashi  _et al._ (2022). With such a spin and a mass ratio, tidal disruption
of the neutron star with $M_{\rm NS}=1.35M_{\odot}$ takes place for both EOSs.
The initial orbital angular velocity $\Omega_{0}$ is set to be
$m_{0}\Omega_{0}=0.056$ for $Q=4$ where $m_{0}=M_{\rm BH,0}+M_{\rm
NS}=6.75M_{\odot}$. With this initial setup, the binary spends about three
orbits before the merger.

We initially superimpose a poloidal or toroidal magnetic field confined in the
neutron star.  Figure 1 shows the magnetic-field lines for two types of
initial field configurations adopted in this work. For the poloidal field
case, following our previous work Kiuchi  _et al._ (2015), the magnetic field
is given in terms of the vector potential as

$\displaystyle A_{j}$ $\displaystyle=$
$\displaystyle\\{-y_{\mathrm{NS}}\delta_{j}^{~{}x}+x_{\mathrm{NS}}\delta_{j}^{~{}y}\\}$
(1) $\displaystyle\times A_{b}\max(P/P_{\mathrm{max}}-10^{-3},0)^{2}.$

For the toroidal field case, the vector potential is given as

$\displaystyle A_{j}$ $\displaystyle=$
$\displaystyle\\{(x_{\mathrm{NS}}(z_{\mathrm{NS}}^{2}-R_{\mathrm{NS}}^{2}))\delta_{j}^{~{}x}+(y_{\mathrm{NS}}(z_{\mathrm{NS}}^{2}-R_{\mathrm{NS}}^{2}))\delta_{j}^{~{}y}$
(2)
$\displaystyle-(z_{\mathrm{NS}}(x_{\mathrm{NS}}^{2}+y_{\mathrm{NS}}^{2}-R_{\mathrm{NS}}^{2}))\delta_{j}^{~{}z}\\}$
$\displaystyle\times A_{b}(1+\cos(r_{\mathrm{NS}}/0.95R_{\mathrm{NS}}))$
$\displaystyle\qquad\qquad(r_{\mathrm{NS}}<0.95R_{\mathrm{NS}}).$

Here, $(x_{\mathrm{NS}},y_{\mathrm{NS}},z_{\mathrm{NS}})$ denote the
coordinates with respect to the neutron-star center (location of the maximum
rest-mass density), $r_{\mathrm{NS}}$ is the  radial coordinate with respect
to the neutron-star center, and $R_{\mathrm{NS}}$ is the coordinate radius of
the neutron-star. $P$ is the pressure, $P_{\mathrm{max}}$ is the maximum
pressure, and $j=x,y$, and $z$. $A_{b}$ is a constant and is chosen so that
the initial maximum magnetic-field strength $b_{\mathrm{0,max}}$ is $3\times
10^{15}~{}\mathrm{G}$ or $5\times 10^{16}~{}\mathrm{G}$. These values are
chosen to obtain a strong magnetic field in the remnant disk formed shortly
after the tidal disruption of the neutron star. The strong magnetic field in
the remnant disk is required to resolve the fastest growing mode of the MRI
Balbus and Hawley (1991, 1998) in the limited grid resolution, because its
wavelength is proportional to the magnetic-field strength. Although strong
fields we prepare are not realistic for orbiting neutron stars, the resulting
turbulent state in the accretion disk established by the MRI is not likely to
depend strongly on the initial magnetic-field strength. Thus, it would be
reasonable to suppose that the turbulent state with a strong magnetic field
will be established even for the case with much weaker initial magnetic-field
strength if the grid resolution is sufficient.  As is often done in this
research field, we implicitly assume that the magnetic-field strength would be
increased by the MRI and magnetic winding, and a turbulent state would be
eventually established even if we started a simulation from low magnetic-field
strengths. This is just an assumption, but the result of a simulation, which
is started with a low magnetic-field strength, indicates certain evidence for
this (see Sec. III). We also note that even with $b_{\mathrm{0,max}}=5\times
10^{16}$ G, the electromagnetic energy (of order $10^{49}$ erg) is much
smaller than the internal energy and gravitational potential energy (of order
$10^{53}$ erg) of the neutron star, and thus, the inspiral and tidal-
disruption stages are not affected by the strong field significantly.

We only consider the magnetic field confined in the neutron star initially,
and do not consider a pulsar-like dipole magnetic field extending to the
outside of the neutron star. This is because only the magnetic field confined
in the neutron star has a significant effect on the subsequent evolution of
the system  for the realistic magnetic-field strength. In terms of the
accretion disk evolution including the post-merger mass ejection, only the
magnetic field in the disk, which originates from the magnetic field inside
the neutron star, plays an important role. In terms of the magnetosphere
formation, a dipolar magnetic field initially located outside the neutron star
may be amplified linearly due to winding. However in the disk, the magnetic
field is amplified exponentially by the MRI, and the amplified magnetic field
flux is ejected from the disk to the polar region by the MRI dynamo and
subsequently forms the magnetosphere of a high field strength. For the
realistic initial magnetic-field strength lower than $10^{12}$ G, the magnetic
field amplified by the MRI should dominantly come into play.

Table 1:  Key parameters and quantities for the initial conditions together with the parameters of grid setup for our numerical simulations. $b_{\rm 0,max}$: the initial maximum magnetic-field strength, $\Delta x_{i_{\rm max}}$: the grid spacing for the finest refinement level, $L_{1}$: the location of the outer boundaries along each axis, and the values of $N$ and $i_{\rm max}$. For all the models, the neutron-star mass is $1.35M_{\odot}$, the initial black-hole mass is $5.4M_{\odot}$, the initial dimensionless spin of the black hole is 0.75, and the initial ADM mass $M_{\rm ADM,0}$ is $6.679M_{\odot}\approx 0.9894m_{0}$. The models from our previous paper (Q4B5L and Q4B5H) are also shown for comparison.  model name | EOS | $b_{\rm 0,max}~{}[{\rm G}]$ |  $b_{\rm 0}$ configuration | plane sym. | $\Delta x_{i_{\rm max}}~{}[{\rm m}]$ | $L_{1}~{}[{\rm km}]$ |  $N$ |  $i_{\rm max}$
---|---|---|---|---|---|---|---|---
Q4B3e15 | DD2 | $3\times 10^{15}$ | poloidal | yes | 400 | $6.98\times 10^{4}$ | 170 | 11
Q4B5tn | DD2 | $5\times 10^{16}$ | toroidal | no | 400 | $6.98\times 10^{4}$ | 170 | 11
Q4B5n | DD2 | $5\times 10^{16}$ | poloidal | no | 400 | $6.98\times 10^{4}$ | 170 | 11
SFHoQ4B5 | SFHo | $5\times 10^{16}$ | poloidal | yes | 250 | $3.10\times 10^{4}$ | 243 | 10
Q4B5L Hayashi  _et al._ (2022) | DD2 | $5\times 10^{16}$ | poloidal | yes | 400 | $1.74\times 10^{4}$ | 170 | 9
Q4B5H Hayashi  _et al._ (2022) | DD2 | $5\times 10^{16}$ | poloidal | yes | 270 | $1.62\times 10^{4}$ | 234 | 9

We do not consider the effect of the neutrino viscosity on the MRI supposing
that the magnetic-field strength could be enhanced to be $\gtrsim 10^{14}$ G
due to the rapid winding in the main region of the accretion disk (see Sec.
III.2) even if the early growth of the MRI could be suppressed Masada and Sano
(2008); Guilet  _et al._ (2017).

The simulation is performed using a fixed-mesh refinement (FMR) algorithm. The
$i$-th refinement level covers a half or full cubic box of
$[-L_{i}:L_{i}]\times[-L_{i}:L_{i}]\times[0:L_{i}]$ or
$[-L_{i}:L_{i}]\times[-L_{i}:L_{i}]\times[-L_{i}:L_{i}]$, where $L_{i}=N\Delta
x_{i}$ and $\Delta x_{i}$ is the grid spacing for the $i$-th level. For the
half-cubic box case, the plane-symmetric boundary condition on the $z=0$ plane
(equatorial plane) is imposed. The grid spacing for each level is determined
by $\Delta x_{i}=2\Delta x_{i+1}$ ($i=1,2,\cdots,i_{\mathrm{max}}-1$) with
$\Delta x_{i_{\mathrm{max}}}=400\,\mathrm{m}$ for the DD2 models and $\Delta
x_{i_{\mathrm{max}}}=250\,\mathrm{m}$ for the SFHo model. That is, we perform
lower-resolution simulations for the DD2 EOS models, because high-resolution
simulations require an extremely high computational cost for simulating the
seconds-long merger processes. In addition, as we showed in our previous paper
Hayashi  _et al._ (2022), the results for the lower-resolution runs are
quantitatively similar to those for the corresponding high-resolution runs
(with $\Delta x_{i_{\mathrm{max}}}=270\,\mathrm{m}$), and hence, we consider
that a fair convergence would be achieved even with the present choice.
$i_{\rm max}$ is chosen to be 11 for the DD2 models and 10 for the SFHo model.
The values of $N$ are $170$ for the DD2 models and 243 for the SFHo model,
respectively (cf. Table 1).

In this paper, we perform four new simulations varying the EOS, the value of
$b_{\mathrm{0,max}}$, the magnetic-field configuration, and the equatorial-
plane symmetry. The parameters and quantities for the four models are
summarized in Table 1. We compare the new results with those obtained in our
previous paper Hayashi  _et al._ (2022) for the models with the same values of
$M_{\rm BH,0}$, $M_{\mathrm{NS}}$, and $\chi$.

As we already mentioned in our previous paper Hayashi  _et al._ (2022), during
the merger stage, the black hole is kicked mainly by the back reaction of the
dynamical mass ejection and the resulting velocity is $v_{\rm kick}=200$–400
km/s in our present setting. For the SFHo model for which the ejecta mass is
smaller, the kick velocity is lower. In order to avoid the black hole from
running into the FMR boundary, we control the shift vector with the
prescription proposed in our previous paper Hayashi  _et al._ (2022).

Following our previous work Hayashi  _et al._ (2022), we stop the time
evolution of the gravitational field at a certain moment after the ratio of
the rest mass of the remnant disk to the black-hole mass drops below
$10^{-2}$. This prescription is reasonable because the self-gravity of the
matter located outside the black hole can be safely neglected and the
gravitational field is approximately stationary in such a low-mass disk stage.
For stabilizing magnetohydrodynamics computation, we introduce floor density
of $\rho=10^{3}\,{\rm g/cm}^{3}$ for $r\leq r_{0}\approx 10^{2}$ km and
max$(10^{3}(r_{0}/r)^{3},0.17)\,{\rm g/cm}^{3}$ for $r>r_{0}$ where
$0.17\,{\rm g/cm^{3}}$ is the lowest value of the rest-mass density in our EOS
table.

##  III Results

###  III.1 Overview

Table 2:  The first two columns list $M_{\mathrm{eje,dyn}}$: the dynamical ejecta mass evaluated at $t=20$ ms, and $M_{\mathrm{eje,pm}}$: the lower bound of the post-merger ejecta mass in units of $M_{\odot}$. Since the mass of the post-merger ejecta is still increasing at the termination of all the runs, we here list the lower bound for it (see footnote 1). The next two columns list $Y_{e\ \mathrm{eje,dyn}}$: the averaged electron fraction of the dynamical ejecta evaluated at $t=20$ ms, and $Y_{e\ \mathrm{eje,pm}}$: that of the post-merger ejecta. The last two columns list $v_{\mathrm{eje,dyn}}$: the averaged velocity of the dynamical ejecta evaluated at $t=20$ ms, and $v_{\mathrm{eje,pm}}$: that of the post-merger ejecta. model name | $M_{\mathrm{eje,dyn}}$ | $M_{\mathrm{eje,pm}}$ | $Y_{e\ \mathrm{eje,dyn}}$ | $Y_{e\ \mathrm{eje,pm}}$ | $v_{\mathrm{eje,dyn}}$ | $v_{\mathrm{eje,pm}}$
---|---|---|---|---|---|---
Q4B3e15 | 0.045 | 0.030 | 0.062 | 0.24 | 0.17 | 0.044
Q4B5tn | 0.045 | 0.030 | 0.062 | 0.21 | 0.17 | 0.051
Q4B5n | 0.046 | 0.033 | 0.062 | 0.22 | 0.17 | 0.053
SFHoQ4B5 | 0.013 | 0.019 | 0.051 | 0.27 | 0.16 | 0.048
Q4B5L Hayashi  _et al._ (2022) | 0.046 | 0.035 | 0.061 | 0.23 | 0.17 | 0.051
Q4B5H Hayashi  _et al._ (2022) | 0.046 | 0.028 | 0.056 | 0.25 | 0.17 | 0.047

Figure 2: The three-dimensional snapshots for model Q4B5tn in a domain with
the length scale of $\sim 10^{4}$ km at $t\approx 150$, $1200$, $1700$, and
$2900$ ms. For each time slice, the left panel shows the ejecta, which is
colored for the electron fraction $Y_{\mathrm{e}}$; the middle panel shows the
rest-mass density $\rho$ $\mathrm{(g/cm^{3})}$ (contours) with magnetic-field
lines (pink lines), unbound outflow (white color) and its velocity (green
arrows); the right panel shows the magnetic-field strength
$b=\sqrt{b^{\mu}b_{\mu}}$ $\mathrm{(G)}$. See also the following link for the
animation:
%****␣main.tex␣Line␣275␣****https://www2.yukawa.kyoto-u.ac.jp/~kota.hayashi/Q4B5tn-3D.mp4.

Figure 3: The 2D snapshots for model Q4B5tn on the $y$-$z$ plane with a region
of $[-2000~{}{\rm km}:2000~{}{\rm km}]$ for both $y$ and $z$ at $t\approx
150$, $350$, $550$, $800$, and $1700$ ms. For each time slice, the first,
second, third, and fourth panels show the rest-mass density $\rho$
$\mathrm{(g/cm^{3})}$, the electron fraction $Y_{\mathrm{e}}$, the entropy per
baryon $s$ in units of $k$, and $-hu_{t}-h_{\mathrm{min}}$, respectively. In
the fourth panel, unbound matter is (non-black) colored and bound matter is
colored by black. See also the following link for the animation:
https://www2.yukawa.kyoto-u.ac.jp/~kota.hayashi/Q4B5tn-2Dyz.mp4.

Figure 4: The 2D snapshots for model Q4B5tn on the $x$-$y$ plane with
$[-200~{}{\rm km}:200~{}{\rm km}]$ for both $x$ and $y$ at $t\approx 30$ ms,
and with $[-2000~{}{\rm km}:2000~{}{\rm km}]$ at $t\approx 250$ and $800$ ms.
For each time slice, the first, second, third, and fourth panels show the
rest-mass density $\rho$ $\mathrm{(g/cm^{3})}$, the electron fraction
$Y_{\mathrm{e}}$, the entropy per baryon $s$ in units of $k$, and the inverse
of the plasma beta $b^{2}/8\pi P$, respectively. See also the following link
for the animation:
https://www2.yukawa.kyoto-u.ac.jp/~kota.hayashi/Q4B5tn-2Dxy.mp4.

First, we show that the overall evolution process in the merger and post-
merger stages are qualitatively the same as our previous findings in Ref.
Hayashi  _et al._ (2022) irrespective of initial magnetic field strengths,
configurations, neutron-star EOSs, and equatorial-plane symmetry. Table 2
shows a summary of the outcome quantities such as the mass, the averaged
electron fraction and the velocity of the dynamical and post-merger ejecta.

Figure 2 displays the three-dimensional snapshots for model Q4B5tn in a domain
with the length scale of $\sim 10^{4}$ km at $t\approx 150$, $1200$, $1700$,
and $2900$ ms. For this model, we initially provide a toroidal magnetic field
with its maximum strength of $5\times 10^{16}$ G, and do not impose the
equatorial-plane symmetry. For each time slice, the left panel shows the
ejecta, which is colored for the electron fraction $Y_{\mathrm{e}}$; the
middle panel shows the rest-mass density $\rho$ $\mathrm{(g/cm^{3})}$
(contours) with magnetic-field lines (pink lines), unbound outflow (white
color) and its velocity (green arrows); the right panel shows the magnetic-
field strength defined by $b=\sqrt{b^{\mu}b_{\mu}}$ $\mathrm{(G)}$ where
$b^{\mu}$ is the magnetic field in the frame comoving with fluid (see Eq. (10)
for the definition of the energy-momentum tensor in our work). Figure 3
displays the two-dimensional (2D) snapshots for the same model (Q4B5tn) on the
$y$-$z$ plane with a region of $[-2000~{}{\rm km}:2000~{}{\rm km}]$ for both
$y$ and $z$ at $t\approx 150$, $350$, $550$, $800$, and $1700$ ms. For each
time slice, the first, second, third, and fourth panels show the rest-mass
density, the electron fraction, the entropy per baryon $s$ in units of the
Boltzmann constant $k$, and $-hu_{t}-h_{\mathrm{min}}$, respectively. Here,
$u_{t}$ is the lower time component of the four velocity, $h$ is the specific
enthalpy, and $h_{\mathrm{min}}$ is the minimum specific enthalpy for a given
electron fraction. In the fourth panel, unbound matter, which is identified by
$-hu_{t}-h_{\mathrm{min}}>0$, is (non-black) colored and bound matter is
colored by black. Figure 4 displays the 2D snapshots for the same model
(Q4B5tn) on the $x$-$y$ plane at $t\approx 30$, $250$, and $800$ ms. For each
time slice, the first, second, third, and fourth panels show the rest-mass
density, the electron fraction, the entropy per baryon, and the inverse of the
plasma beta, $b^{2}/8\pi P$, respectively.

These figures show that the merger and post-merger processes are qualitatively
quite similar to those found in our previous paper Hayashi  _et al._ (2022).
The neutron star is tidally disrupted but about 80% of the neutron-star matter
falls immediately into the black hole. On the other hand, a part of the
neutron-star matter becomes unbound from the system to be the dynamical ejecta
and the rest of the matter forms an accretion disk around the remnant black
hole. In the accretion disk which has shear layer due to the presence of the
differential rotation, the magnetic field is amplified by the winding, MRI,
and Kelvin-Helmholtz instability. Here the Kelvin-Helmholtz instability plays
an important role for the region in which spiral arms collide with each other
(see, e.g., the first snapshots of Fig. 4). After a sufficient amplification
of the magnetic field, a turbulent state is developed, and as a result, the
effective viscosity associated with the turbulent state is established. Then,
the angular momentum is efficiently transported outward, and the density and
temperature of the disk decrease with time. In an early stage of the disk
evolution, the temperature is so high that the neutrino cooling is efficient,
and hence, the thermal energy generated by the effective viscous processes is
liberated by the neutrino emission. However, after the disk temperature drops,
the neutrino cooling is inefficient, and thus, the viscous heating can be used
to eject a part of the matter from the disk, inducing the post-merger mass
ejection. After the amplification of the magnetic field in the disk, the
magnetic field along the spin axis of the black hole is also enhanced (see
Fig. 2) forming a magnetosphere along the spin axis of the black hole, in
which the density is low (see, e.g., Fig. 3) and the magnetization parameter
$b^{2}/4\pi\rho$ is above unity.

Although the overall qualitative evolution process is mostly independent of
the initial and numerical setups, we find some differences in the post-merger
mass ejection and the evolution of the magnetosphere. In the absence of the
equatorial-plane symmetry, the distribution of the post-merger ejecta is not
equatorially symmetric (see, e.g., the third and fourth snapshots of Fig. 3).
We also observe that the accretion disk and the magnetosphere tilt with time
in the late evolution stage of $t>2$ s (see, e.g., the fifth snapshots of Fig.
3). In addition, another model with no equatorial-plane symmetry (model Q4B5n)
shows that magnetic-field polarity in the magnetosphere reverses (see Sec.
III.3.1). These differences are likely due to the fact that the post-merger
evolution is determined by the stochastic turbulent processes in the disk. The
properties of the magnetosphere mentioned here play a role in the decreasing
stage of the isotropic-equivalent Poynting luminosity (see Sec. III.3.3).

In the following, thus, we present the results for a more detailed analysis of
the system paying attention to the quantitative dependence on the numerical
setting. In addition, the anisotropic part of the Maxwell and the Reynolds
stresses are carefully evaluated to discuss the qualitative aspect of the
effective viscosity induced by the magnetohydrodynamical turbulence (see Sec.
III.2.3). We also propose a new method to assess whether the black hole has
the ability to form a magnetosphere and launch a jet by evaluating
magnetohydrodynamics properties near the horizon (see Sec. III.3.2).

###  III.2 The evolution of the accretion disk and post-merger mass ejection

Figure 5: The time evolution of the rest mass of the matter located outside
the apparent horizon (dotted curves) and the accretion-disk mass (solid
curves) for models Q4B3e15 (top panel), Q4B5n and Q4B5tn (middle panel), and
SFHoQ4B5 (bottom panel). The results for models Q4B5L and Q4B5H of our
previous paper Hayashi  _et al._ (2022) are also shown for comparison (in grey
color).

Figure 6: The same as Fig. 5 but for the time evolution of the rest mass of
the unbound matter.

Figure 7: The same as Fig. 5 but for the time evolution of the total neutrino
luminosity. The post-merger mass ejection sets in at $t\sim 300$–600 ms at
which $L_{\nu}\sim 10^{51}$–$10^{51.5}~{}{\rm erg/s}$.

####  III.2.1 Disk evolution and ejecta

In this subsection, we present details on the evolution of the accretion disk
and on the post-merger mass ejection  focusing particularly on the
quantitative dependence on the models. Figure 5 shows the time evolution of
the rest mass of the matter located outside the apparent horizon
$M_{>\mathrm{AH}}$ (dashed curves) and the accretion disk mass $M_{\rm disk}$
(solid curves). Figure 6 shows the time evolution of the rest mass of the
unbound matter (ejecta) $M_{\mathrm{eje}}$. For both figures, the new results
obtained in this paper are compared with the previous results of Ref. Hayashi
_et al._ (2022).

The quantities shown in Figs. 5 and 6 are defined by

$\displaystyle M_{>{\rm AH}}$ $\displaystyle:=$ $\displaystyle\int_{r>r_{\rm
AH}}\rho_{*}d^{3}x+M_{\rm esc},$ (3) $\displaystyle M_{\rm eje}$
$\displaystyle:=$ $\displaystyle\int_{-hu_{t}>h_{\rm min},r>r_{\rm
AH}}\rho_{*}d^{3}x+M_{\rm esc},$ (4) $\displaystyle M_{\rm disk}$
$\displaystyle:=$ $\displaystyle M_{>{\rm AH}}-M_{\rm eje},$ (5)

where $\rho_{*}:=\rho\sqrt{-g}u^{t}$ with $g$ the determinant of the spacetime
metric, $g_{\mu\nu}$, $u^{t}$ the time component of the four velocity,
$u^{\mu}$, and $r_{\rm AH}$ denotes the coordinate radius of the apparent
horizon with respect to the position of the black-hole puncture. $M_{\rm esc}$
denotes the rest mass of the matter escaping from the computational domain
calculated by

$\displaystyle\dot{M}_{\rm esc}$ $\displaystyle:=$
$\displaystyle\oint\rho\sqrt{-g}u^{i}dS_{i},$ (6) $\displaystyle M_{\rm esc}$
$\displaystyle:=$ $\displaystyle\int^{t}\dot{M}_{\rm esc}dt.$ (7)

The surface integral is performed at the outer boundaries of the computational
domain. The ejecta component is identified by considering the Bernoulli
criterion; we identify the matter located outside the apparent horizon that
satisfies $hu_{t}<-h_{\rm min}$ as the unbound component.

Our new results for $M_{>\mathrm{AH}}$, $M_{\mathrm{disk}}$, and
$M_{\mathrm{eje}}$ show that the dependence of these quantities on the initial
magnetic field setup and equatorial symmetry imposed is weak. Also, the
results for the SFHo EOS model are qualitatively similar to those for the DD2
model, although for the SFHo model, $M_{>\mathrm{AH}}$, $M_{\mathrm{disk}}$,
and $M_{\mathrm{eje}}$ are smaller than those for the DD2 model reflecting the
smaller neutron-star radius for the SFHo model.

$M_{>\mathrm{AH}}$ decreases steeply at $\sim 10$ ms at which the merger
occurs and $\sim 80$–$85\%$ of the neutron-star matter plunges into the black
hole. After that, $M_{>\mathrm{AH}}$ continues to decrease gradually due to
the matter accretion into the black hole induced by the angular-momentum
transport from the magnetohydrodynamics effect. $M_{\mathrm{disk}}$ right
after the merger is $M_{\mathrm{disk,0}}\approx 0.28M_{\odot}$ and $\approx
0.22M_{\odot}$ for the DD2 and SFHo models, respectively.

The evolution of $M_{\mathrm{eje}}$ clearly shows that two distinct components
of the ejecta exist. One is the dynamical ejecta, for which $M_{\mathrm{eje}}$
steeply increases right after the merger spending only for a few ms. The rest
mass for this component is $\approx 0.046M_{\odot}$ and $\approx
0.013M_{\odot}$ for DD2 and SFHo models, respectively. After this increase by
the dynamical mass ejection, $M_{\mathrm{eje}}$ remains approximately constant
for several hundred ms. Then, $M_{\mathrm{eje}}$ starts increasing again at
$t\sim 300$–$600$ ms. This component is the post-merger ejecta driven by the
heating associated with MRI turbulence after the neutrino luminosity decreases
below the heating rate (see Fig. 7), i.e., $L_{\nu}$ decreases below
$10^{51}$–$10^{51.5}$ erg/s. The rest mass for this component is $\approx
0.030M_{\odot}$ and $\approx 0.019M_{\odot}$ for the DD2 and SFHo models at
the termination of the simulations, respectively. 111Strictly speaking we can
provide only the lower bound of the post-merger ejecta mass because at the
termination of the simulations, the mass is still increasing slightly. Note,
however, that the possible additional increase is less than $0.01M_{\odot}$
because the disk mass at the termination of the simulation is less than
$0.01M_{\odot}$. These values are about 10% of $M_{\mathrm{disk,0}}$. The
result for the DD2 model shows good agreement with our previous results
Hayashi  _et al._ (2022).

Compared with the simulations starting from a black hole-disk system such as
Refs. Siegel and Metzger (2018); Fernández  _et al._ (2019); Christie  _et
al._ (2019), the ratio of the post-merger ejecta mass to the remnant disk mass
for the present simulations is smaller by a factor of several. One reason for
this is that the early-post-merger mass ejection (which could occur within
100–200 ms after the onset of the merger) is not found in the present
simulations. This is because the strong poloidal magnetic field, which is the
key for this, is not developed in the disk remaining after the merger. Another
reason is that the mass accretion onto the black hole is enhanced due to the
highly non-axisymmetric structure of the remnant disk from black hole-neutron
star mergers.

The only significant quantitative dependence on the computational setting is
found on the onset time of the post-merger mass ejection. For example, for
model Q4B3e15 which has the low initial-magnetic-field strength, the onset
time of the post-merger mass ejection is $t\sim 600$ ms, i.e., $\sim 200$ ms
behind a high initial-magnetic-field strength model Q4B5L. The reason for this
is that for the model with the low initial field strength it takes a longer
time until the magnetic-field strength is enhanced enough for the disk to be
in the equipartition state and for numerical computation to resolve the
fastest growing mode of the MRI. It results in the delay of the development of
the MRI-induced turbulence in the accretion disk. Note, however, that this
delay  is likely due to the insufficient grid resolution and may not be
present in the realistic case, in which the MRI is resolved even for the weak
fields. Besides this difference,  the quantities of the post-merger mass
ejection such as the ejecta mass (see Table 2) are similar to those found in
our previous paper Hayashi  _et al._ (2022).

####  III.2.2 Magnetic-field evolution

Figure 8: The time evolution of the electromagnetic energy evaluated for
models Q4B3e15 (top panel), Q4B5n and Q4B5tn (middle panel), and SFHoQ4B5
(bottom panel). The results for models Q4B5L and Q4B5H of our previous paper
Hayashi  _et al._ (2022) are also shown for comparison (in grey color).

Figure 9: The same as Fig. 8 but for the time evolution of the ratio of the
electromagnetic energy to the internal energy evaluated for the outside of the
apparent horizon.

Figures 8 and 9 show the time evolution of the electromagnetic energy, $E_{\rm
B}$, and the ratio of the electromagnetic energy to the internal energy,
$E_{\rm int}$, respectively. Here, $E_{\rm B}$ and $E_{\rm int}$ are defined,
respectively, by

$\displaystyle E_{\rm B}$ $\displaystyle:=$
$\displaystyle\frac{1}{8\pi}\int_{r>r_{\rm
AH}}u^{t}\sqrt{-g}\,b_{\mu}b^{\mu}d^{3}x,$ (8) $\displaystyle E_{\rm int}$
$\displaystyle:=$ $\displaystyle\int_{r>r_{\rm AH}}\rho_{*}\varepsilon
d^{3}x,$ (9)

where $\varepsilon$ denotes the specific internal energy. In this paper the
energy-momentum tensor for the ideal magnetohydrodynamics is written as

$\displaystyle T_{\mu\nu}$ $\displaystyle=$ $\displaystyle\rho
hu_{\mu}u_{\nu}+Pg_{\mu\nu}$ (10) $\displaystyle+$ $\displaystyle{1\over
4\pi}\left(b^{\alpha}b_{\alpha}u_{\mu}u_{\nu}+{1\over
2}b^{\alpha}b_{\alpha}g_{\mu\nu}-b_{\mu}b_{\nu}\right).$

A quantitative difference is found for models Q4B3e15 and Q4B5nt at $t\lesssim
100$ ms, which corresponds to the stage where the disk is not yet in the
equipartition state. For model Q4B3e15, the initial magnetic field is weak,
and hence,  the MRI is not well resolved and moreover the magnetic-field
amplification by winding is insufficient to reach the saturation in the early
stage after the merger. For this case, the magnetic-field energy keeps
increasing until $t\sim 60$ ms at which the disk settles eventually into an
equipartition state, and thus, the magnetic-field amplification saturates. The
main contributor to this amplification is clearly understood to be the
magnetic winding because $E_{B}$ increases in proportion to $t^{2}$
approximately for $t\approx 20$–$60$ ms. In this model, the MRI is not well
resolved until $t\sim 200$ ms because of the insufficient grid resolution
and/or the insufficient magnetic-field strength (see below), so that we can
only find the effect of winding in the magnetic-field amplification. However,
for a real system (or in an ideal computation with an infinite grid
resolution), the MRI should take place significantly from an earlier stage
after the merger. We may expect that in such realistic cases, the disk would
achieve the equipartition state earlier.  It is worthy to note that the onset
time of the post-merger mass ejection numerically found correlates with the
time at which the equipartition is established (cf. Fig. 6). Thus, the present
results illustrate that it is not easy to strictly identify the onset time of
the post-merger mass ejection by numerical simulation (see also a discusion in
Sec. III.2.3).

For model Q4B5tn, we find another remarkable behavior right after the merger:
During the merger stage $E_{\mathrm{B}}$ is amplified up to $\sim 4\times
10^{51}$ $\mathrm{erg}$, but it rapidly drops by an order of magnitude to
$\sim 4\times 10^{50}$ $\mathrm{erg}$. We do not see this drop for other
models. Our interpretation for this drop is that the magnetic-field
dissipation by reconnection near the equatorial plane occurs. For this model,
we initially embed a strong magnetic field with opposite polarities across the
equatorial plane. This magnetic-field configuration is the source of the
efficient magnetic-field reconnection. After this drop, $E_{\mathrm{B}}$
starts increasing again by winding, although we do not find the clear power
law proportional to $t^{2}$ because the magnetic-field energy is close to
saturation. At $t\approx 40$ ms $E_{\mathrm{B}}$ reaches saturation and starts
decreasing again. After the saturation, $E_{\mathrm{B}}/E_{\mathrm{int}}$
approaches asymptotically $\sim 10^{-2}$ as in other models.

The evolution process after the magnetic field saturates is qualitatively
identical irrespective of the initial magnetic-field strength, configuration,
neutron-star EOS, and equatorial-plane symmetry. Thus, we conclude that the
evolution process shown here is the universal one for black hole-neutron star
mergers that experience tidal disruption of neutron stars.

Figure 10 shows the evolution of an MRI quality factor, defined by

$\displaystyle Q_{z}:=\left\langle\left|\lambda_{z}/\Delta
x\right|\right\rangle_{\rm ave},$ (11)
$\displaystyle{\lambda_{z}}:=\frac{b_{z}}{\sqrt{4{\pi}{\rho}h+b^{\mu}b_{\mu}}}\frac{2{\pi}}{\Omega},$
(12)

where $\Omega$ is the local angular velocity, $z$-direction is the direction
of the rotation axis, and $\langle\cdots\rangle_{\rm ave}$ denotes the spatial
average with the weight of the rest-mass density for the region with $\rho\geq
10^{6}\,{\rm g/cm^{3}}$. For $Q_{z}>10$, we interpret that the MRI is
numerically well-resolved. It is found that for most of the models, $Q_{z}>10$
is achieved for $t\gtrsim 20$–$50$ ms, while for model Q4B3e15, $Q_{z}>10$ is
achieved only for $t\gtrsim 250$ ms.222A steep increase of $Q_{z}$ takes place
at $t\sim 200$ ms, and thus, the MRI activity is also partly visible already
for $t\sim 200$ ms in this model. This illustrates that for the model with
lower initial magnetic-field strengths, it takes a longer time until the
fastest growing mode of the MRI can be well resolved. As we already remarked,
this is an artifact due to the insufficient grid resolution in numerical
computation, and hence, in real systems, the MRI turbulence would be developed
from an earlier stage.

####  III.2.3 Effective viscosity

Figure 10: The time evolution of the MRI quality factor for models Q4B3e16
(left panel) and Q4B5n and Q4B5tn (right panel). The results for model Q4B5L
of our previous paper Hayashi  _et al._ (2022) are also shown for comparison
(in grey color).

Figure 11: The time evolution of the ratio of the magnetohydrodynamical
anisotropic stress (left) and Maxwell stress (right) to the pressure,
$\alpha_{ij}$ and $\alpha^{\mathrm{M}}_{ij}$, respectively. The results for
$r\varphi$, $xy$, $xz$, and $yz$ components are shown. The top two panels show
the results for model Q4B3e15, and the bottom two panels show the results for
models Q4B5n and Q4B5tn.

Associated with the development of the MRI turbulence and resulting dynamo
action, effective viscosity is enhanced in the accretion disk. We here analyze
an effective viscosity tensor by evaluating the ratio of the
magnetohydrodynamical anisotropic stress to the pressure, which is defined by

$\displaystyle\alpha_{ij}:=\left\langle\left|{1\over P}\left(\rho
h\hat{u}_{i}\hat{u}_{j}-{1\over
4\pi}b_{i}b_{j}\right)\right|\right\rangle_{\rm ave}.$ (13)

We also evaluate the ratio of the Maxwell stress to the pressure defined by

$\displaystyle\alpha^{\mathrm{M}}_{ij}:=\left\langle\left|{1\over
P}\left(-{1\over 4\pi}b_{i}b_{j}\right)\right|\right\rangle_{\rm ave}.$ (14)

Here, $i\not=j$ ($i$, $j=x,y,z$) and $\langle\cdots\rangle_{\rm ave}$ denotes
the spatial average with the weight of the rest-mass density for the region
with $\rho\geq 10^{7}\,{\rm g/cm^{3}}$. $\hat{u}_{i}$ is defined by
$u_{i}-\langle u_{i}\rangle_{t,{\rm ave}}$ where $\langle u_{i}\rangle_{t,{\rm
ave}}$ denotes the local time average of $u_{i}$. The time average needs to be
subtracted from $u_{i}$ to eliminate the contribution of coherent motion (not
random motion;  e.g., the orbital motion around the black hole) for evaluating
the anisotropic stress associated with the turbulent motion.

Figure 11 plots the time evolution of the off-diagonal components of
$\alpha_{ij}$ and $\alpha^{\mathrm{M}}_{ij}$. This shows that for $t\lesssim
100$ ms $\alpha_{ij}>O(0.1)$, but we interpret that this is not a physical
value, nor associated with the magnetohydrodynamics effect: The remnant matter
shows the non-axisymmetric structure for $\lesssim 100$ ms after the merger,
and in such a case, the value of $\alpha_{ij}$, specifically the contribution
from the Reynolds stress part, cannot be evaluated properly. Thus we focus
only on the stage for $t\gtrsim 100$ ms for which the non-axisymmetric
structure is not very appreciable and the MRI turbulence is developed.

When the disk is in a MRI turbulent stage, we find that the $r\varphi$ and
$xy$ components are $\approx 0.05$–$0.1$, and $xz$ and $yz$ components are
$\approx 0.02$–$0.05$. Hence, the order of the magnitude of $\alpha_{ij}$
agrees with the often-used value of the alpha viscous parameter for the
accretion disk Shakura and Sunyaev (1973), although the magnitude for each
component of $\alpha_{ij}$ has anisotropy. Our interpretation for this
anisotropy is that not only the MRI turbulence but also the effects by the
global magnetic fields such as magneto-centrifugal effects Blandford and Payne
(1982) contribute to the angular momentum transport because $\alpha_{ij}$ for
$r\varphi$ and $xy$ components are larger than the others. We also note that
the dominant part of $\alpha_{ij}$ stems from the Maxwell stress; the
contribution of the Reynolds stress, which originates from the fluid turbulent
motion, is $\sim 0.01$ irrespective of the model and component. This trend is
universally found for all the models.

Model Q4B3e15, which has a low initial magnetic-field strength, shows a factor
of $\sim 2$ smaller values of $\alpha_{ij}$ than for the other models, but
this is relatively minor compared to the difference in the initial magnetic
field strength (the initial magnetic field is smaller than the other models by
a factor of $\sim 17$). This result suggests that for this model, the fastest
growing mode of the MRI might be only partly resolved. Indeed, for this model,
a delay with $\sim 200$ ms in the magnetic-field amplification as well as in
all the processes of the disk evolution is found. This indicates that the disk
expansion due to the angular momentum transport is delayed due to the weaker
turbulence viscosity. As a result of this delay, the drop of the disk
temperature and neutrino luminosity delays, and thus, the post-merger mass
ejection is delayed by $\sim 200$ ms. Note, however, that the delay could be
the artifact due to the insufficient grid resolution.

Similar analyses were done for the simulations starting from a black hole-disk
system to evaluate the effective viscosity $\alpha_{r\varphi}$ in Ref. De and
Siegel (2021), and the value of $\alpha_{r\varphi}$ found in their analysis
were comparable to our present results. This is because the magnetic field
amplification is saturated as the disk reaches to equipartition state. We
note, however, that the result from the present simulation has a slightly
higher value. The existence of the additional magnetic field amplification
mechanism due to the non-axisymmetric structure of the disk and the minor
difference in the definition of $\alpha_{ij}$ are interpreted as the reasons.

####  III.2.4 Property of ejecta

Figure 12: Mass histogram as a function of the electron fraction of the ejecta
for models Q4B3e15 (top panel), Q4B5n and Q4B5tn (middle panel), and SFHoQ4B5
(bottom panel). The results for models Q4B5L and Q4B5H of our previous paper
are also shown for comparison (in grey color). Note that the vertical axis is
normalized by the total ejecta mass.

Figure 13: The same as Fig. 12 but for the mass histogram as a function of the
velocity of the ejecta.

Now, we turn our attention to the properties of the ejecta. Figures 12 and 13
show the rest-mass histogram as a function of the electron fraction
$Y_{\mathrm{e}}$ and velocity $v$ for the ejecta component, respectively. The
ejecta velocity $v$ is defined by

$\displaystyle v$ $\displaystyle:=$
$\displaystyle\sqrt{1-\Gamma_{\infty}^{-2}},$ $\displaystyle\Gamma_{\infty}$
$\displaystyle:=$ $\displaystyle-hu_{t}/h_{\mathrm{min}},$ (15)

where $\Gamma_{\infty}$ is interpreted as the terminal Lorentz factor under
the assumption that the internal energy is converted completely to the kinetic
energy of the ejecta in the far region.

There are two distinctive $Y_{\mathrm{e}}$ components for the ejecta as found
in our previous work Hayashi  _et al._ (2022) , where the comparison with the
results in the literature is also made. One is the dynamical ejecta for which
$Y_{\mathrm{e}}\approx 0.03$–$0.07$ irrespective of the simulation setups.
However, the range of $Y_{\mathrm{e}}$ for the dynamical ejecta depends
slightly on the EOS: For model SFHoQ4B5, $Y_{\mathrm{e}}\approx 0.03$–$0.06$,
while for model Q4B5H with DD2 EOS, $Y_{\mathrm{e}}\approx 0.03$–$0.07$
Kyutoku  _et al._ (2018). The electron fraction of the dynamical ejecta
directly reflects the neutron richness of the neutron star because the
dynamical ejecta is affected only weakly by the thermal and weak-interaction
process in the merger and post-merger stages. Thus, the difference in the
distribution of $Y_{\mathrm{e}}$ for the dynamical ejecta between models
SFHoQ4B5 and Q4B5H comes directly from the difference in the neutron star
EOSs.  This difference is reflected in the product of the heavy elements in
the $r$-process nucleosynthesis Wanajo  _et al._ (2022).

For the post-merger ejecta, the electron fraction is as higher as $0.1\lesssim
Y_{\mathrm{e}}\lesssim 0.4$. Irrespective of the initial magnetic-field setups
and equatorial-plane symmetry, similar distribution is found for the DD2
models, while the difference in the EOS makes a quantitative difference.
Comparing the high-resolution models SFHoQ4B5 and Q4B5H, we find that model
SFHoQ4B5 has a distribution with higher $Y_{\mathrm{e}}$ values: e.g, the
maximum value of $Y_{\mathrm{e}}$ is smaller than 0.4 for all the DD2 models
while it is $\approx 0.42$ for the SFHo model. Our speculation for the reason
to this is as follows: The neutron star modeled by the SFHo EOS is tidally
disrupted at an orbit closer to the black hole because it has a smaller
neutron-star radius. Then, the matter that forms the one-armed structure right
after the merger and subsequently forms an accretion disk experiences stronger
compression and shock heating between the inner and the outer spiral arms. As
a result, the disk temperature is enhanced right after the tidal disruption
and the value of $Y_{\mathrm{e}}$ is also increased (see also Ref. Kyutoku
_et al._ (2018)).

There are also two components in the mass histogram as a function of the
ejecta velocity (see Fig. 13). The high-velocity component with $v/c\gtrsim
0.1$ stems primarily from the dynamical ejecta, while the low-velocity
component stems primarily from the post-merger ejecta. For the DD2 models, the
velocity histogram has approximately identical distributions irrespective of
the setups, and it also agrees with previous results Hayashi  _et al._ (2022).
For the SFHo model, the low-velocity component has a higher fraction than that
for the DD2 models, because the dynamical ejecta mass is smaller while the
post-merger ejecta mass is comparable with that of the DD2 models.

###  III.3 Magnetic field in the funnel region and the relation to the short-
hard gamma-ray burst

####  III.3.1 Poynting luminosity and magnetic field in the funnel

Figure 14: The time evolution of isotropic-equivalent Poynting luminosity
$L_{\mathrm{iso}}$ for models Q4B3e15 (top panel), Q4B5n and Q4B5tn (middle
panel), and SFHoQ4B5 (bottom panel). The results for models Q4B5L and Q4B5H of
our previous paper Hayashi  _et al._ (2022) is also shown for comparison (in
grey color). The characters “u” and “l” denote the upper and lower
hemispheres, respectively.

Figure 15: The profile of the average toroidal magnetic field along the polar
direction ($\theta$) at $r\approx 150~{}{\rm km}$ as a function of time for
models Q4B3e15 (top panel), Q4B5n (middle), and Q4B5tn (bottom).  Figure 16:
The snapshots of the toroidal magnetic field (color profile) together with the
poloidal magnetic-field lines (curves) on the $x$-$z$ plane at selected time
slices for model Q4B5n. See also the following link for an animation:
https://www2.yukawa.kyoto-u.ac.jp/~kota.hayashi/Q4B5n-mf.mp4.

As in our previous simulations Hayashi  _et al._ (2022), in addition to the
dynamical and post-merger ejecta, we find a launch of the collimated outflow
of matter and Poynting flux in the magnetosphere, which is established after
the magnetic tower effect associated with the black-hole spin develops a
large-scale helical magnetic fields around the spin axis of the black hole.

Figure 14 shows the time evolution of $L_{\rm iso}$: an isotropic-equivalent
Poynting luminosity, which we define using the Poynting luminosity for
$\theta<10^{\circ}$ and $r\approx 1500\,{\rm km}$ as

$\displaystyle L_{\rm
iso}:=\frac{2}{1-\cos(10^{\circ})}L_{\theta<10^{\circ},r\approx 1500\,{\rm
km}},$ (16)

where

$\displaystyle L_{\theta<10^{\circ},r\approx 1500\,{\rm km}}$
$\displaystyle:=-\int_{\theta<10^{\circ},r\approx 1500\,{\rm
km}}{T^{\mathrm{(EM)}}}^{~{}r}_{t}\sqrt{-g}dS_{r}.$ (17)

$T^{\mathrm{(EM)}}_{~{}\mu\nu}$ denotes the electromagnetic part of the
energy-momentum tensor,  i.e., the second line of Eq. (10).  The polar angle
is defined by $\theta=\tan^{-1}(\sqrt{x^{2}+y^{2}}/z)$. A particular value
($10^{\circ}$) is chosen for the surface integral because the opening angle of
the funnel region is initially as narrow as $\sim 10^{\circ}$. We always
assume an observer located along the $z$-axis for evaluating
$L_{\mathrm{iso}}$.

For models Q4B5n and Q4B5tn, the typical maximum value of $L_{\rm iso}$ is
$O(10^{50})\,{\rm erg/s}$. The high-Poynting luminosity stage, which is
designated by $L_{\mathrm{iso}}\gtrsim 3\times 10^{49}$ erg/s in this paper,
is identified for $t\approx 300$–2500 ms for model Q4B5n. For model Q4B5tn,
the high-Poynting luminosity stage for the upper and lower hemispheres is
identified for $t\sim 500$–1500 ms and $\sim 1000$–1600 ms, respectively.
During the high-Poynting luminosity stage, $L_{\mathrm{iso}}$ varies with time
by more than an order of magnitude for these models, reflecting the variation
of the magnetic-field strength and configuration. The reason for this
variability is that the magnetic fields with high field strengths are often
provided from the disk to the black hole  and the polar region by the MRI
dynamo activity. Figure 15 shows the average value of the toroidal field
$b_{\bar{\varphi}\mathrm{,ave}}$ as a function of time and polar angle
$\theta$ for models Q4B3e15, Q4B5n, and Q4B5tn. Here, $x$, $y$, and $z$ are
defined with respect to the black-hole center (the location of the puncture).
The toroidal field is defined by
$b_{\bar{\varphi}}=(xb_{y}-yb_{x})/\sqrt{x^{2}+y^{2}}$. The average is taken
with respect to the azimuthal angle $\varphi=\tan^{-1}(y/x)$ at the selected
radius of $r:=\sqrt{x^{2}+y^{2}+z^{2}}\approx 150$ km. Figure 15 displays the
so-called butterfly structure Brandenburg and Subramanian (2005); the polarity
of the toroidal magnetic field flips in a quasi-periodic manner with the
period of $\sim 20$ local orbital periods. Also, the strong magnetic-field
fluxes continuously ascend from the equatorial plane to the surface of the
accretion disk.

Although this butterfly structure induced by the MRI dynamo was already found
in our previous models Hayashi  _et al._ (2022), we find an interesting new
feature in our new computation; the MRI dynamo activity in the accretion disk
determines the magnetic-field structure in the magnetosphere. In the dynamo
activity, the magnetic fields often ascend from the disk to the vertical and
polar regions. For most of the cases, they do not cancel out or alternate the
fields originally stayed in the polar region, and thus, the polarity of the
magnetosphere is unchanged. However, for exceptional cases, the inversion of
the polarity is achieved. For model Q4B5n, this occurs at $t\sim 1.1$ s and
$1.4$ s (see the middle panel of Fig. 15). For $t\lesssim 1.1$ s, the polarity
of $b_{\bar{\varphi}\mathrm{,ave}}$ is positive in the polar region of the
upper hemisphere. Then at $t\sim 1.1$ s, the polarity flips to negative,
following the polarity flip at an inner region of the accretion disk.
Subsequently, at $t\sim 1.4$ s, the polarity flips back to positive. For
$L_{\mathrm{iso}}$ of this model, there are three characteristic peaks at
$t\sim 0.3$ s, $1.3$ s, and $1.5$ s. These peaks reflect the variation of the
butterfly structure at the polar region: During the polarity flips in
progress, the intensity of the outgoing Poynting flux and $L_{\mathrm{iso}}$
naturally drop, because the magnetic field in the polar region is not aligned
and the magnetosphere loses a coherency with respect to the magnetic-fields
lines.

This polarity flip in the magnetosphere is also found in the snapshots of the
magnetic-field structure. Figure 16 shows the toroidal and poloidal magnetic
field structures on the $x$-$z$ plane at $t\sim 0.3$ s, $1.3$ s, $1.5$ s, and
$2.7$ s for model Q4B5n. The first three panels correspond to the snapshots at
which $L_{\mathrm{iso}}$ is at local peaks. The figure clearly shows the
polarity flip of both poloidal and toroidal magnetic fields in the
magnetosphere. Around these time ranges, the magnetic fields ascending from
the disk reconnect the originally-existing fields in the magnetosphere, and
subsequently, the polarity is changed.

The middle and bottom panels of Fig. 15 also show the magnetic-field polarity
flip near the polar region for $t>2.5$ s of model Q4B5n and for $t\sim 1.0$ s
on the upper hemisphere of model Q4B5tn. However, for these stages, no peak in
$L_{\mathrm{iso}}$ is found and its typical value is lower than $\sim 10^{49}$
erg/s, which is an order of magnitude lower than the peak luminosity. Our
interpretation for this is that the magnetic fields ascending from the disk
due to the MRI dynamo activity disturb or deform the magnetosphere, but are
not strong enough or aligned enough to replace the polarity of the field
completely and reform the magnetosphere that can launch a high-intensity
Poyting flux with $L_{\mathrm{iso}}\sim 10^{50}$ erg/s.

We note that the polarity flip was already reported in magnetohydrodynamics
simulations for the accretion disks around a spinning black hole Christie  _et
al._ (2019); Liska  _et al._ (2020). For these simulations, the authors also
found a turbulent state of the accretion disks. Thus, the polarity flip is
likely to occur often, if magnetic fluxes with high field strengths are
ejected from the inner region of the disks in a turbulent state.

For both models Q4B3e15 and SFHoQ4B5, the maximum value of $L_{\mathrm{iso}}$
is $\sim 5\times 10^{49}$ erg/s, which is slightly lower than those of other
models for which typically $L_{\mathrm{iso}}\sim 10^{50}$ erg/s in the bright
stages. Our interpretation is that this is due to the lower magnetic-field
strength in the magnetosphere. For these models, the rest-mass density in the
disk at the time when the magnetosphere is formed is lower. The reason for
model Q4B3e15 is that it takes a longer time to form the magnetosphere than
for the other DD2 models and the reason for model SFHoQ4B5 is that the disk
mass is smaller than for the DD2 models. As we already mentioned, the field
strength in the magnetosphere is determined by the field strength of the disk
at which the equipartition state is achieved. For model Q4B3e15, the
equipartition is achieved in the relatively late stage, at which the rest-mass
density and internal energy of the disk are relatively low. This leads to a
lower magnetic-field strength in the magnetosphere for these two models. As a
result of the lower field strength, the Poynting luminosity, which is powered
by the Blandford-Znajek mechanism, becomes lower. The lower maximum value of
$L_{\mathrm{iso}}$ for model SFHoQ4B5 is understood as the physical result,
while that for model Q4B3e15 could be the numerical artifact due to the
insufficient grid resolution.

The high-Poynting luminosity stage for model Q4B3e15 starts at $\sim 2600$ ms
and lasts for $\sim 400$ ms, entering the fading stage at $\sim 3000$ ms. For
model SFHoQ4B5, the high-Poynting luminosity stage starts at $\sim 1100$ ms.
We do not find a clear fading stage for this model, but $L_{\mathrm{iso}}$
appears to gradually decrease to $\sim 10^{49}$ erg/s at the termination of
the simulation. We indeed find for this model that the opening angle increases
with a timescale of a few seconds, and thus, we expect that $L_{\mathrm{iso}}$
will eventually drop in this timescale.

The isotropic-equivalent Poynting luminosity of $L_{\mathrm{iso}}\sim 10^{50}$
erg/s together with the opening angle of $\theta\sim 10^{\circ}$ (cf. Fig. 2)
fairly agrees with those for short-hard gamma-ray bursts (or at least for low-
luminosity short-hard gamma-ray bursts) in the assumption that the conversion
efficiency of the Poynting flux to the gamma-ray radiation is sufficiently
high (i.e., close to unity) Nakar (2007); Berger (2014).

The jet efficiency, which is defined by the ratio of the Poynting luminosity
to mass accretion rate on the apparent horizon is $\lesssim 10^{-2}$ in the
present simulation. Although the value is increasing, it is still much smaller
than that for Ref. Christie  _et al._ (2019), which increases up to $\sim 1$.
This difference is mostly because the infalling mass flux near the equatorial
plane is still high in our simulation. Around the polar region, however, the
ratio of the Poynting flux to the infalling mass flux is $\sim 1$ locally, and
through such a region, the black hole is capable of powering outgoing Poynting
luminosity into the magnetosphere (see also the discussion of the next
subsection).

####  III.3.2 MADness parameter

Figure 17: The time evolution of the MADness parameter for models Q4B3e15 (top
panel), Q4B5n and Q4B5tn (middle panel), and SFHoQ4B5 (bottom panel). The
results for models Q4B5L and Q4B5H of our previous paper Hayashi  _et al._
(2022) are also shown for comparison (in grey color). The characters “u” and
“l” denote the upper and lower hemispheres, respectively.

Figure 18: The same as Fig. 17 but for the time evolution of the magnetic flux
evaluated on the apparent horizon.

Figure 19: The profile of the local MADness parameter
$\phi_{\mathrm{AH,local}}$ along the polar direction ($\theta$) as a function
of time for models Q4B3e15 (top panel), Q4B5n (middle), and Q4B5tn (bottom).

Figure 20: The angular distribution of the Poynting flux per steradian on a
sphere of $r\approx 1500$ km for model Q4B5tn at selected time slices. The
left and right panels display the upper and lower hemispheres, respectively.
The bright color displayed in the polar region stems from the Blandford-Znajek
effect, while for other regions, the magnetic field accompanied by the
outflowing matter contributes mainly to the Poynting flux. The region for
which the Poynting flux is intense moves in the direction of the $y$-axis from
the vicinity of the pole. At the same time, the opening angle of the Poynting
flux gradually increases. See also the following link for the animation:
https://www2.yukawa.kyoto-u.ac.jp/~kota.hayashi/Q4B5tn-pf.mp4.

Figure 17 shows the time evolution of the so-called MADness parameter which is
defined by Tchekhovskoy  _et al._ (2011)

$\displaystyle\phi_{\mathrm{AH}}:=\frac{\Phi_{\mathrm{AH}}}{\sqrt{\dot{M}_{\mathrm{AH}}4\pi
r_{\mathrm{AH}}^{2}c}},$ (18)

where

$\displaystyle\dot{M}_{\mathrm{AH}}$ $\displaystyle:=$
$\displaystyle\oint_{\mathrm{AH}}\rho\sqrt{-g}u^{i}dS_{i}$ (19)
$\displaystyle\approx$
$\displaystyle\oint_{r=\max{(r_{\mathrm{AH}}})}\rho_{*}v^{r}r^{2}\sin{\theta}d\theta
d\varphi,$

and

$\displaystyle\Phi_{\mathrm{AH}}$ $\displaystyle:=$
$\displaystyle\oint_{\mathrm{AH}}B^{i}\sqrt{\gamma}dS_{i}$ (20)
$\displaystyle\approx$
$\displaystyle\oint_{r=\max{(r_{\mathrm{AH}}})}\mathcal{B}^{r}r^{2}\sin{\theta}d\theta
d\varphi,$

which denotes the magnetic flux penetrating the apparent horizon. It is found
that for all the models $\phi_{\mathrm{AH}}$ remains smaller than 10 within
the simulation time. Thus, the accretion disks in our simulations do not
satisfy the often-referred condition for the magnetically-arrested disk (MAD),
$\phi_{\mathrm{AH}}\gtrsim 50$ Tchekhovskoy  _et al._ (2011).333 The present
results are quantitatively similar to the model BT in Ref. Christie  _et al._
(2019) as we discussed in our previous work Hayashi  _et al._ (2022). However,
as we have already described, a high-intensity Poynting flux is generated even
if the condition for the MAD is not satisfied. Thus, the MADness parameter
$\phi_{\mathrm{AH}}$ may not be suitable for assessing whether the jet is
launched in the context of neutron-star mergers.

Figure 18 shows the time evolution of $\Phi_{\mathrm{AH}}$. For models Q4B5n
and Q4B5tn, the magnetic fluxes penetrating the upper and lower hemisphere of
the apparent horizon are shown separately, although these two components
approximately agree with each other. For most of the models, the magnetic flux
on the apparent horizon reaches its peak in the timescale similar to that for
achieving the equipartition in the accretion disk. Since $\Phi_{\mathrm{AH}}$
does not increase after the peak is reached, the decrease of the accretion
rate is the only path for the increase of the MADness parameter for such
models. For model Q4B5tn, which initially has the toroidal magnetic field in
the neutron star, $\Phi_{\mathrm{AH}}$ is low for the early post-merger stage
($t\lesssim 50$ ms) and it increases gradually in the entire simulation time.
For this case, field lines starting from a point in each hemisphere always end
in the same hemisphere for the early stage. However, due to the MRI turbulence
and dynamo, the poloidal field is developed and subsequently penetrates the
black-hole horizon, resulting in the increase of $\Phi_{\mathrm{AH}}$. By
contrast, for the pure poloidal initial field, the black hole horizon is
penetrated by the poloidal field from the early stage.

Because the MADness parameter defined on the entire horizon surface might not
be a good indicator for assessing the launch of the strong Poynting flux,
instead of it, we propose another parameter based on the local quantities. The
point is that the magnetic-field lines that generate the strong Poynting flux
do not penetrate the black-hole horizon in the vicinity of the equator, at
which dense matter infalling from the accretion disk is always present and the
(low-beta) magnetosphere is not formed. This suggests that focusing on the
polar region on the apparent horizon for evaluating the MADness-like quantity
would be a better strategy. Thus, we introduce a local MADness parameter,
$\phi_{\mathrm{AH,local}}$, which is defined by

$\displaystyle\phi_{\mathrm{AH,local}}:=\frac{\mathcal{B}^{r}}{\sqrt{\rho_{*}v^{r}c}}.$
(21)

Figure 19 shows the azimuthal-average value of the local MADness parameter
$\phi_{\mathrm{AH,local}}$ as a function of time and polar angle $\theta$ for
models Q4B3e15, Q4B5n, and Q4B5tn. For all the models, we find time intervals
with $L_{\mathrm{iso}}\gtrsim 3\times 10^{49}$ erg/s, and for such time
intervals we always find $\phi_{\mathrm{AH,local}}\gtrsim 100$ at polar region
of the apparent horizon. This suggests that the black hole has the ability to
form a magnetosphere and launch a jet if the value of
$\phi_{\mathrm{AH,local}}$ at the polar region exceeds $100$, even if the
value of $\phi_{\mathrm{AH}}$ is smaller than 50.

We note here that $\phi_{\mathrm{AH,local}}$ only gives us the necessary
condition for the launch of a jet with high-Poynting luminosity. The
disturbance or the deformation of the magnetosphere far from the horizon
associated with the evolution of the accretion disk could result in a low
value of $L_{\mathrm{iso}}$. For example the local MADness parameter exceeds
100 at $t\gtrsim 1700$ ms but $L_{\mathrm{iso}}$ falls below $10^{49}$ erg/s
for Q4B5tn model. This is due to the deformation (tilt) of the magnetosphere
induced by the post-merger mass ejection. The details of this behavior are
given in the next subsection.

Focusing on the polarity, the local MADness parameter for model Q4B5n shows
interesting behavior. As we already pointed out in this subsection, this model
shows a clear butterfly structure of $b_{\bar{\varphi}\mathrm{,ave}}$
extending to the polar region due to the flip of the magnetic-field polarity,
and this flip occurs at the peaks of $L_{\mathrm{iso}}$. A similar flip is
also observed for the local MADness parameter for $t=1000$–$1500$ ms (see Fig.
19). Thus, we conclude that the complete flip of the magnetic-field polarity
in the magnetosphere is the result of this polarity flip on the apparent
horizon. Just like in the MRI dynamo and butterfly structure of
$b_{\bar{\varphi}\mathrm{,ave}}$ in the accretion disk, the polarity flip
starts in the vicinity of the equatorial plane and propagates to the polar
region. This is because the value of $\phi_{\mathrm{AH,local}}$ is lower than
$10$ near the equatorial plane and the fluid dynamics dominates over the
magnetic-field dynamics. This feature enables the matter accretion from the
disk to occasionally carry the magnetic field with opposite polarity. Then,
the magnetic tower effect enhances the magnetic-field strength along the polar
direction. In this process, the preexisting magnetic field near the pole is
dissipated away due to the reconnection by the magnetic field with opposite
polarity ascending from the equatorial region, which replaces the polarity of
the field penetrating the polar region of the horizon. Once the magnetic tower
effect is in action, the magnetic field is amplified by the winding
(associated with the black-hole spin) and the matter is pushed outward. As a
result, a high-$\phi_{\mathrm{AH,local}}$ region is realized near the pole,
where the magnetic-field strength is high and the rest-mass accretion rate is
low. In the high-$\phi_{\mathrm{AH,local}}$ region, the magnetic-field
dynamics dominates the fluid dynamics, and hence, the polarity flip cannot
start in the polar region.

####  III.3.3 Time duration for high Poynting luminosity

The stage with a high value of $L_{\mathrm{iso}}$ continues for $\sim
400$–$2200$ ms and subsequently starts decreasing. This is particularly clear
for models Q4B3e15, Q4B5n, and Q4B5tn. For these models, we confirmed that the
value of $L_{\mathrm{iso}}$ decreases by nearly two orders of magnitude in the
fading stage. This is consistent with the duration of the short-hard gamma-ray
bursts, whose typical duration is $\sim 1$ s Nakar (2007).

In our previous paper Hayashi  _et al._ (2022), we discussed that the decrease
of $L_{\mathrm{iso}}$ is due to the increase in the opening angle of the
funnel region and the resulting decrease of the magnetic-flux density in the
magnetosphere. The opening angle of the strong Poynting-flux region increases
from $\lesssim 10^{\circ}$ to $\sim 20^{\circ}$. This is directly related to
the accretion disk evolution. The location of the funnel wall is determined by
the balance between the gas pressure of the thick torus and the magnetic
pressure at the funnel wall. In the seconds-long evolution of the torus
(disk), the rest-mass density and the gas pressure at the funnel wall
gradually decrease due to the post-merger mass ejection and matter accretion
onto the black hole. On the other hand, the magnetic pressure in the
magnetosphere and at the funnel wall does not decrease significantly, in
particular for the late stage of the evolution. Thus the magnetic pressure can
eventually exceed the gas pressure at the original position of the funnel
wall, resulting in the gradual expansion of the funnel region. As discussed
previously, this mechanism could be one of the ingredients that determine the
time duration of short-hard gamma-ray bursts.

For the SFHoQ4B5 model, the opening angle also increases with time but with a
longer timescale. As a result, we do not find the clear fading stage of
$L_{\mathrm{iso}}$ for this model in our simulation time. Due to the high
computational cost, we terminated the simulation at $t\sim 2500$ ms, but if we
evolve the system longer, the fading stage is likely to be present.

From the results for models Q4B3e15, Q4B5n, and Q4B5tn, we find two
additionally possible mechanisms for the fade-away of the Poynting luminosity,
which could also explain the time duration of short-hard gamma-ray bursts. For
both mechanisms, the non-trivial evolution of the magnetosphere associated
with the evolution of the accretion disk is essential.

For models Q4B3e15 and Q4B5n, in the very late stage of our simulation, the
aligned magnetic field is dissipated away, and as a result, the magnetosphere
with aligned magnetic fields disappears. Figure 16 shows that for $t\lesssim
1550$ ms the poloidal magnetic-field lines in the polar region are aligned
approximately with the black-hole spin axis ($z$-axis). In the magnetosphere
along the $z$-axis, the magnetic-field lines are clearly helical and maintain
the high-intensity Poynting flux. However, it is found that at $t\approx 2730$
ms the magnetic field in the polar region is not aligned anymore. Moreover,
the magnetic-field dynamics cannot govern the fluid dynamics and the clear
magnetosphere disappears. In the absence of the well-ordered magnetic field,
the system cannot maintain the high-intensity Poynting flux.

Figure 15 shows that for model Q4B5n with $t\sim 1500$–$2500$ ms for which the
high Poynting luminosity is maintained (see Fig. 14), the polarity of the
poloidal magnetic field at the polar region remains to be preserved and does
not reverse. However, for $t\gtrsim 2500$ ms, the Poynting luminosity
decreases. For this late stage, the polarity in the polar region frequently
reverses in response to the polarity reversal in the disk near the equatorial
plane. The local MADness parameter in Fig. 19 also shows similar behavior. We
interpret that the disappearance of the magnetosphere and the decrease of the
Poynting flux stem from the MRI dynamo activity in the accretion disk. Due to
the MRI dynamo activity, the polarity of the magnetic field in the accretion
disk is reversed quasi-periodically and the magnetic flux continuously ascends
toward the polar region from the equatorial region. If the magnetic field
ascending from the disk has the polarity different from that in the
magnetosphere and its field strength is high enough, the magnetic field that
is originally located in the magnetosphere could be dissipated away by the
magnetic-field reconnection. This is what is observed for $t\gtrsim 2500$ ms
of model Q4B5n, and for $t\gtrsim 3000$ ms of Q4B3e15 model. If the strength
of the ascending magnetic field is even stronger, then it could replace the
polarity completely and reform the magnetosphere again. This is what is
observed for $t\approx 1000$–$1500$ ms of model Q4B5n. However, we expect that
the magnetosphere formation will not occur again after the disappearance of
the strong magnetic fields for models Q4B3e15 and Q4B5n because as already
discussed, the magnetic-field strength in the disk is determined by the rest-
mass density (and thus internal energy density) achieved for the equipartition
relation. Since the rest-mass density decreases as a result of the disk
expansion, the post-merger mass ejection, and the mass accretion onto the
black hole, the magnetic-field strength in the disk also decreases. Therefore,
the revival of a magnetosphere with strong magnetic fields is not possible in
the late stage of the disk. In the absence of strong magnetic fields, the
Poynting luminosity is low because the Blandford-Znajek Poynting luminosity is
proportional to the square of the magnetic-field strength Blandford and Znajek
(1977).

Model Q4B5tn shows the second mechanism. In this mechanism, the magnetosphere
still exists, but the Poynting luminosity becomes apparently low if we measure
along the direction of the black-hole spin. The mechanism of this is the tilt
of the accretion disk and magnetosphere. Figure 20 shows the angular
distribution of the Poynting flux per steradian defined by
$-T^{~{}r}_{t}\sqrt{-g}/\sin{\theta}$ on a sphere of $r\approx 1500$ km for
model Q4B5tn. For $t\lesssim 1300$ ms the strong Poynting-flux region is
approximately aligned with the polar direction ($z$-axis direction and black-
hole spin direction). However, at $t\sim 2000$ ms the strong Poynting-flux
region starts deviating from the polar region clearly, and at $t\approx 2900$
ms it is tilted by $25^{\circ}$–$40^{\circ}$ from the $z$-axis direction in
the $y$-$z$ plane. The isotropic-equivalent luminosity is defined by the
integration for $0^{\circ}<\theta<10^{\circ}$ and
$170^{\circ}<\theta<180^{\circ}$ on the upper and lower hemispheres,
respectively, and thus, it decreases significantly by the tilt. We note that
the black-hole spin axis is aligned with the $z$-axis during the entire
simulation time. As Fig. 3 shows, the system approximately has the equatorial-
plane symmetry in the early stage of the post-merger evolution. However,
later, the system loses the symmetry; the accretion disk misaligns with the
equatorial plane and the major axis of the magnetosphere also deviates from
the $z$-axis direction. This is due to the asymmetric nature of the post-
merger mass ejection resulting from the turbulent state of the accretion
disk.444 This asymmetric mass ejection is always present in the absence of
equatorial symmetry, and the origin is a small turbulence in the disk. The
asymmetry is enhanced subsequently. It is difficult to mention what happens in
the long-term evolution process, but one important point is that there is no
mechanism that suppresses the enhancement of the asymmetry. Thus, once the
equatorial asymmetry is induced in the matter distribution, the asymmetry is
always enhanced, and as a result, the asymmetric mass ejection takes place.
Specifically, in this model, the post-merger ejecta carries a large amount of
the $y$-component of the angular momentum and this breaks the symmetry of the
accretion disk. Indeed, the remnant accretion disk has the $y$-component of
the angular momentum, and as a result, the system tilts in the $y$-$z$ plane.
555We note that the system is still capable of launching a jet in an off-axis
direction. Indeed we still find the strong Poynting-flux region even after
$L_{\mathrm{iso}}$ (in our definition) drops.

We find not only the tilt of the high Poynting-flux region, but also the
widening of the opening angle in Fig. 20 like in the models of our previous
simulations Hayashi  _et al._ (2022). For model Q4B5tn in Fig. 20, at
$t\approx 850$ ms, the opening angle is $\lesssim 10^{\circ}$ but it increases
to $\sim 30^{\circ}$ at $t\approx 2900$ ms. Thus, the intensity of the
Poynting flux decreases simultaneously with tilting. This widening of the
opening angle results from the decrease of the gas pressure at the funnel wall
due to the post-merger mass ejection and matter accretion onto the black hole,
as we already described.

We have found the two new possible mechanisms for the decrease of the Poynting
luminosity in addition to one mechanism that we already found in our previous
paper Hayashi  _et al._ (2022). For all the three mechanisms, the evolution of
the accretion disk and the post-merger mass ejection are critical processes.
Our results show that irrespective of the mechanisms, the timescale of $\sim
1$ s for the high Poynting-luminosity stage is determined by the evolution
timescale of the accretion disk, which is determined by the neutrino cooling
and magnetohydrodynamics turbulence that control the post-merger mass
ejection.

A word of caution is appropriate here. The system, specifically the accretion
disk, is in a turbulent state by the MRI. That is, the evolution of the disk
and MRI dynamo activity are determined by a stochastic process. This implies
that we cannot precisely predict the strength of the magnetic field that
penetrates the black hole and forms the magnetosphere and the angular momentum
that is carried away by the post-merger ejecta. It is also not easy to predict
whether the black hole will always be penetrated by the magnetic field strong
enough to form a magnetosphere and launch the Poynting flux that can explain
typical short-hard gamma-ray bursts. It is also not easy to predict by which
mechanisms the high Poynting-luminosity stage is terminated. However, our
simulation results show that once the magnetosphere is formed, its subsequent
evolution is determined by the global properties of the post-merger ejecta and
the accretion disk.

##  IV Conclusion

We have reported the new results of general-relativistic neutrino-radiation
magnetohydrodynamics simulation for seconds-long black hole-neutron star
mergers. We employed a variety of setups that are different from those in our
previous paper Hayashi  _et al._ (2022), while the mass of the black hole and
neutron star and the black-hole spin are identical with those of the previous
setting so that the neutron star is disrupted by the tidal force of the black
hole. The difference in the setup for the present work from the previous work
is summarized as follows: (1) The initial maximum magnetic-field strength in
the neutron star is assumed to be $3\times 10^{15}$ G, which is by a factor 17
weaker than in the previous simulations. (2) The toroidal magnetic field is
assumed for the initial magnetic-field configuration in the neutron star. (3)
The equatorial-plane symmetry is removed from the simulation. (4) The SFHo EOS
in addition to the DD2 EOS is employed to model the neutron star. We performed
four new simulations for $\sim 2.5$–$6$ s in order to self-consistently follow
the dynamical mass ejection, remnant disk evolution, post-merger mass
ejection, and generation of collimated Poynting flux in the magnetosphere
which may drive a short-hard gamma-ray burst.

We found, irrespective of the difference in the setups listed above, that the
essential part of the merger and post-merger processes is unchanged. First,
the matter with its mass $\approx 0.046M_{\odot}$ (DD2 models) or $\approx
0.014M_{\odot}$ (SFHo model) is dynamically ejected right after the tidal
disruption of the neutron star in the timescale of a few ms. Then the
accretion disk with the initial rest mass $M_{\mathrm{disk},0}\approx
0.28M_{\odot}$ (DD2 models) or $\approx 0.22M_{\odot}$ (SFHo model) is formed
around the remnant black hole. In the accretion disk, the magnetic field is
amplified by magnetohydrodynamics effects such as the MRI, winding, and
Kelvin-Helmholtz instability. The MRI turbulence and dynamo action induce
effective viscosity that enhances the angular-momentum transport. As a result,
the mass accretion and the expansion of the disk are induced. In addition, a
portion of the matter with the strong magnetic fields is outflowed vertically
from the disk. This plays a role in the disk expansion to the vertical
direction. In the turbulent process, thermal energy is generated, but in the
first several hundred ms, the generated heat is dissipated primarily by the
neutrino emission and it does not have a significant effect on the post-merger
mass ejection.

The disk expands gradually due to the angular momentum transport effects. As a
result, the maximum temperature of the disk drops below $\sim 3$ MeV, and the
neutrino luminosity drops below $\sim 10^{51}$–$10^{51.5}$ erg/s. Then the
neutrino emission cannot carry away an appreciable fraction of the thermal
energy generated by the turbulent process from the disk, and the generated
thermal energy induces the convective motion in the disk, which carries the
thermal energy generated around the inner edge of the disk to the outer
region. This convective motion contributes to the heating for the outer part
of the disk, and eventually, induces the post-merger mass ejection. In our
present setups, the mass of the post-merger ejecta is $\gtrsim 0.030M_{\odot}$
(DD2 model), or $\gtrsim 0.018M_{\odot}$ (SFHo model), which are $\sim 10\%$
of $M_{\mathrm{disk},0}$.

As in our previous paper Hayashi  _et al._ (2022), we found, irrespective of
the models, that there are two components in the electron fraction
distribution for the ejected matter. One is a low-electron fraction component
($Y_{\mathrm{e}}<0.1$) produced by the dynamical ejecta and the other is a
mildly neutron-rich component ($0.1\lesssim Y_{\mathrm{e}}\lesssim 0.4$)
produced by the post-merger ejecta.  We found that the electron fraction
distribution depends weakly on the EOS, and the dependence is reflected in the
elemental pattern of the $r$-process nucleosynthesis Wanajo  _et al._ (2022).
Also, there are two components in the velocity distribution. One is a fast
component (up to $v\sim 0.4c$) produced by the dynamical ejecta and the other
is a relatively slow component ($v<0.1c$) produced by the post-merger ejecta.
These distributions are suitable for reproducing an elemental abundance
pattern similar to the solar abundance and those of the metal-poor stars
Wanajo  _et al._ (2022).

We analyzed the ratio of the anisotropic stress to the pressure to evaluate
the alpha viscous tensor, $\alpha_{ij}$. It is found that all the components
of $\alpha_{ij}$ have a value between $10^{-2}$ and $10^{-1}$. However, the
value of the $r\varphi$ or $xy$ components of this tensor is larger than the
other components. This reflects that not only the magnetically-induced viscous
effect but also intrinsic magnetohydrodynamics effects play a role in the
momentum transport. The large value of $\alpha_{r\varphi}$ suggests that a
magnetohydrodynamics effect associated with global magnetic fields such as the
magneto-centrifugal effects Blandford and Payne (1982) play an important role
in the angular momentum transport.

We found quantitative differences between the results for the models with low
and high initial magnetic-field strengths. Because it takes a longer time to
amplify the magnetic field up to saturation and to achieve an equipartition
state in the disk from the low initial field strength (with a limited grid
resolution), the evolution of the accretion disk and post-merger mass ejection
are delayed. However, essentially no differences are found in the properties
of the ejected matter.  We note that even with the low initial magnetic-field
strength employed here, the initial field strength is still several orders
higher than the realistic value. The simulation starting from a neutron star
endowed with $\sim 10^{12}$ G is ultimately necessary.

Accompanying the turbulent disk formation, a funnel-shaped magnetosphere with
the low rest-mass density and the aligned helical magnetic-field lines are
formed near the spin axis of the black hole. This magnetosphere is a
magnetically dominated region and is in an approximate force-free state. The
helical magnetic field lines that form the magnetosphere penetrate the black
hole and it extracts the rotational kinetic energy of the rapidly spinning
black hole by the Blandford-Znajeck mechanism Blandford and Znajek (1977).
Then, the collimated outgoing Poynting flux is generated with the opening
angle of $\sim 10^{\circ}$, and its isotropic-equivalent luminosity is $\sim
10^{50}$ erg/s. The high Poynting luminosity stage continues for $\sim 0.5$–2
s, and the luminosity subsequently decreases. These properties are consistent
with typical short-hard gamma-ray bursts Nakar (2007); Berger (2014).

For the model with the SFHo EOS as well as for previous models Hayashi  _et
al._ (2022), the Poynting luminosity is likely to drop due to the spreading of
the funnel wall and the decrease of the magnetic-field strength. The spreading
of the funnel wall is caused by the decrease of the gas pressure from the
torus at the funnel wall. However, we also found other two processes that
result in the decrease of the Poynting luminosity in the late stage of the
post-merger evolution. First, for the model with a low initial magnetic-field
strength and the model with no equatorial-plane symmetry, the Poynting
luminosity for a given observer drops due to the disappearance of the
magnetosphere stemming from the reconnection of the magnetic-field lines. This
is caused by the MRI dynamo activity in the accretion disk, which enforces the
magnetic flux with a variety of the polarity to be ejected quasi-periodically
from the disk to the polar region. When strong magnetic fields with the
polarity opposite to that in the magnetosphere emerge from the disk, they
pair-annihilate by the reconnection, and the magnetosphere temporarily
disappears. If the magnetic field emerging from the disk is strong enough by
any chance, it replaces the magnetic field in the magnetosphere,  resulting in
the polarity flip. However, this replacement is only found in the model with
no equatorial-plane symmetry, and the disappearance occurs for the models with
low initial magnetic-field strength and the model with no equatorial plane
symmetry. Second, for the model with the initially toroidal magnetic field in
the neutron star, the Poynting luminosity drops due to the tilt of the
magnetosphere. Because the post-merger ejecta occasionally carries the angular
momentum component not parallel to the black hole spin axis, the accretion
disk is enforced to tilt, in particular in the late-time evolution. Then, it
results in the tilt of the magnetosphere because the funnel structure of the
magnetosphere is determined by the gas pressure from the disk (torus).
Irrespective of these mechanisms, the evolution process of the accretion disk
does determine the evolution process of the magnetosphere. These three
mechanisms all include stochastic processes, and it is not feasible to
precisely predict which mechanisms determine the evolution process of the
magnetosphere. However, these three could be plausible mechanisms to make
short-hard gamma-ray bursts as short as $\sim 0.5$–$2$ s.

###### Acknowledgements.

We thank Sho Fujibayashi, Kunihito Ioka, and Shinya Wanajo for useful
discussions. Numerical simulations were performed on Sakura, Cobra, and Raven
clusters at Max Planck Computing and Data Facility, Yukawa-21 at Yukawa
Institute for Theoretical Physics of Kyoto University, and Cray XC50 at CfCA
of National Astronomical Observatory of Japan. This work was in part supported
by Grant-in-Aid for Scientific Research (grant Nos. 19K14720, 20H00158,
22K03617, and 23H04900) of Japanese MEXT/JSPS. KH was supported by JST SPRING
(grant No. JPMJSP2110).

## References

  * Hayashi  _et al._ (2022) K. Hayashi, S. Fujibayashi, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Phys. Rev. D 106, 023008 (2022).
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016).
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 11, 021053 (2021a).
  * Abbott  _et al._ (2021b) R. Abbott _et al._ , (2021b), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017a).
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ , Astrophys. J. 848, L12 (2017b).
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ , Astrophys. J. 848, L13 (2017c).
  * Abbott  _et al._ (2021c) R. Abbott _et al._ , Astrophys. J. 915, L5 (2021c).
  * Shibata (2015) M. Shibata,  _Numerical Relativity_ , Vol. 1 (World Scientific, 2015).
  * Kyutoku  _et al._ (2021) K. Kyutoku, M. Shibata, and K. Taniguchi, Living Rev. Rel. 24, 5 (2021).
  * Lattimer and Schramm (1974) J. M. Lattimer and D. N. Schramm, Astrophys. J. Lett. 192, L145 (1974).
  * Eichler  _et al._ (1989) D. Eichler, M. Livio, T. Piran, and D. N. Schramm, Nature (London) 340, 126 (1989).
  * Nakar (2007) E. Nakar, Phys. Rep. 442, 166 (2007), the Hans Bethe Centennial Volume 1906-2006.
  * Berger (2014) E. Berger, Annu. Rev. Astron. Astrophys. 52, 43 (2014).
  * Li and Paczyński (1998) L.-X. Li and B. Paczyński, Astrophys. J. 507, L59 (1998).
  * Metzger  _et al._ (2010) B. D. Metzger, G. Martínez-Pinedo, S. Darbha, E. Quataert, A. Arcones, D. Kasen, R. Thomas, P. Nugent, I. V. Panov, and N. T. Zinner, Mon. Not. R. Astron. Soc. 406, 2650 (2010).
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ , Living Rev. Rel. 23, 3 (2020).
  * Gardner  _et al._ (2006) J. P. Gardner _et al._ , Space Sci. Rev. 123, 485 (2006).
  * Ivezić  _et al._ (2019) Ž. Ivezić _et al._ , Astrophys. J. 873, 111 (2019).
  * Shibata and Uryū (2006) M. Shibata and K. Uryū, Phys. Rev. D 74, 121503 (2006).
  * Shibata and Uryū (2007) M. Shibata and K. Uryū, Classical Quantum Gravity 24, S125 (2007).
  * Shibata and Taniguchi (2008) M. Shibata and K. Taniguchi, 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, 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, Phys. Rev. D 78, 104015 (2008).
  * Shibata  _et al._ (2009) M. Shibata, K. Kyutoku, T. Yamamoto, and K. Taniguchi, Phys. Rev. D 79, 044030 (2009).
  * Etienne  _et al._ (2009) Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baumgarte, Phys. Rev. D 79, 044024 (2009).
  * Chawla  _et al._ (2010) S. Chawla, M. Anderson, M. Besselman, L. Lehner, S. L. Liebling, P. M. Motl, and D. Neilsen, 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, Classical Quantum Gravity 27, 114106 (2010).
  * Kyutoku  _et al._ (2010) K. Kyutoku, M. Shibata, and K. Taniguchi, Phys. Rev. D 82, 044049 (2010).
  * Kyutoku  _et al._ (2011) K. Kyutoku, H. Okawa, M. Shibata, and K. Taniguchi, Phys. Rev. D 84, 064018 (2011).
  * Foucart  _et al._ (2011) F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukolsky, Phys. Rev. D 83, 024005 (2011).
  * Foucart  _et al._ (2012) F. Foucart, M. D. Duez, L. E. Kidder, M. A. Scheel, B. Szilagyi, and S. A. Teukolsky, Phys. Rev. D 85, 044015 (2012).
  * Etienne  _et al._ (2012a) Z. B. Etienne, Y. T. Liu, V. Paschalidis, and S. L. Shapiro, Phys. Rev. D 85, 064029 (2012a).
  * Etienne  _et al._ (2012b) Z. B. Etienne, V. Paschalidis, and S. L. Shapiro, Phys. Rev. D 86, 084026 (2012b).
  * Kyutoku  _et al._ (2013) K. Kyutoku, K. Ioka, and M. Shibata, Phys. Rev. D 88, 041503 (2013).
  * Kyutoku  _et al._ (2015) K. Kyutoku, K. Ioka, H. Okawa, M. Shibata, and K. Taniguchi, Phys. Rev. D 92, 044028 (2015).
  * Foucart  _et al._ (2013) 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, Phys. Rev. D 87, 084006 (2013).
  * Lovelace  _et al._ (2013) G. Lovelace, M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Classical Quantum Gravity 30, 135004 (2013).
  * 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, Astrophys. J. 776, 47 (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, Phys. Rev. D 90, 024026 (2014).
  * Paschalidis  _et al._ (2015) V. Paschalidis, M. Ruiz, and S. L. Shapiro, Astrophys. J. 806, L14 (2015).
  * Kawaguchi  _et al._ (2015) K. Kawaguchi, K. Kyutoku, H. Nakano, H. Okawa, M. Shibata, and K. Taniguchi, Phys. Rev. D 92, 024014 (2015).
  * Kiuchi  _et al._ (2015) K. Kiuchi, Y. Sekiguchi, K. Kyutoku, M. Shibata, K. Taniguchi, and T. Wada, Phys. Rev. D 92, 064034 (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, Classical Quantum Gravity 34, 044002 (2017).
  * Kyutoku  _et al._ (2018) K. Kyutoku, K. Kiuchi, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D 97, 023009 (2018).
  * 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, Phys. Rev. D 98, 063009 (2018).
  * Ruiz  _et al._ (2018) M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D 98, 123017 (2018).
  * Foucart  _et al._ (2019a) 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, Phys. Rev. D 99, 044008 (2019a).
  * Foucart  _et al._ (2019b) F. Foucart, M. D. Duez, L. E. Kidder, S. M. Nissanke, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 99, 103025 (2019b).
  * Hinderer  _et al._ (2019) 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, Phys. Rev. D 100, 063021 (2019).
  * Hayashi  _et al._ (2021) K. Hayashi, K. Kawaguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Phys. Rev. D 103, 043007 (2021).
  * Foucart  _et al._ (2021) F. Foucart, A. Chernoglazov, M. Boyle, T. Hinderer, M. Miller, J. Moxon, M. A. Scheel, N. Deppe, M. D. Duez, F. Hébert, L. E. Kidder, W. Throwe, and H. P. Pfeiffer, Phys. Rev. D 103, 064007 (2021).
  * Most  _et al._ (2021) E. R. Most, L. J. Papenfort, S. D. Tootle, and L. Rezzolla, Astrophys. J. 912, 80 (2021).
  * Chaurasia  _et al._ (2021) S. V. Chaurasia, T. Dietrich, and S. Rosswog, Phys. Rev. D 104, 084010 (2021).
  * Most  _et al._ (2021) E. R. Most, L. J. Papenfort, S. D. Tootle, and L. Rezzolla, Mon. Not. R. Astron. Soc. 506, 3511 (2021).
  * Fernández and Metzger (2013) R. Fernández and B. D. Metzger, Mon. Not. R. Astron. Soc. 435, 502 (2013).
  * Metzger and Fernández (2014) B. D. Metzger and R. Fernández, Mon. Not. R. Astron. Soc. 441, 3444 (2014).
  * Just  _et al._ (2015) O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely, and H.-T. Janka, Mon. Not. R. Astron. Soc. 448, 541 (2015).
  * Fernández  _et al._ (2015) R. Fernández, E. Quataert, J. Schwab, D. Kasen, and S. Rosswog, Mon. Not. R. Astron. Soc. 449, 390 (2015).
  * Fernández  _et al._ (2017) R. Fernández, F. Foucart, D. Kasen, J. Lippuner, D. Desai, and L. F. Roberts, Classical Quantum Gravity 34, 154001 (2017).
  * Siegel and Metzger (2018) D. M. Siegel and B. D. Metzger, Astrophys. J. 858, 52 (2018).
  * Fernández  _et al._ (2019) R. Fernández, A. Tchekhovskoy, E. Quataert, F. Foucart, and D. Kasen, Mon. Not. R. Astron. Soc. 482, 3373 (2019).
  * Janiuk (2019) A. Janiuk, Astrophys. J. 882, 163 (2019).
  * Christie  _et al._ (2019) I. M. Christie, A. Lalakos, A. Tchekhovskoy, R. Fernández, F. Foucart, E. Quataert, and D. Kasen, Mon. Not. R. Astron. Soc. 490, 4811 (2019).
  * Miller  _et al._ (2019) J. M. Miller, B. R. Ryan, J. C. Dolence, A. Burrows, C. J. Fontes, C. L. Fryer, O. Korobkin, J. Lippuner, M. R. Mumpower, and R. T. Wollaeger, Phys. Rev. D 100, 023008 (2019).
  * Fujibayashi  _et al._ (2020a) S. Fujibayashi, M. Shibata, S. Wanajo, K. Kiuchi, K. Kyutoku, and Y. Sekiguchi, Phys. Rev. D 101, 083029 (2020a).
  * Fujibayashi  _et al._ (2020b) S. Fujibayashi, M. Shibata, S. Wanajo, K. Kiuchi, K. Kyutoku, and Y. Sekiguchi, Phys. Rev. D 102, 123014 (2020b).
  * Li and Siegel (2021) X. Li and D. M. Siegel, Phys. Rev. Lett. 126, 251101 (2021).
  * Fernández  _et al._ (2020) R. Fernández, F. Foucart, and J. Lippuner, Mon. Not. R. Astron. Soc. 497, 3221 (2020).
  * Just  _et al._ (2022) O. Just, S. Goriely, H. T. Janka, S. Nagataki, and A. Bauswein, Mon. Not. R. Astron. Soc. 509, 1377 (2022).
  * Shibata  _et al._ (2021) M. Shibata, S. Fujibayashi, and Y. Sekiguchi, Phys. Rev. D 104, 063026 (2021).
  * Blandford and Znajek (1977) R. D. Blandford and R. L. Znajek, Mon. Not. R. Astron. Soc. 179, 433 (1977).
  * Balbus and Hawley (1991) S. A. Balbus and J. F. Hawley, Astrophys. J. 376, 214 (1991).
  * Balbus and Hawley (1998) S. A. Balbus and J. F. Hawley, Rev. Mod. Phys. 70, 1 (1998).
  * Banik  _et al._ (2014) S. Banik, M. Hempel, and D. Bandyopadhyay, Astrophys. J. 214, 22 (2014).
  * Steiner  _et al._ (2013) A. W. Steiner, M. Hempel, and T. Fischer, Astrophys. J. 774, 17 (2013).
  * Timmes and Swesty (2000) F. X. Timmes and F. D. Swesty, Astrophys. J. 126, 501 (2000).
  * Miller  _et al._ (2019) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, Astrophys. J. Lett. 887, L24 (2019).
  * Masada and Sano (2008) Y. Masada and T. Sano, Astrophys. J. 689, 1234 (2008).
  * Guilet  _et al._ (2017) J. Guilet, A. Bauswein, O. Just, and H.-T. Janka, Mon. Not. R. Astron. Soc. 471, 1879 (2017).
  * Shakura and Sunyaev (1973) N. Shakura and R. Sunyaev, Astron. Astrophys. 24, 337 (1973).
  * Blandford and Payne (1982) R. D. Blandford and D. G. Payne, Mon. Not. R. Astron. Soc. 199, 883 (1982).
  * De and Siegel (2021) S. De and D. M. Siegel, Astrophys. J. 921, 94 (2021).
  * Wanajo  _et al._ (2022) S. Wanajo, S. Fujibayashi, K. Hayashi, K. Kiuchi, Y. Sekiguchi, and M. Shibata, (2022), arXiv:2212.04507 [astro-ph.HE] .
  * Brandenburg and Subramanian (2005) A. Brandenburg and K. Subramanian, Phys. Rep. 417, 1 (2005).
  * Liska  _et al._ (2020) M. Liska, A. Tchekhovskoy, and E. Quataert, Mon. Not. R. Astron. Soc. 494, 3656 (2020).
  * Tchekhovskoy  _et al._ (2011) A. Tchekhovskoy, R. Narayan, and J. C. McKinney, Mon. Not. R. Astron. Soc. 418, L79 (2011).