<figure><img src="image_1.png" /></figure>

# Magnetorotational instability in neutron star mergers: impact of neutrinos

Jérôme Guilet,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />* Andreas Bauswein,<fnref n="4" /> Oliver Just<fnref n="1" />,<fnref n="2" /> and Hans-Thomas Janka<fnref n="1" />

Accepted 2017 July 10. Received 2017 July 10; in original form 2016 October 26

## ABSTRACT

The merger of two neutron stars may give birth to a long-lived hypermassive neutron star. If it harbours a strong magnetic field of magnetar strength, its spin-down could explain several features of short gamma-ray burst afterglows. The magnetorotational instability (MRI) has been proposed as a mechanism to amplify the magnetic field to the required strength. Previous studies have, however, neglected neutrinos that may have an important impact on the MRI by inducing a viscosity and drag. We investigate the impact of these neutrinos effects on the linear growth of the MRI by applying a local stability analysis to snapshots of a neutron star merger simulation. We find that neutrinos have a significant impact inside the hypermassive neutron star, but have at most a marginal effect in the torus surrounding it. Inside the hypermassive neutron star, the MRI grows in different regimes depending on the radius and on the initial magnetic-field strength. For magnetic fields weaker than $10^{13}-10^{14}$ G, the growth rate of the MRI is significantly reduced due to the presence of neutrinos. We conclude that neutrinos should be taken into account when studying the growth of the MRI from realistic initial magnetic fields. Current numerical simulations, which neglect neutrino viscosity, are only consistent, i.e. in the adopted ideal regime, if they start from artificially strong initial magnetic fields above ~$10^{14}$ G. One should be careful when extrapolating these results to lower initial magnetic fields, where the MRI growth is strongly affected by neutrino viscosity or drag.

**Key words:** MHD – gamma-ray burst: general – stars: magnetars – stars: magnetic field – stars: neutron.

## 1 INTRODUCTION

Neutron star (NS) mergers are the leading candidates to power short gamma-ray bursts (short GRBs; Paczynski 1986; Eichler et al. 1989; Barthelmy et al. 2005) and are one of the most promising targets for the new generation of gravitational wave detectors (LIGO, VIRGO; Abadie et al. 2010). The nature of the compact object formed during the merger and powering short GRBs is debated. It could be either a black hole (BH) or a rapidly rotating, massive, strongly magnetized NS (Usov 1992; Metzger et al. 2011). Such a supermassive magnetar is invoked to explain the extended emission and the X-ray light curves observed in a number of short GRBs (Metzger, Quataert & Thompson 2008; Bucciantini et al. 2012; Gompertz et al. 2013; Rowlinson et al. 2013; Gompertz, O’Brien & Wynn 2014; Lü et al. 2015) and could also lead to late UV, visible and infrared emission (Fan et al. 2013; Metzger & Piro 2014).

The evolution and impact of the magnetic field in NS mergers have been explored through numerical simulations (Anderson

et al. 2008; Liu et al. 2008; Giacomazzo, Rezzolla & Baiotti 2011; Rezzolla et al. 2011; Kiuchi et al. 2014, 2015; Giacomazzo et al. 2015; Palenzuela et al. 2015; Kawamura et al. 2016). An amplification of the magnetic field to magnetar strength in the hypermassive NS (HMNS) has been suggested to occur due to several potential amplification mechanisms. First, the shear layer that forms where the surfaces of the NSs touch is subject to the Kelvin–Helmholtz instability and the turbulence it generates can lead to a small-scale dynamo (e.g. Price & Rosswog 2006; Obergaulinger, Aloy & Müller 2010; Kiuchi et al. 2015). Secondly, the magnetorotational instability (MRI; e.g. Balbus & Hawley 1991) has also been argued to be at work as suggested by numerical simulations (Duez et al. 2006; Siegel et al. 2013; Kiuchi et al. 2014). The wavelength of the MRI, which is proportional to the magnetic field, can only be resolved by current global numerical simulations if the initial magnetic field is already extremely strong ($10^{15}$ G and higher). Observed binary pulsars have a surface dipole magnetic-field strength in the range of $10^{10}-10^{12}$ G (Lorimer 2008). It therefore remains to be investigated whether the growth of the MRI would be significantly different for the weaker initial fields as suggested by these observations. In this paper, we investigate the so-far neglected impact of

<footnote n="1">¹Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany</footnote>

<footnote n="2">²Max Planck/Princeton Center for Plasma Physics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany</footnote>

<footnote n="3">³Laboratoire AIM, CEA/DSM-CNRS-Université Paris 7, Irfu/Service d'Astrophysique, CEA-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴Heidelberger Institut für Theoretische Studien, Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany</footnote>

<footnote>* E-mail: jerome.guilet@cea.fr</footnote>
neutrinos on the MRI in NS mergers. We find that it is significant for realistic initial magnetic fields and that the MRI then grows in regimes different from the high magnetic field case considered so far.

In Section 2, we describe the numerical model of an NS merger that we use. In Section 3, analytical predictions for the growth of the MRI in the presence of neutrino radiation are applied to this numerical model. In Section 4, the impact of neutrinos on the Kelvin-Helmholtz instability is discussed. Finally, in Section 5, we discuss and conclude on the consequences of our results for NS mergers.

## 2 NUMERICAL MODEL

In order to estimate the physical quantities relevant for the growth of the MRI, we use the result of a numerical simulation of an NS merger. The computations were performed with a relativistic smoothed particle hydrodynamics (SPH) code, which imposes the conformal flatness condition on the spatial metric for solving the Einstein equations (Isenberg & Nester 1980; Wilson, Mathews & Marronetti 1996). Details of the numerical and physical model can be found in Oechslin, Janka & Marek (2007) and Bauswein, Oechslin & Janka (2010). For this study, we employed a temperature-dependent, microphysical high-density equation of state, specifically we used the DD2 model (Hempel & Schaffner-Bielich 2010), which is based on the relativistic mean-field parametrization proposed by Typel et al. (2010). This equation of state is moderately stiff and results in an NS radius of 13.2 km for a star with a gravitational mass of $1.35 M_\odot$. We focus on a simulation of a binary with two stars of $1.35 M_\odot$, which can be considered to be representative with respect to the observed distribution of masses in NS binaries (see, e.g. the compilation in Lattimer 2012). This binary set-up leads to the formation of a rapidly rotating NS merger remnant, which is stable until the end of the simulation (several ten milliseconds after merging). The formation of a stable stellar remnant is expected for a large variety of equation of state models in the total binary mass range of about 2.4–3.0 $M_\odot$, which is representative for the observed binaries (Bauswein, Baumgarte & Janka 2013). From the hydrodynamical data of the post-merger object, we extract radial profiles of various quantities by averaging over the azimuthal direction and we restrict our analysis to the equatorial plane for simplicity. Fig. 1 shows the radial profiles at the two times we consider in the analysis: 5 and 15 ms after the collision of the NSs. At the early time, the HMNS is not yet axisymmetric such that this procedure only gives a rough estimate of the conditions locally felt by the MRI. This should be less problematic at the later time we consider. Numerical diffusion inherent to any numerical simulation probably leads to significant uncertainties in particular in two aspects. First, the flattening of the rotation profile in the central 5–10 km of the HMNS observed at the later time may be accelerated by numerical viscosity. Secondly, the temperature is quite uncertain as it is determined by the conversion of kinetic to thermal energy, which is also sensitive to the numerical viscosity. In order to estimate how this uncertainty impacts our results, we also considered temperature profiles 20 per cent higher or lower than the output of the simulation (see Section 3.2). This has a significant quantitative impact on our results but we expect the qualitative picture to remain robust. Similarly, we anticipate that models with different binary masses or a different equation of state lead to quantitative differences but leave the qualitative results unchanged. Moreover, we caution that a merger remnant initially shows a complex velocity field, which significantly de-

<figure><img src="image_2.png" /><figcaption>Figure 1. Azimuthally averaged radial profiles in two snapshots taken 5 ms (black lines) and 15 ms (red lines) after the collision of the NSs. The density (top), temperature (middle) and rotation frequency (bottom) are shown. The exponential fit to the angular frequency profile at the early time is shown with a black dashed line (equation 1). The blue solid line in the bottom panel shows the rotation frequency profile defined by equation (16), while the dashed blue line shows the power-law profile $\Omega = \Omega_1(r/r_1)^{-3/2}$ for comparison. The rotation profile of equation (16) approximates the behaviour of the models with the highest maximum of $\Omega$ in fig. 11 of Hanauske et al. (2016, their models APR4-M125 and APR4-M135).</figcaption></figure>

viates from axisymmetry, while we consider an azimuthally averaged angular velocity and neglect gauge effects as well (Kastaun & Galeazzi 2015).

The somewhat noisy rotation profile at the early time is problematic for the analysis because the results depend critically on the radial derivative of the angular frequency. To avoid this problem, we replace the numerical profile by the following fit of the angular frequency $\Omega$ as a function of radius $r$ (dashed black line in Fig. 1):

$$ \Omega = \Omega_0 \exp(-r/r_0), \quad (1) $$

with $\Omega_0 = 1.12 \times 10^3 \text{ s}^{-1}$ and $r_0 = 20 \text{ km}$. The relation of these angular momentum profiles to results in the literature and the influence on our analysis will be discussed in Section 3.4.
## 3 MRI GROWTH IN THE PRESENCE OF NEUTRINO RADIATION

In this section, we apply analytical estimates obtained in Guilet, Müller & Janka (2015) for the growth of the MRI in the presence of neutrino radiation. These analytical results were obtained under a number of assumptions, which we quickly summarize (see Guilet et al. 2015, for a more detailed discussion). We use a local dispersion relation, which is valid if the wavelength of the modes is much shorter than the typical length-scales of the object. This is well justified for moderate magnetic fields (see Section 3.3). For simplicity, the magnetic field is assumed to be perpendicular to the equatorial plane, which is the most classical and favourable configuration for MRI growth. Only axisymmetric modes are considered as these are the most unstable ones under this magnetic-field geometry. The resistivity is neglected because it is expected to be extremely small (Thompson & Duncan 1993). We also neglect the impact of the (presumably stable) stratification connected to entropy and composition gradients. A stable stratification could in principle have a significantly stabilizing influence on the MRI, but neutrino diffusion is expected to alleviate this impact (Menou, Balbus & Spruit 2004; Masada, Sano & Shibata 2007; Guilet & Müller 2015). We chose to neglect it for simplicity because the composition and entropy gradients are not perfectly axisymmetric and quite uncertain (see Section 2), which hinders a robust conclusion on their impact on the MRI. Finally, all the calculations are non-relativistic because no relativistic linear analysis of the MRI including neutrinos has been performed so far. This is most likely not a major limitation, because general relativistic analyses in ideal MHD suggest that the MRI is not critically changed with respect to a Newtonian approximation even for a Kerr metric (Araya-Góchez 2002; Gammie 2004; Yokosawa & Inui 2005).

### 3.1 MRI with neutrino viscosity (diffusive regime)

On length-scales longer than the mean free path of neutrinos, neutrino diffusion induces a viscosity. We estimate it by using the approximate analytical expression obtained by Keil, Janka & Müller (1996),

$$v = 1.2 \times 10^{10} \left( \frac{T}{10 \text{ MeV}} \right)^2 \left( \frac{\rho}{10^{13} \text{ g cm}^{-3}} \right)^{-2} \text{cm}^2 \text{s}^{-1}, \quad (2)$$

where v is the kinematic viscosity, ρ is the density and T is the temperature. This expression was successfully compared to a direct calculation of the neutrino viscosity using the neutrino distribution from a numerical simulation of protoneutron star formation with elaborate neutrino transport (Guilet et al. 2015). The viscosity we estimate from the simulation outputs is shown in the upper panel of Fig. 2; it increases strongly with radius from ~10^7 cm^2 s⁻<fnref n="1" /> in the inner core of the HMNS to ~10^11–10^12 cm^2 s⁻<fnref n="1" /> near the surface of the HMNS.<fnref n="1" />

The effect of viscosity on the linear growth of the MRI is controlled by the viscous Elsasser number $E_ν \equiv \frac{v_\Lambda}{v_\Omega}$ (e.g. Pessah & Chan 2008; Longaretti & Lesur 2010; Guilet et al. 2015), where $v_\Lambda \equiv B/\sqrt{4\pi\rho}$ is the Alfvén speed with the magnetic-field strength B. For $E_ν < 1$, viscosity significantly affects the growth of the MRI; as a result the growth rate is decreased, and the wavelength of

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Top panel: radial profile of neutrino viscosity in the equatorial plane. Bottom panel: wavelength of the most unstable MRI mode in the viscous regime computed using equation (7) (dashed line) compared to the mean free path of heavy lepton neutrinos (solid lines). The two snapshots are shown in black (5 ms after the collision of the NSs) and red (15 ms after the collision of the NSs).</figcaption></figure>

the most unstable mode becomes longer. Typical conditions inside the HMNS lead to the following estimate of the viscous Elsasser number at a radius of ~5–10 km

$$E_ν \sim 8 \times 10^{-4} \left( \frac{B}{10^{12} \text{ G}} \right)^2 \left( \frac{\rho}{4 \times 10^{14} \text{ g cm}^{-3}} \right)^{-1} \times \left( \frac{\Omega}{6000 \text{ s}^{-1}} \right)^{-1} \left( \frac{v}{4 \times 10^7 \text{ cm}^2 \text{ s}^{-1}} \right)^{-1}. \quad (3)$$

Viscosity therefore has a large impact on the linear growth of the MRI, unless the magnetic field is significantly stronger than the surface field of normal pulsars or another mechanism significantly amplifies an initially weak field. The critical strength of the magnetic field below which viscous effects become important (at which $E_ν = 1$) is

$$\begin{aligned} B_{\text{visc}} &= \sqrt{4\pi\rho\nu\Omega} \\ &= 3.5 \times 10^{13} \left( \frac{\rho}{4 \times 10^{14} \text{ g cm}^{-3}} \right)^{1/2} \left( \frac{v}{4 \times 10^7 \text{ cm}^2 \text{ s}^{-1}} \right)^{1/2} \\ &\quad \times \left( \frac{\Omega}{6000 \text{ s}^{-1}} \right)^{1/2} \text{G}. \end{aligned} \quad (4)$$

The dispersion relation governing the growth of the MRI in the presence of viscosity reads (Lesur & Longaretti 2007; Masada & Sano 2008; Pessah & Chan 2008)

$$[(\sigma + k^2 v)\sigma + k^2 v_\Lambda^2]^2 + \kappa^2 (\sigma^2 + k^2 v_\Lambda^2) - 4\Omega^2 k^2 v_\Lambda^2 = 0, \quad (5)$$

<footnote n="1">¹ We define the surface of the HMNS as the location beyond which the rotational support exceeds the pressure support. It is located at a radius of about 20 km.</footnote>
where $\sigma$ is the growth rate, $k$ is the wavenumber and $\kappa$ is the epicyclic frequency defined by $\kappa^2 = \frac{1}{r^3} \frac{d(r^4\Omega^2)}{dr}$.

In the asymptotic limit $E_v \ll 1$, the growth rate and wavelength of the most unstable mode can then be expressed as (Masada & Sano 2008; Pessah & Chan 2008; Masada et al. 2012)

$$ \sigma = \left( \frac{qE_v}{\tilde{\kappa}} \right)^{1/2} \Omega = 90 \left( \frac{B}{10^{12} \text{ G}} \right) \left( \frac{\rho}{4 \times 10^{14} \text{ g cm}^{-3}} \right)^{-1/2} \times \left( \frac{\nu}{4 \times 10^7 \text{ cm}^2 \text{ s}^{-1}} \right)^{-1/2} \left( \frac{\Omega}{6000 \text{ s}^{-1}} \right)^{1/2} \text{ s}^{-1}, \quad (6) $$

$$ \lambda = 2\pi \left(\frac{\nu}{\kappa}\right)^{1/2} = 390 \left(\frac{\kappa}{6000 \text{ s}^{-1}}\right)^{-1/2} \left(\frac{\nu}{4 \times 10^7 \text{ cm}^2 \text{ s}^{-1}}\right)^{1/2} \text{ cm}, \quad (7) $$

where $\tilde{\kappa} = \kappa/\Omega = \sqrt{2(2-q)}$ is the dimensionless epicyclic frequency and $q = -d \log \Omega/d \log r = 0.5$ in the numerical estimates. In contrast to the inviscid case, the wavelength of the fastest growing mode is independent of the magnetic-field strength and it is longer than it would be without viscosity. The growth rate on the other hand is reduced and is proportional to the magnetic-field strength. The MRI therefore requires a minimum magnetic-field strength in order to grow on a given time-scale. The minimum magnetic field necessary for the MRI to grow at a minimum growth rate $\sigma_{min}$ is (Guilet et al. 2015)

$$ \begin{aligned} B_{\min} &= \left( \frac{4\pi\rho\tilde{\kappa}\nu}{q\Omega} \right)^{1/2} \sigma_{\min}, & (8) \\ &= 10^{12} \left( \frac{\sigma_{\min}}{100 \text{ s}^{-1}} \right) \left( \frac{\rho}{4 \times 10^{14} \text{ g cm}^{-3}} \right)^{1/2} \\ & \times \left( \frac{\nu}{4 \times 10^7 \text{ cm}^2 \text{ s}^{-1}} \right)^{1/2} \left( \frac{\Omega}{6000 \text{ s}^{-1}} \right)^{-1/2} \text{ G}. & (9) \end{aligned} $$

The viscous regime is valid if the wavelength of the fastest growing mode is larger than the neutrino mean free path. Fig. 2 (lower panel) shows that this is indeed the case inside the HMNS but not in the surrounding torus. In that figure, the mean free path of heavy lepton neutrinos (which give the most restrictive constraint among the neutrino species) is estimated using the scaling relation $l_v = 10^4 (\rho/10^{13} \text{ g cm}^{-3})^{-1} (T/10 \text{ MeV})^{-2} \text{ cm}$ (Guilet et al. 2015).

## 3.2 MRI with neutrino drag

At wavelengths shorter than their mean free path, neutrinos induce a drag force of the form $-\Gamma\nu$, where $\Gamma$ is a damping rate depending on the neutrino distribution and opacity (Jedamzik, Katalinić & Olinto 1998; Thompson, Quataert & Burrows 2005; Guilet, Müller & Janka 2015). We compute the neutrino drag damping rate from the simulation outputs using the scaling obtained by Guilet et al. (2015)

$$ \Gamma = 6 \times 10^3 (T/10 \text{ MeV})^6 \text{ s}^{-1}. \quad (10) $$

The dispersion relation governing the growth of the MRI in the presence of neutrino drag is (Guilet et al. 2015)

$$ [(\sigma + \Gamma)\sigma + k^2 v_A^2]^2 + \kappa^2 (\sigma^2 + k^2 v_A^2) - 4\Omega^2 k^2 v_A^2 = 0. \quad (11) $$

The neutrino drag significantly impacts the MRI when $\Gamma \gtrsim \Omega$. Fig. 3 shows that this is the case at all radii inside the HMNS ($r \lesssim 20 \text{ km}$) but not in the surrounding torus ($r \gtrsim 20 \text{ km}$). In the asymptotic limit

<figure><img src="image_4.png" /><figcaption>Figure 3. Top panel: radial profile of the neutrino drag damping rate (solid lines) compared to the angular frequency (dashed lines). Bottom panel: growth rate of the MRI including (solid lines) or neglecting (dashed lines) the impact of neutrinos. The black and red lines correspond to the early and late times, respectively.</figcaption></figure>

$\Gamma \gg \Omega$, the growth rate and wavenumber of the fastest growing mode are (Guilet et al. 2015)

$$ \sigma = \frac{q \Omega^2}{2 \Gamma}, \quad (12) $$

and

$$ k = \sqrt{q/2} \frac{\Omega}{v_A}, \quad (13) $$

respectively. Compared to the ideal MHD case, the maximum growth rate is reduced by a factor $\Gamma/\Omega$ while the wavelength of the fastest growing mode is not drastically changed. Importantly, and contrary to the viscous regime, the growth rate is independent of the magnetic-field strength. The growth rate taking into account the neutrino drag is compared to the ideal MHD case in the bottom panel of Fig. 3. As expected from the above arguments, the growth of the MRI is significantly slowed down by the neutrino drag inside the HMNS. Because of the steep dependence of the neutrino drag damping rate on the temperature, the growth rate in this regime is somewhat uncertain. This is illustrated by Fig. 4 where the impact of varying the temperature profile by 20 per cent is shown to have quantitatively significant consequences but without changing the qualitative picture that neutrino drag reduces the MRI growth rate significantly inside the HMNS but impacts it at most marginally in the surrounding torus.

The neutrino drag regime is valid if the wavelength of the fastest growing MRI mode is shorter than the neutrino mean free path.
<figure><img src="image_5.png" /><figcaption>Figure 4. Growth rate of the MRI in the drag regime at early time ($t = 5$ ms) when the temperature is scaled up (red line) or down (blue line) by 20 per cent. The temperature uncertainty has a significant quantitative impact but leaves the qualitative conclusions unaffected.</figcaption></figure>

Using equation (13), this condition can be expressed as an upper limit on the magnetic-field strength,

$$B < \sqrt{\rho q / 2\pi\Omega l_v} \quad (14)$$

This critical magnetic field is shown as a dashed line in Figs 5 and 6 with the mean free path of electron neutrinos estimated as $l_v = 2 \times 10^3 (\rho/10^{13} \text{ g cm}^{-3})^{-1} (T/10 \text{ MeV})^{-2} \text{ cm}$.

### 3.3 Growth rate and wavelength

In this section, we determine the relevant regime of MRI growth as a function of radius and magnetic-field strength. This then allows us to compute the growth rate and wavelength of the most unstable mode (Figs 5 and 6). At any radius, both the viscous and the drag regimes can in principle be relevant, and one has to compare the growth rate in each regime in order to obtain the fastest growing mode. Since the growth rate in the viscous regime is proportional to the magnetic-field strength, the growth is faster in the viscous regime (dark blue in top panel of Figs 5 and 6) than in the drag regime (light blue) above a critical magnetic-field strength, which can be obtained by combining equations (6) and (12)

$$B_{\text{visc-drag}} = \sqrt{q\pi\rho v\kappa \frac{\Omega}{\Gamma}} \quad (15)$$

This expression is valid when $E_v \ll 1$ and $\Gamma/\Omega \gg 1$ and the result is shown by the solid line going from the bottom left to upper right in the middle and bottom panels of Figs 5 and 6. There is also a mixed regime of MRI growth (green) where some neutrino species are in the viscous regime while some are in the drag regime. There is no analytical description of this regime so far and since it is quite limited in parameter space, we do not attempt to derive one. Instead, we adopt the practical and very rough approach of using the average of the values obtained in the viscous and drag regimes for the growth rate and wavelength of the most unstable mode.

The growth rate and wavelength of the fastest growing modes are obtained by numerically solving the dispersion relation in each relevant regime (equations 5 and 11). The result is shown as a function of radius and magnetic-field strength in Fig. 5 for the early snapshot (5 ms after the collision of the NSs) and Fig. 6 for the late snapshot (15 ms). The MRI can grow unabated on time-scales of 1 ms or shorter for magnetic fields stronger than ~$10^{14}$ G inside

<figure><img src="image_6.png" /><figcaption>Figure 5. Properties of the MRI in a snapshot 5 ms after the collision of the NSs. Top panel: different regimes of MRI growth as a function of radius and magnetic-field strength. The regime labelled slow corresponds to growth time-scales longer than 10 ms. Middle panel: growth rate of the MRI. Bottom panel: wavelength of the most unstable MRI mode. The vertical dotted lines show the boundary between the pressure supported HMNS and the rotation supported torus surrounding it. The lines separating the different regimes are given by equation (4) (giving the limit between viscous and ideal regimes), equation (14) (delineating the boundary between ‘drag’ and ‘mixed’ regimes, and shown as a dashed line in middle and bottom panels) and equation (15) (solid line going from bottom left to top right). The vertical solid line separating the drag and ideal regimes shows the radius at which $\Gamma = \Omega$.</figcaption></figure>

the HMNS and for any magnetic field in the surrounding torus. Inside the HMNS and for weaker magnetic fields, the MRI growth is slowed down by neutrinos in either of the two regimes: viscous or drag. The growth time-scale is longer than 10 ms for densities roughly above $10^{14} \text{ g cm}^{-3}$ and magnetic fields smaller than $10^{12} \text{ G}$. In the viscous regime, the wavelength of the fastest growing mode is of the order of a few metres, but can be much shorter in the drag regime (<1 cm). These short length-scales highlight the difficulty of resolving the MRI in a global numerical simulation. Our analysis also stresses the fact that if one avoids this problem by starting with a strong initial magnetic field, the MRI then grows in a very different physical regime.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Same as Fig. 5 but at a later time, 15 ms after the collision of the NSs. The range of radii has been limited to $r > 6$ km to exclude the MRI-stable inner core in solid-body rotation.</figcaption></figure>

## 3.4 Dependence on the angular velocity profile

It is important to compare our rotation profiles to the results published for simulations performed with grid-based codes (Shibata, Taniguchi & Uryū 2005; Shibata & Taniguchi 2006; Kastaun & Galeazzi 2015; Endrizzi et al. 2016; Hanauske et al. 2016; Kastaun, Ciolfi & Giacomazzo 2016a; Kastaun et al. 2016b). Overall, these simulations suggest a rotation frequency profile that has a maximum at the centre at early times (for a few milliseconds after the merger). This profile then evolves to one with a slowly rotating core, where the rotation frequency increases with the radius up to a maximum [in the range $\Omega \sim (6-10) \times 10^3 \text{ s}^{-1}$ at a radius of $r \sim 5-15$ km depending on the simulations] and decreases outwards. Our early-time profile is qualitatively consistent with these results at early times. Our late-time profile, on the other hand, is flat in the centre with a maximum value of $\Omega \simeq 7.5 \times 10^3 \text{ s}^{-1}$ up to a radius of $r \simeq 10$ km and then decreases outwards again. This is different to most published grid-based simulations but similar to one simulation published by Shibata & Taniguchi (2006) and recent simulations including a viscous modelling of turbulence (Radice 2017; Shibata & Kiuchi 2017). A plausible explanation of the similarities and dissimilarities of our SPH results and published $\Omega$-profiles may be related to a larger numerical viscosity of our SPH simulations and

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as Fig. 6 but with the rotation profile as defined by equation (16). The range of radii has been limited to $r > 7$ km to exclude the MRI-stable inner core where $\Omega$ increases with the radius.</figcaption></figure>

of the early grid-based simulations. In any case, we stress that this difference is restricted to the inner core, which is anyway stable to the MRI in either way (whether $\Omega$ is flat or increasing). For our analysis of MRI growth, the only relevant region is therefore the outer part of the HMNS, where the frequency decreases outwards. In this region, our late-time rotation profile is qualitatively but also quantitatively in the ballpark of the published results (compare to fig. 11 of Hanauske et al. 2016). In order to determine the sensitivity of our results to the rotation profile, we introduce a rotation profile alternative to the ones used in Sections 3.1–3.3 (equation 1 and our late-time profile). This profile (shown by the solid blue line in Fig. 1) has the following functional form:

$$ \Omega(r) = \Omega_1 \frac{(r/r_1)^{-3/2}}{1 + (r/r_{\text{max}})^{-3}}, \quad (16) $$

with $\Omega_1 = 2\pi \times 10^3 \text{ s}^{-1}$, $r_1 = 15$ km and $r_{\text{max}} = 7$ km. It approximates the behaviour of the models with the highest maximum of $\Omega$ in fig. 11 of Hanauske et al. (2016, their models APR4-M125 and APR4-M135). We repeated our analysis of the late-time profile but with this alternative rotation law and show the result in Fig. 7. Compared to Fig. 6, the regime of slow MRI is somewhat less extended due to the larger values of $\Omega$. Overall, the differences are none the less small and the qualitative picture remains unchanged.
Overall, our set of three rotation velocity profiles is representative of the full range of profiles observed in fig. 11 of Hanauske et al. (2016): a low-Ω curve like their model GNH3-M135 is close to our exponential fit at radii larger than 5–10 km, the profile given by equation (16) is close to high-Ω profiles (their models APR4-M125 and APR4-M135), and finally our late-time profile lies in-between. Our results demonstrate that the differences in the Ω(r) behaviour do not have a crucial influence on our main conclusions, the differences in the MRI growth rate and wavelength being noticeable but not major overall (compare Figs 5–7).

It is more difficult to compare the temperature profiles since the aforementioned papers have not published any azimuthally averaged profiles, but we can make the following discussion based on the colour figures of Hanauske et al. (2016) and Kastaun et al. (2016b). Similarly to the rotation frequency profile, the late-time temperature profile shows a low-temperature core surrounded by larger temperatures, the maximum temperature being reached at a radius close to the maximum of the rotation frequency. The inner core of these simulations is significantly cooler than ours, but again this difference is irrelevant to our analysis as the inner core is stable to the MRI. On the other hand, the maximum temperature (with values of $T \sim 30$ MeV) seems somewhat larger than ours. If this higher temperature is more realistic than ours, it would only strengthen our result that neutrinos slow down the MRI growth (see Section 3.2 for a discussion of the impact of the temperature uncertainty).

## 3.5 Accretion torus around BH

We have also applied the same analysis to models of an accretion torus around a BH, which are representative of the remnants expected from NS–BH mergers or from NS–NS mergers leading to prompt or delayed BH formation. The numerical models were evolved by Just et al. (2015a), using an M1-type neutrino transport scheme (Just, Obergaulinger & Janka 2015b), with Newtonian dynamics but including a pseudo-Newtonian gravitational potential (Artemova, Bjoernsson & Novikov 1996), and by employing the $\alpha$-viscosity treatment by Shakura & Sunyaev (1973) to parametrize turbulent angular momentum transport. More details regarding the evolution scheme and model set-up can be found in Just et al. (2015a). For this study, we consider models M3A8m1a2, M3A8m03a2 and M4A8m3a2. The two former models, which were initialized with BH masses of $M_{BH} = 3 M_{\odot}$ and torus masses of $m_{torus} = 0.1$ and $0.03 M_{\odot}$, respectively, are representative for remnants of binary NS mergers. Model M4A8m3a2 with $M_{BH} = 4 M_{\odot}$ and $m_{torus} = 0.3 M_{\odot}$ represents the remnant of a BH–NS merger. We note that $M_{BH}$ in this model is rather low in view of the expected BH mass distribution (e.g. Dominik et al. 2012). A less massive and thus smaller BH typically leads to higher densities and temperatures in the torus. As a consequence, in our model the (already small) impact of neutrinos on the MRI tends to be overestimated with respect to more likely cases of higher BH mass. For all three models, the $\alpha$-parameter for the viscosity is 0.02, the BH spin parameter is 0.8 and the considered snapshot is 20 ms after the start of the simulations.

In all three models, the neutrino mean free path is comparable to the vertical thickness of the torus (i.e. the optical depth is of order unity) and as a result the drag regime is the relevant one. We find that the neutrino drag is negligible ($\Gamma \ll \Omega$) in models of low or moderate mass of the accretion torus (Fig. 8, black and blue lines) and has at most a marginal impact ($\Gamma \lesssim \Omega$) in the case of a massive torus (Fig. 8, red line). This confirms the results of Foucart et al. (2015).

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Neutrino drag damping rate in three models of post-merger BH accretion tori (solid lines). The black and blue lines correspond to tori of masses 0.03 and $0.1 M_{\odot}$, respectively, with a $3 M_{\odot}$ BH, while the red lines correspond to a torus of mass $0.3 M_{\odot}$ with a $4 M_{\odot}$ BH. The dashed lines show the corresponding angular frequencies (with the black line being nearly identical to the blue line).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of the neutrino damping rate (red line) with the growth rate of the Kelvin–Helmholtz instability when neutrinos are neglected (black line). The Kelvin–Helmholtz instability is not significantly modified by neutrinos since the neutrino damping rate is always smaller than the growth rate. The neutrino damping rate shown here corresponds to $\rho = 10^{14} \text{ g cm}^{-3}$ and $T = 30$ MeV.</figcaption></figure>

## 4 KELVIN–HELMHOLTZ INSTABILITY IN THE PRESENCE OF NEUTRINOS

In this section, we give a rough estimate of the impact of neutrinos on the Kelvin–Helmholtz instability that develops at the contact interface between the two NSs. To determine whether neutrinos have a significant impact, we compare the growth rate of the instability (black line in Fig. 9) to the damping rate due to neutrinos (red line in Fig. 9) at different wavelengths. The damping rate due to neutrinos is $k^2\nu$ in the viscous regime (where the viscosity $\nu$ is given by equation 2), while it is $\Gamma$ as given by equation (10) in the drag regime. The transition between the two regimes takes place at a wavelength close to the neutrino mean free path. Initially, the shear layer is very thin and we will first assume that it can be approximated by a discontinuity with a total velocity jump $\Delta\nu \sim 6 \times 10^9 \text{ cm s}^{-1}$ (see fig. 5 of Kiuchi et al. 2015). The growth rate of the Kelvin–Helmholtz instability can then be written as $k\Delta\nu/2$ (black line in Fig. 9). Fig. 9 shows, for physical conditions typical of NS mergers, that the damping rate due to neutrinos is significantly smaller than the growth rate of the Kelvin–Helmholtz instability for all
wavelengths. It can be shown (see the appendix) that the criterion for neutrinos to have an impact requires a high temperature and a low density, which are unlikely to be obtained:

$$ \left(\frac{T}{30 \text{ MeV}}\right)^8 \left(\frac{\rho}{10^{14} \text{ g cm}^{-3}}\right)^{-2} > 25. \qquad (17) $$

As time passes, the width $\Delta x$ of the shear layer is expected to increase. The impact of the smoothness of the shear layer on the instability is to stabilize modes with a wavelength shorter than $\Delta x$. The wavelength and growth rate of the fastest growing mode are then a fraction of $\Delta x$ and of the maximum vorticity $\sim \Delta v/\Delta x$, respectively, and therefore are close to the idealized curve plotted on Fig. 9. We therefore expect our conclusion to remain valid for a smooth velocity profile.

We conclude that the neutrinos do not change the growth rate and length-scale of the Kelvin–Helmholtz instability. This stands in contrast to the MRI, the difference being attributable to the shorter growth time-scale of the Kelvin–Helmholtz instability.

## 5 DISCUSSION AND CONCLUSION

We have studied the growth of the MRI in compact object mergers, taking into account the impact of neutrino viscosity and drag. The main findings can be summarized as follows.

(i) Inside the HMNS, the growth rate of the MRI is significantly reduced by neutrino effects if the magnetic-field strength is lower than $\sim 10^{13}-10^{14}$ G. Depending on the radius and the magnetic field, the MRI grows in different regimes: the drag regime at wavelengths shorter than the neutrino mean free paths for weak magnetic fields, the viscous regime at wavelengths longer than the mean free path for stronger magnetic fields or a mixed regime at intermediate fields (see Figs 5 and 6). The growth time-scale becomes longer than 10 ms for densities above roughly $10^{14} \text{ g cm}^{-3}$ and magnetic fields weaker than $10^{12}$ G.

(ii) In the rotationally supported torus around the HMNS, the MRI is at most marginally affected by neutrinos. This conclusion also applies to tori in cases of a central BH, thus confirming the results of Foucart et al. (2015).

(iii) The growth rate of the Kelvin–Helmholtz instability at the contact interface between the two NSs is not significantly influenced by neutrinos.

Even though we believe that these results are qualitatively robust, we caution that our analysis suffers from significant quantitative uncertainties. In particular, while we consider an azimuthally averaged angular velocity field, the merger remnant initially shows a more complex flow pattern, which significantly deviates from axisymmetry.

Our analysis demonstrates that inside the HMNS and for initial magnetic fields as expected for typical NSs, the MRI is significantly impacted by neutrinos. We stress that current global numerical simulations, starting with an unrealistically strong magnetic field, are in a different physical regime of MRI growth. Results within that regime can therefore not easily be extrapolated to lower initial magnetic fields. Given the extreme difficulties in numerically resolving the small length-scales at which the MRI grows, local simulations are an indispensable complementary method to study the non-linear phase of the MRI in the presence of neutrinos. The impact of neutrino drag on the non-linear MRI dynamics, for example, is unknown and could be studied with such an approach.

The slower growth of the MRI that we predict inside the HMNS is likely to lead to a delay of several tens of milliseconds between

the collision of binary NSs and the formation of a magnetic field of magnetar strength inside the HMNS. In the scenario where the central engine of a short GRB is a magnetar, we speculate that this could contribute to a delay between the emission of gravitational waves and the formation of the relativistic jet responsible for the GRB. Alternatively, the slow MRI growth may imply that another mechanism such as the Kelvin–Helmholtz instability at the NS–NS collision interface is a more efficient agent during the first phase of magnetic-field amplification. The viability of this mechanism has been debated (Obergaulinger et al. 2010; Giacomazzo et al. 2011; Dionysopoulou, Alic & Rezzolla 2015) but recent high-resolution numerical simulations support its efficiency (Kiuchi et al. 2015). It remains to be clarified whether the magnetic field can pervade the whole volume of the HMNS and be coherent on large scales as is implicitly assumed in models of magnetar-powered GRBs.

## ACKNOWLEDGEMENTS

At Garching, this work was supported by the Max-Planck–Princeton Center for Plasma Physics (MPPC), by the Deutsche Forschungs-gemeinschaft through Excellence Cluster ‘Universe’ (EXC 153) and by the European Research Council through grant ERC-AdG no. 341157-COCO2CASA. AB acknowledges support by the Klaus Tschira Foundation. JG acknowledges support from the European Research Council (grant MagBURST number 715368).

## REFERENCES

Abadie J. et al., 2010, Class. Quantum Grav., 27, 173001

Anderson M., Hirschmann E. W., Lehner L., Liebling S. L., Motl P. M., Neilsen D., Palenzuela C., Tohline J. E., 2008, Phys. Rev. Lett., 100, 191101

Araya-Góchez R. A., 2002, MNRAS, 337, 795

Artemova I. V., Bjoernsson G., Novikov I. D., 1996, ApJ, 461, 565

Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214

Barthelmy S. D. et al., 2005, Nature, 438, 994

Bauswein A., Oechslin R., Janka H.-T., 2010, Phys. Rev. D, 81, 024012

Bauswein A., Baumgarte T. W., Janka H.-T., 2013, Phys. Rev. Lett., 111, 131101

Bucciantini N., Metzger B. D., Thompson T. A., Quataert E., 2012, MNRAS, 419, 1537

Dionysopoulou K., Alic D., Rezzolla L., 2015, Phys. Rev. D, 92, 084064

Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, ApJ, 759, 52

Duez M. D., Liu Y. T., Shapiro S. L., Shibata M., Stephens B. C., 2006, Phys. Rev. Lett., 96, 031101

Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126

Endrizzi A., Ciolfi R., Giacomazzo B., Kastaun W., Kawamura T., 2016, Class. Quantum Grav., 33, 164001

Fan Y.-Z., Yu Y.-W., Xu D., Jin Z.-P., Wu X.-F., Wei D.-M., Zhang B., 2013, ApJ, 779, L25

Foucart F. et al., 2015, Phys. Rev. D, 91, 124021

Gammie C. F., 2004, ApJ, 614, 309

Giacomazzo B., Rezzolla L., Baiotti L., 2011, Phys. Rev. D, 83, 044014

Giacomazzo B., Zrake J., Duffell P. C., MacFadyen A. I., Perna R., 2015, ApJ, 809, 39

Gompertz B. P., O’Brien P. T., Wynn G. A., Rowlinson A., 2013, MNRAS, 431, 1745

Gompertz B. P., O’Brien P. T., Wynn G. A., 2014, MNRAS, 438, 240

Guilet J., Müller E., 2015, MNRAS, 450, 2153

Guilet J., Müller E., Janka H.-T., 2015, MNRAS, 447, 3992

Hanauske M., Takami K., Bovard L., Rezzolla L., Font J. A., Galeazzi F., Stöcker H., 2016, preprint (arXiv:e-prints)

Hempel M., Schaffner-Bielich J., 2010, Nucl. Phys. A, 837, 210
Isenberg J., Nester J., 1980, in Held A., ed., General Relativity and Gravitation, Vol. 1. Plenum Press, New York, p. 23

Jedamzik K., Katalinić V., Olinto A. V., 1998, Phys. Rev. D, 57, 3264

Just O., Bauswein A., Pulpillo R. A., Goriely S., Janka H.-T., 2015a, MNRAS, 448, 541

Just O., Obergaulinger M., Janka H.-T., 2015b, MNRAS, 453, 3386

Kastaun W., Galeazzi F., 2015, Phys. Rev. D, 91, 064027

Kastaun W., Ciolfi R., Giacomazzo B., 2016a, Phys. Rev. D, 94, 044060

Kastaun W., Ciolfi R., Endrizzi A., Giacomazzo B., 2016b, preprint (arXiv:e-prints)

Kawamura T., Giacomazzo B., Kastaun W., Ciolfi R., Endrizzi A., Baiotti L., Perna R., 2016, Phys. Rev. D, 94, 064012

Keil W., Janka H., Müller E., 1996, ApJ, 473, L111

Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., Wada T., 2014, Phys. Rev. D, 90, 041502

Kiuchi K., Cerdá-Durán P., Kyutoku K., Sekiguchi Y., Shibata M., 2015, Phys. Rev. D, 92, 124034

Lattimer J. M., 2012, Ann. Rev. Nucl. Part. Sci., 62, 485

Lesur G., Longaretti P.-Y., 2007, MNRAS, 378, 1471

Liu Y. T., Shapiro S. L., Etienne Z. B., Taniguchi K., 2008, Phys. Rev. D, 78, 024012

Longaretti P.-Y., Lesur G., 2010, A&A, 516, A51

Lorimer D. R., 2008, Living Rev. Relativ., 11

Lü H.-J., Zhang B., Lei W.-H., Li Y., Lasky P. D., 2015, ApJ, 805, 89

Masada Y., Sano T., 2008, ApJ, 689, 1234

Masada Y., Sano T., Shibata K., 2007, ApJ, 655, 447

Masada Y., Takiwaki T., Kotake K., Sano T., 2012, ApJ, 759, 110

Menou K., Balbus S. A., Spruit H. C., 2004, ApJ, 607, 564

Metzger B. D., Piro A. L., 2014, MNRAS, 439, 3916

Metzger B. D., Quataert E., Thompson T. A., 2008, MNRAS, 385, 1455

Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., Quataert E., 2011, MNRAS, 413, 2031

Obergaulinger M., Aloy M. A., Müller E., 2010, A&A, 515, A30

Oechslin R., Janka H.-T., Marek A., 2007, A&A, 467, 395

Paczynski B., 1986, ApJ, 308, L43

Palenzuela C., Liebling S. L., Neilsen D., Lehner L., Caballero O. L., O'Connor E., Anderson M., 2015, Phys. Rev. D, 92, 044045

Pessah M. E., Chan C.-k., 2008, ApJ, 684, 498

Price D. J., Rosswog S., 2006, Science, 312, 719

Radice D., 2017, ApJ, 838, L2

Rezzolla L., Giacomazzo B., Baiotti L., Granot J., Kouveliotou C., Aloy M. A., 2011, ApJ, 732, L6

Rowlinson A., O'Brien P. T., Metzger B. D., Tanvir N. R., Levan A. J., 2013, MNRAS, 430, 1061

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shibata M., Kiuchi K., 2017, preprint (arXiv:e-prints)

Shibata M., Taniguchi K., 2006, Phys. Rev. D, 73, 064027

Shibata M., Taniguchi K., Uryü K., 2005, Phys. Rev. D, 71, 084021

Siegel D. M., Ciolfi R., Harte A. I., Rezzolla L., 2013, Phys. Rev. D, 87, 121302

Thompson C., Duncan R. C., 1993, ApJ, 408, 194

Thompson T. A., Quataert E., Burrows A., 2005, ApJ, 620, 861

Typel S., Röpke G., Klähn T., Blaschke D., Wolter H. H., 2010, Phys. Rev. C, 81, 015803

Usov V. V., 1992, Nature, 357, 472

Wilson J. R., Mathews G. J., Marronetti P., 1996, Phys. Rev. D, 54, 1317

Yokosawa M., Inui T., 2005, ApJ, 631, 1051

## APPENDIX: KELVIN-HELMHOLTZ INSTABILITY IN THE PRESENCE OF NEUTRINOS

For the relevant physical conditions, neutrinos can have a significant impact on the Kelvin-Helmholtz instability (i.e. the neutrino damping rate is larger than the growth rate) in an intermediate range of wavelengths between $\lambda_1$ and $\lambda_2$. $\lambda_1$ can be obtained from $k\Delta v/2 = \Gamma$ (drag regime), while $\lambda_2$ follows from $k\Delta v/2 = k^2\nu$ (viscous regime). Neutrinos are therefore important for wavenumbers

$$ \frac{\Delta v}{2\nu} < k < \frac{2\Gamma}{\Delta v}. \qquad (A1) $$

We can then write the condition for this range to exist and therefore for neutrinos to have an impact on the Kelvin-Helmholtz instability as

$$ \frac{\lambda_2}{\lambda_1} = \frac{4\Gamma\nu}{\Delta\nu^2} > 1. \qquad (A2) $$

Using equations (2) and (10), we find that this condition can be marginally satisfied only for high temperatures and low densities

$$ \frac{\lambda_2}{\lambda_1} \approx 0.04 \left( \frac{T}{30 \text{ MeV}} \right)^8 \left( \frac{\rho}{10^{14} \text{ g cm}^{-3}} \right)^{-2}. \qquad (A3) $$

This paper has been typeset from a TeX/Latex file prepared by the author.