# Tayler-Spruit dynamo in binary neutron star merger remnants

Alexis Reboul-Salze, Paul Barrère, Kenta Kiuchi, Jérôme Guilet, Raphaël Raynaud, Sho Fujibayashi, Masaru Shibata

► To cite this version:

Alexis Reboul-Salze, Paul Barrère, Kenta Kiuchi, Jérôme Guilet, Raphaël Raynaud, et al. Tayler-Spruit dynamo in binary neutron star merger remnants. Astronomy & Astrophysics - A&A, 2025, 699, pp.A4. 10.1051/0004-6361/202453126. hal-04843892

HAL Id: hal-04843892

https://hal.science/hal-04843892v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Tayler-Spruit dynamo in binary neutron star merger remnants

Alexis Reboul-Salze<fnref n="1" />,*, Paul Barrère<fnref n="2" />,*, Kenta Kiuchi<fnref n="1" />,<fnref n="5" />,*, Jérôme Guilet<fnref n="3" />, Raphaël Raynaud<fnref n="4" />,*, Sho Fujibayashi<fnref n="6" />,<fnref n="7" />,<fnref n="1" />,†, and Masaru Shibata<fnref n="1" />,<fnref n="5" />,‡

Received 22 November 2024 / Accepted 9 April 2025

## ABSTRACT

**Context.** In binary neutron star mergers, the remnant can be stabilized by differential rotation before it collapses into a black hole. Therefore, the angular momentum transport mechanisms are crucial for predicting the lifetime of the hypermassive neutron star. One such mechanism is the Tayler-Spruit dynamo, and recent simulations have shown that it could grow in proto-neutron stars that formed during supernova explosions.

**Aims.** We aim to investigate whether hypermassive neutron stars with high neutrino viscosity could be unstable to the Tayler-Spruit dynamo and study how magnetic fields would evolve in this context.

**Methods.** Using a one-zone model based on the result of a 3D GRMHD simulation, we investigate the time evolution of the magnetic fields generated by the Tayler-Spruit dynamo. In addition, we analyze the dynamics of the 3D GRMHD simulation to determine whether the dynamo is present.

**Results.** Our one-zone model predicts that the Tayler-Spruit dynamo can increase the toroidal magnetic field to $\ge 10^{17}$ G and the dipole field to amplitudes $\ge 10^{16}$ G. The dynamo's growth timescale depends on the initial large-scale magnetic field right after the merger. In the case of a long-lived hypermassive neutron star, an initial magnetic field of $\ge 10^{12}$ G would be enough for the magnetic field to be amplified in a few seconds. However, we show that the resolution of the current GRMHD simulations is insufficient to resolve the Tayler-Spruit dynamo due to high numerical dissipation at small scales.

**Conclusions.** We find that the Tayler-Spruit dynamo could occur in hypermassive neutron stars and shorten their lifetime, which would have consequences on multi-messenger observations.

**Key words.** dynamo – gravitational waves – magnetohydrodynamics (MHD) – binaries: general – gamma-ray burst: general – stars: neutron

## 1. Introduction

The association of GW170817 and GRB 170817A confirmed that binary neutron star (BNS) mergers are the progenitors of some short gamma-ray bursts (Abbott et al. 2017). However, the nature of GW170817’s remnant remains uncertain, and various physical parameters can impact the lifetime of the merger remnant before its collapse into a black hole. Depending on the hot and dense matter equation of state and the masses of the neutron stars, different fates await the remnant after the merger: if the remnant mass *M* remains lower than the maximum mass of a nonrotating neutron star *M*<sub>TOV</sub>, it stays indefinitely stable. Otherwise, the remnant can promptly collapse into a black hole, or it can form a hypermassive neutron star (HMNS) (Baumgarte et al. 2000), which is a neutron star that is stabilized by differential rotation while having a mass higher than the maximum mass of a uniformly rotating cold neutron star. An HMNS can either (i) collapse into a black hole when differential rotation is transported (Duez et al. 2006; Shibata et al. 2006), or (ii) be stabilized by solid body rotation, increasing *M*<sub>TOV</sub> by around ~20% (Cook et al. 1994a; Lasota et al. 1996; Breu & Rezzolla 2016; Musolino et al. 2024). In the latter case, the remnant is called

a *supramassive* neutron star (Cook et al. 1994b) that collapses after spinning down by magnetic dipole radiation on a timescale ranging from minutes to hours, or even longer. In the following, we call the long-lived remnant case when the remnant survives for *O*(1) s and it often corresponds to the *supramassive* neutron star case.

After the merger of two neutron stars, the nature of the remaining object can significantly affect the electromagnetic counterparts due to various factors. In the case of an HMNS, these include powerful neutrino luminosity, strong magnetic fields ($\ge 10^{15}$ G), and rotational kinetic energy reservoirs of order $10^{53}$ erg. Several studies have examined the impact of these factors on post-merger events (e.g., Metzger et al. 2008; Bucciantini et al. 2012; Gao et al. 2013; Metzger & Piro 2014; Gompertz et al. 2014; Sarin et al. 2022). Magnetic fields, along with fast rotation, play a critical role in neutron star mergers as they can extract rotational kinetic energy and drive powerful relativistic outflows (e.g., Shibata et al. 2021; Combi & Siegel 2023; Kiuchi et al. 2024). Therefore, understanding the evolution of the magnetic field and the rotation profile of the HMNS is crucial for interpreting future multi-messenger observations.

During the first few milliseconds after the merger, the magnetic field of the HMNS is amplified by the Kelvin-Helmholtz instability due to the presence of the shear

<footnote n="1">¹ Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam, Germany</footnote>

<footnote n="2">² Observatoire de Genève, Université de Genève, 51 Ch. Pegasi, 1290 Versoix, Switzerland</footnote>

<footnote n="3">³ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴ Université Paris Cité, Université Paris-Saclay, CNRS, CEA, AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="5">⁵ Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="6">⁶ Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 980-8578, Japan</footnote>

<footnote n="7">⁷ Astronomical Institute, Graduate School of Science, Tohoku University, Sendai 980-8578, Japan</footnote>

<footnote>* Corresponding author: alexis.reboul-salze@aei.mpg.de</footnote>
layer at the contact interface between the two neutron stars (Rasio & Shapiro 1999; Price & Rosswog 2006; Kiuchi et al. 2014, 2015; Aguilera-Miret et al. 2024). This generates a strong, small-scale magnetic field. However, another mechanism is necessary to get a large-scale magnetic field capable of driving an outflow. Both the magnetorotational instability (MRI) and the Tayler-Spruit (TS) dynamo could possibly fulfill this role, but act in different regions of the HMNS. Indeed, the remnant from the BNS merger has a rotation profile with an angular velocity that increases with radius until 10 km and then decreases afterwards (Shibata & Taniguchi 2006). Therefore, this rotation profile is unstable to the MRI at radii larger than 10 km. This instability is studied in the case of long-lived remnants, for which it is shown that it could drive a luminous relativistic outflow (Combi & Siegel 2023; Kiuchi et al. 2024). By contrast, the differential rotation is such that the flow is stable to the MRI at radii smaller than 10 km, and core regions could instead be unstable to a Tayler-instability driven dynamo (Tayler 1973; Spruit 2002).

The Tayler instability is a type of magnetic instability in flows of electrically conducting, stably stratified fluids where a strong azimuthal magnetic field can become unstable, causing disruptions due to the Lorentz forces that arise from the interaction of currents along the axis of symmetry and magnetic fields (Tayler 1973; Pitts & Tayler 1985). The first model of a dynamo driven by the Tayler instability was proposed by Spruit (2002) for stably stratified, differentially rotating regions. Fuller et al. (2019) have proposed an alternative description to address criticisms of Spruit's model (see Denissenkov & Pinsonneault 2007; Zahn et al. 2007). For a long time, numerical simulations were unable to capture the TS dynamo, but recent numerical simulations have finally provided numerical evidence that this dynamo exists, both in the presence of negative shear (Petitdemange et al. 2023, 2024; Daniel et al. 2023) and positive shear (Barrère et al. 2023, 2025). The former set-up is more relevant to model stellar radiative interiors, where the Tayler instability is proposed to explain the angular momentum evolution of stellar interiors (Fuller et al. 2019) and the magnetic dichotomy observed in intermediate-mass and massive stars (Rüdiger & Kitchatinov 2010; Szklarski & Arlt 2013; Bonanno et al. 2020; Jouve et al. 2020). On the other hand, a positive shear (with a faster-rotating surface) can be encountered in proto-neutron stars (PNSs). In this context, the TS dynamo plays a central role in a new scenario explaining magnetar formation from slow-rotating progenitors (Barrère et al. 2022).

The TS is also proposed as a mechanism for angular momentum transport potentially leading to their collapse into black holes (Margalit et al. 2022). This study estimates the angular momentum transport timescale due to the transport by the saturated magnetic fields from analytical scalings but neglects the duration of the dynamo growth. However, neutrino viscosity could make the TS dynamo stable, or it could require more time to reach saturation than the lifetime of the HMNS.

Therefore, this paper aims to study whether the conditions in BNS merger remnants could allow the TS dynamo to grow and impact its evolution. In particular, this study applies a similar one-zone model to that of Barrère et al. (2022), within the framework of BNS mergers. The paper is organized as follows: in Sect. 2, we describe the modelling of the HMNS and the equations for a one-zone model of the TS dynamo. The results of the one-zone models are presented in Sect. 3, followed by a comparison to the recent super-high resolution GRMHD simulation from Kiuchi et al. (2024) in Sect. 4. Finally, we discuss the validity of our assumptions in Sect. 5 and conclude in Sect. 6.

## 2. A one-zone model of the Tayler-Spruit dynamo in the hypermassive neutron star

### 2.1. Modeling of the hypermassive neutron star

For our HMNS model, we use the data from a 3D GRMHD simulation with a stiff equation of state, DD2 (Hempel & Schaffner-Bielich 2010), where the remnant lasts for more than 0.15 s (Kiuchi et al. 2024) (see Sect. 4 for more detail). To assess whether the TS dynamo is relevant for this HMNS model, we look at the rotation profile of the HMNS at different times (see top panel of Fig. 1). The remnant displays a cylindrical rotation profile with a large region in positive differential rotation between a cylindrical radius of $R = 5$ km and $R = 10$ km<fnref n="1" />, where the TS dynamo could grow. We stress that this rotation profile is typical of HMNSs in many different simulations (Shibata & Taniguchi 2006; Kastaun et al. 2016, 2017). To characterize the differential rotation, we use the shear rate

$$q \equiv \frac{d \ln \Omega}{d \ln R},$$

which means that a Keplerian shear corresponds to $q = -1.5$. In the following, we use the cylindrical radius $R = 7$ km as a typical radius for our thermodynamical quantities. At this radius, we have $\Omega = 5468$ s⁻<fnref n="1" /> and a shear rate of $q = 1.12$.

Since the TS dynamo requires a stable stratification, we test whether the core of the BNS remnant is convectively stable (Radice & Bernuzzi 2023). For that, we use the relativistic Ledoux criterion, given by

$$C_L = \frac{\partial \rho(1 + \epsilon)}{\partial r} - \frac{1}{c_s^2} \frac{\partial P}{\partial r}, \quad (1)$$

where $\rho$, $\epsilon$, $c_s$ and $P$ are the density, specific internal energy, sound speed and pressure, respectively. We find that the HMNS is indeed stably stratified, and we use the general relativistic formula for the Brunt-Väisälä (BV) frequency $N$ derived in Müller et al. (2013):

$$N^2 = -c^2 \frac{\alpha C_L}{\rho h \phi^4} \frac{\partial \alpha}{\partial r}, \quad (2)$$

where $\alpha$, $h$ and $\phi$ are the lapse function, specific enthalpy, and conformal factor for the spatial metric: $N$ is presented in the bottom panel of Fig. 1. Since it is rapidly rotating, the HMNS tends to be cylindrical. This means that the HMNS is more compact at the poles than at the equator, which leads to spherical radial gradients being stronger at the poles and a small increase in the Brunt-Väisälä frequency. This effect is more pronounced at larger radii as the frequency at $r = 7$ km stays relatively constant in theta. This formula comes from the hydrostatic equilibrium written as $G = -c^2 \nabla \ln \alpha$, where, in the Newtonian limit, $G$ is the gravitational acceleration. Due to the rapid rotation, we prefer to use this definition of hydrostatic equilibrium rather than the one based on the pressure gradient for our background. With this formula, we find typical values of $N = 4970$ s⁻<fnref n="1" /> at a radius of $R_{TI} = 7$ km for the Brunt-Väisälä frequency. The density is equal to $\rho = 3.7 \times 10^{14}$ g cm⁻³. Note that both the shear rate and the Brunt-Väisälä frequency defined here are not coordinate invariant.

In this context, we cannot assume that the rotational angular velocity $\Omega$ is much lower than the Brunt-Väisälä frequency. We

<footnote n="1">¹ In the following, the cylindrical radius is denoted by $R$ and the spherical one by $r$.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Top: Azimuthal averaged rotation profile of the remnant HMNS from the 3D GRMHD simulation as a function of the cylindrical radius, at different times t and vertical position z. Bottom: Averaged Brunt-Väisälä frequency of the remnant in the TS-unstable region as a function of the cylindrical radius, at different times t and different θ angles.</figcaption></figure>

also assume that the Alfvén frequency remains small compared to the rotation and Brunt-Väisälä frequencies,

$$ \omega_A \equiv \frac{B}{\sqrt{4\pi\rho r^2}} \ll \Omega \text{ and } \omega_A \ll N \Rightarrow B \ll 2 \times 10^{17} \text{ G}, \quad (3) $$

for the one-zone model. Note that this hypothesis may no longer be valid when we analyze the core dynamics in the GRMHD simulation, as the toroidal magnetic field can be winded to higher values. When these assumptions hold, the growth rate of the Tayler instability is $\sigma_{TI} = \omega_A^2/\Omega$. In the case of slow rotation $\Omega \ll \omega_A$, the growth rate is $\sigma_{TI} = \omega_A$.

## 2.2. Newtonian Tayler-Spruit picture

As in Barrère et al. (2022), we decompose the velocity and magnetic field components into axisymmetric components, denoted as $B_i$ or $v_i$ and non-axisymmetric components, denoted as $\delta B_i$ and $\delta v_i$. We justify the Newtonian approach in the Appendix.

The shear rate evolution due to angular momentum transport by Maxwell stresses is written as

$$ \partial_t q = -\gamma_{AM} q = \frac{T_{\text{merger}}^{\text{MAX}}}{4\pi\rho\Omega R_{\text{TI}}^2} = -\frac{B_R B_\phi}{4\pi\rho\Omega R_{\text{TI}}^2}, \quad (4) $$

where $T_{ij}^{\text{MAX}}$ is the Maxwell stress tensor. For the evolution of the angular velocity, we consider that the initial value decreases due to the magnetic torque. We also express it in terms of the Maxwell stress tensor, which gives

$$ \partial_t \Omega = \frac{R_{\text{TI}}^3 T_{\text{merger}}^{\text{MAX}}}{I} = -\frac{R_{\text{TI}}^3 B_R B_\phi}{I}, \quad (5) $$

where $I = 1.667 \times 10^{45} \text{ g cm}^2$ is the moment of inertia of the inner zone below 10 km. This means that in the case of a strong magnetic field, the rotation rate could decrease fast and be in the regime $\Omega \ll \omega_A$. In the following, we keep the growth rate of the Tayler instability $\sigma_{TI}$ that is equal to

$$ \sigma_{\text{TI}} = \begin{cases} \omega_A & \text{if } \Omega \le \omega_A, \\ \omega_A^2/\Omega & \text{if } \omega_A \ge \Omega, \end{cases} \quad (6) $$

in the evolution equations (Eqs. (7), (12) and (14)) of the magnetic field. For the dissipation of the perturbed magnetic field, we consider the same damping rate $\gamma_{cas} \sim \frac{\delta v_A}{r}$ as Fuller et al. (2019), where $\delta v_A \equiv \frac{\delta B_\phi}{\sqrt{4\pi\rho}}$ is the Alfvén velocity of the non-axisymmetric $\delta B_\phi$. It comes from the cascade rate of the perturbed magnetic field towards smaller scales. The evolution equation of $\delta B_\phi$ is therefore

$$ \partial_t \delta B_\phi = (\mathcal{B}_{\text{TI}}\sigma_{\text{TI}} - \gamma_{\text{cas}}) \delta B_\phi = \mathcal{B}_{\text{TI}}\sigma_{\text{TI}}\delta B_\phi - \frac{\delta v_A}{R_{\text{TI}}}\delta B_\phi, \quad (7) $$

where $\mathcal{B}_{\text{TI}}$ is a boolean that describes whether the HMNS is unstable to the Tayler instability (see Eq. (24)). In the case of ideal MHD, the instability is stabilized when the magnetic tension due to the radial magnetic field, $\sim k_{TI}B_R\delta B$, where $k_{TI}$ is the wavenumber of the Tayler instability, overcomes the driving force of the Tayler instability, the magnetic tension due to the toroidal magnetic field $\sim m/rB_\phi\delta B$, with $m = 1$ being the main mode of the Tayler instability (Fuller et al. 2019; Skoutnev & Beloborodov 2024). We have therefore

$$ \mathcal{B}_{\text{TI}} = \begin{cases} 0 & \text{if } B_R/B_\phi > \omega_A/N, \\ 1 & \text{otherwise.} \end{cases} \quad (8) $$

One other difference between the HMNS model and the PNS model of Barrère et al. (2022) is the geometry: the rotation profile in an HMNS is cylindrical while it is spherical for a late stage of the PNS with fallback accretion. In our study, we, therefore, consider the cylindrical radial field $B_R$ instead of the spherical radial field. The derivation of the equations is almost the same as it is just adapted to cylindrical coordinates. With a cylindrical differential rotation, the growth of the toroidal field $B_\phi$ is due to the winding of the cylindrical radial field $B_R$, described by $\partial_t B_\phi|_{\text{winding}} = q\Omega B_R$. The growth rate of the toroidal field is therefore $\sigma_{\text{shear}} = q\Omega B_R/B_\phi$. As in Barrère et al. (2022), we estimate the dissipation rates of $B_R$ and $B_\phi$ by the non-linear magnetic energy dissipation rate $\dot{E}_{\text{mag}}$ following

$$ \gamma_{\text{diss}} \equiv \frac{\dot{E}_{\text{mag}}}{B_{\phi}^2}, \quad (9) $$

where the radial magnetic field is neglected compared to the dominant azimuthal field such that we can estimate $\dot{E}_{\text{mag}} \sim B_\phi\partial_t B_\phi$. The change in axisymmetric magnetic energy is due to the growth of the Tayler instability, and thus, we have

$$ \dot{E}_{\text{mag}} \sim \sigma_{\text{TI}} |\delta B_{\phi}|^2. \quad (10) $$
It finally gives a dissipation rate

$$ \gamma_{\text{diss}} \sim \sigma_{\text{TI}} \left( \frac{\delta B_{\phi}}{B_{\phi}} \right)^2 . \qquad (11) $$

We take the simplification in order of magnitude of $\delta B_{\perp} \sim \delta B_{\phi}$, as $\delta B_{\perp}$ is not well defined for cylindrical coordinates. The evolution equation of $B_{\phi}$ is therefore

$$ \partial_t B_{\phi} = (\sigma_{\text{shear}} - \gamma_{\text{diss}}) B_{\phi} = q \Omega B_{\text{R}} - \sigma_{\text{TI}} \frac{\delta B_{\phi}^2}{B_{\phi}}. \qquad (12) $$

In cylindrical coordinates, the growth rate $\sigma_{\text{NL}}$ of $B_{\text{R}}$ comes from non-linear effects of the perturbed magnetic field and velocity. This corresponds to the toroidal component of the electromagnetic force $\mathcal{E}_{\phi}$ that is given by

$$ \frac{\partial B_{\text{R}}}{\partial t} = -\frac{\partial \mathcal{E}_{\phi}}{\partial z} \sim \frac{\mathcal{E}_{\phi}}{L_z}, \qquad (13) $$

where $L_z$ is the vertical wavelength of the Tayler instability. By definition, in cylindrical coordinates, the toroidal component of the electromotive force is given by $\mathcal{E}_{\phi} = \delta v_{\text{R}} \delta B_z - \delta v_z \delta B_{\text{R}}$. As in Barrère et al. (2022), we assume $\mathcal{E}_{\phi} \sim \delta v_{\text{R}} \delta B_z \sim \delta v_{\text{R}} \delta B_{\perp} \sim \delta v_{\text{R}} \delta B_{\phi}$. By assuming a magnetostrophic balance between the Lorentz and Coriolis forces, which gives $\delta v_{\perp} \Omega \sim \delta v_{\text{A}} \omega_{\text{A}}$, we retrieve the same formula as previous studies in spherical geometry if we assume $L_z \sim R_{\text{TI}}$, which we do for simplicity. For consistency, we also assume that the non-linear growth rate of the axisymmetric radial field is zero when the HMNS is stable to the Tayler instability. The evolution equation of $B_{\text{R}}$ is therefore

$$ \partial_t B_{\text{R}} = (\sigma_{\text{NL}} - \gamma_{\text{diss}}) B_{\text{R}} = \mathcal{B}_{\text{TI}} \frac{\omega_{\text{A}}^2}{N\Omega} \frac{\delta B_{\phi}^2}{\sqrt{4\pi\rho R_{\text{TI}}^2}} - \sigma_{\text{TI}} \left(\frac{\delta B_{\phi}}{B_{\phi}}\right)^2 B_{\text{R}}. \qquad (14) $$

Overall, the magnetic field evolution is governed by equations similar to the equations of Barrère et al. (2022), with some changes to take into account the growth rate and the stabilization of the instability.

## 2.3. Impact of diffusive processes on the Tayler-Spruit dynamo

In an HMNS, the impact of neutrinos on the momentum equation can be modeled as a strong viscosity for scales larger than the neutrino mean free path. In the HMNS, the kinematic shear viscosity can be estimated by the formula (Keil et al. 1996; Guilet et al. 2015)

$$ \nu = 1.2 \times 10^8 \left( \frac{\rho}{10^{14} \text{ g cm}^{-3}} \right)^{-2} \left( \frac{T}{10 \text{ MeV}} \right)^2 \text{cm}^2 \text{s}^{-1}, \qquad (15) $$

which gives $\nu = 7.4 \times 10^7 \text{ cm}^2 \text{s}^{-1}$ at $R = 7$ km with $\rho = 3.7 \times 10^{14} \text{ g cm}^{-3}$ and $T = 29$ MeV. As there are temperature and density gradients in the 5 to 10 km region where the TS dynamo could be active, the viscosity would lie in the range $\nu = 1.3 \times 10^7 - 4.3 \times 10^8 \text{ cm}^2 \text{s}^{-1}$. Note that this formula also applies to the region outside 10 km, which is also stably stratified but where the MRI is most likely the dominant instability (see Section 5.2).

This viscosity would lead to a magnetic Prandtl number of

$$ Pm \equiv \frac{\nu}{\eta} = 5.8 \times 10^{12}, \qquad (16) $$

where $\eta$ is assumed to be the resistivity due to electron-proton scattering (Thompson & Duncan 1993). The thermal diffusion coefficient $\kappa$ due to neutrinos in core-collapse supernovae is higher than that by the viscous effect and the thermal Prandtl number $Pr$ is found to be

$$ Pr \equiv \frac{\nu}{\kappa} = 10^{-3} - 10^{-2}, \qquad (17) $$

(Raynaud et al. 2020; Reboul-Salze et al. 2022), and we assume a similar ratio here.

We checked that the assumption of a magnetostrophic balance in the formalism of Fuller et al. (2019) is verified and not changed to a viscous-Lorentz force balance with the neutrino viscosity (see Appendix B).

The HMNS could be stabilized due to the fast dissipation of the perturbed velocity. We therefore estimate the critical magnetic field $B_{\phi,crit}$ to be unstable by checking the allowed wavenumber $k_{\text{TI}}$ of the instability. The minimum of this wavenumber has been estimated by Spruit (2002) as

$$ k_{\text{TI}} > \frac{N}{R_{\text{TI}}\omega_{\text{A}}}. \qquad (18) $$

We also estimate that the effect of viscosity is not dominant if the growth timescale of the Tayler instability is faster than the viscous timescale, which limits the maximum wavenumber to be smaller than

$$ k_{\text{TI}}^2 < \frac{\sigma_{\text{TI}}}{\nu} = \frac{\omega_{\text{A}}^2}{\Omega\nu}. \qquad (19) $$

In order to satisfy both limits, the magnetic field must be stronger than the critical magnetic field

$$ B_{\phi,crit} = (4\pi\rho R_{\text{TI}}N)^{1/2} (\nu\Omega)^{1/4} \sim 3.2 \times 10^{15} \text{ G}, \qquad (20) $$

This value is much higher than that found in Barrère et al. (2022) because we consider the neutrino viscosity, which is much higher than the fluid viscosity or resistivity. Due to the viscosity variation in the region from $R = 5-10$ km, the critical field estimate lies in the range $B_{\phi,crit} = 1.1 \times 10^{15} - 4.8 \times 10^{15} \text{ G}$. Note that in the case of the slow-rotating limit, the growth rate changes and the critical magnetic field would be given by

$$ B_{\phi,crit} = (4\pi\rho R_{\text{TI}})^{1/2} \left( \frac{\nu N^2}{R_{\text{TI}}^2} \right)^{1/3}. \qquad (21) $$

For simplification, the wavelength of the Tayler instability $k_{\text{TI}}$ is assumed to be equal to $k_{\text{TI,crit}}$ for all the evolution of the one-zone model. This critical magnetic field strength gives the corresponding critical wavenumber

$$ k_{\text{TI,crit}} = \left( \frac{\sigma_{\text{TI,crit}}}{\nu} \right)^{1/2}. \qquad (22) $$

This high critical strength for our model means that the HMNS is stable to the Tayler instability during the time needed for the toroidal field to be winded up by the shear. We note that this critical magnetic field strength is a conservative estimate as the recent study of Skoutnev & Beloborodov (2024) shows that, in the case of a large cylindrical gradient of the toroidal field (corresponding to the geometric criterion of Goossens & Tayler 1980), the criterion depending on diffusivities is

$$ \frac{\omega_{\text{A}}}{2\Omega} \gg \left(\frac{N}{2\Omega}\right)^{1/2} \left(\frac{\eta k_{\theta}^2}{2\Omega}\right)^{1/4} \min\left(1, \left(\frac{\eta}{\kappa}\right)^{1/4}\right) \approx 2.15 \times 10^{-9}, \qquad (23) $$
where $k_\theta$ is the wavelength along the vector $e_\theta$ in spherical coordinates. This gives a lower limit for the magnetic field of $B_{\phi,crit} > 1.1 \times 10^9$ G. This limit is easily met through the winding of the magnetic field, and hence, we keep the viscous criterion to be conservative.

In order to take into account the impact of the viscosity in the evolution equations, we modify the boolean $B_{TI}$ (Eq. (8)) in the evolution equation of $\delta B_\phi$ (Eq. (7)) to include the cases where the HMNS is stable to the Tayler instability. It gives now

$$B_{TI} = \begin{cases} 0 & \text{if } B_{\phi} < B_{\phi,crit} \text{ or } B_{R}/B_{\phi} > \omega_{A}/N, \\ 1 & \text{otherwise.} \end{cases} \qquad (24)$$

The system can therefore be summarized by the following equations

$$\partial_t B_\phi = (\sigma_{\text{shear}} - \gamma_{\text{diss}}) B_\phi = q\Omega B_R - \sigma_{\text{TI}} \frac{\delta B_\phi^2}{B_\phi}, \qquad (25)$$

$$\partial_t \delta B_\phi = (B_{TI} \sigma_{TI} - \gamma_{\text{cas}}) \delta B_\phi = B_{TI} \sigma_{TI} \delta B_\phi - \frac{\delta v_A}{R_{TI}} \delta B_\phi, \qquad (26)$$

$$\partial_t B_R = (\sigma_{\text{NL}} - \gamma_{\text{diss}}) B_R = B_{\text{TI}} \frac{\omega_A^2}{N\Omega} \frac{\delta B_\phi^2}{\sqrt{4\pi\rho R_{\text{TI}}^2}} - \sigma_{\text{TI}} \left(\frac{\delta B_\phi}{B_\phi}\right)^2 B_R. \qquad (27)$$

# 3. Results of the one-zone model

## 3.1. A fiducial case

We first consider a fiducial case, in which the initial large-scale radial magnetic field is assumed to be $B_{R,0} = 10^{13}$ G. The perturbed magnetic field also needs to be initialized, and we consider that the first perturbations have a similar order of magnitude, $\delta B_{\phi,0} = B_{R,0} = 10^{13}$ G. The azimuthal magnetic field is initialized to $B_{\phi,0} = 1$ G but this choice has little influence, as the winding is very efficient at amplifying this component of the magnetic field. Figure 2 shows the time evolution of the magnetic field, which divides into three phases as expected:

(i) First, the shear amplifies the toroidal magnetic field. This winding phase lasts for $\approx 0.15$ s.

(ii) When the toroidal magnetic field is stronger than the critical magnetic field $B_{\phi,crit}$, the non-axisymmetric magnetic field $\delta B_\phi$ grows due to the Tayler instability. The duration of this phase is $\approx 0.25$ s.

(iii) When $\delta B_\phi$ becomes strong enough for the non-linear induction to be relevant, the dynamo sets in and increases all the magnetic field components for the next 0.08 s. The dynamo stops when the angular momentum is transported and the HMNS core settles in a rigid rotation with a low angular velocity (see bottom panel in Fig. 2).

The saturation is reached after several hundred milliseconds. Note that we neglect the contraction of the HMNS due to neutrino cooling on this timescale. This would actually increase the shear rate as the rotation maximum shifts inside, and, for cases with high mass or soft EOS, the core of the HMNS may also collapse before the saturation.

At saturation, the toroidal, radial, and perturbed magnetic fields are equal to

$$B_{\phi}^{\text{sat}} = 1.6 \times 10^{17} \text{ G}, \qquad (28)$$

$$B_{R}^{\text{sat}} = 5.9 \times 10^{16} \text{ G}, \qquad (29)$$

$$\delta B_{\phi}^{\text{sat}} = 1.0 \times 10^{17} \text{ G}. \qquad (30)$$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Top: time evolution of the magnetic field strength for a typical HMNS with thermodynamic quantities taken at $t - t_{merger} = 0.03$ s at a radius of $R_{TI} = 7$ km and with an initial magnetic field strength of $B_{R,0} = 10^{13}$ G. Bottom: time evolution of the rotational angular velocity $\Omega$ and the shear rate $q$.</figcaption></figure>

With the saturated toroidal magnetic field, we have $\frac{\omega_{A,sat}}{\Omega} = 0.58$, meaning that the initial assumption of $\omega_A \ll \Omega$ is verified for most of the dynamo phase until the last rapid growth where the rotational angular velocity decreases to low values due to the magnetic torque. We compare the saturation values of the axisymmetric field $B_R$ and $B_\phi$ obtained by the model to the estimated saturated magnetic field from the following formulas (Fuller et al. 2019)

$$B_{\phi}^{\text{F,sat}} = \sqrt{4\pi\rho R_{\text{TI}}^2 \Omega} \left(\frac{q\Omega}{N}\right)^{1/3} \approx 2.8 \times 10^{17} \text{ G}, \qquad (31)$$

$$B_{r}^{\text{F,sat}} = \sqrt{4\pi\rho R_{\text{TI}}^2 \Omega} \left(\frac{q^2\Omega^5}{N^5}\right)^{1/3} \approx 3.3 \times 10^{17} \text{ G}, \qquad (32)$$

$$\delta B_{\phi}^{\text{F,sat}} = \sqrt{4\pi\rho R_{\text{TI}}^2 \Omega} \left(\frac{q\Omega}{N}\right)^{2/3} \approx 3.0 \times 10^{17} \text{ G}. \qquad (33)$$

These values differ from the results of our models because some of the model's assumptions are not valid before reaching the predicted saturated values. Indeed, both Fuller's and Spruit's predictions give $\omega_A > \Omega$ and $\omega_A > N$ with our values. These values are not reached, as the differential rotation is quenched before. Nonetheless, the saturated values still agree within a factor of 3–5.
<figure><img src="image_4.png" /></figure>

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

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Time evolution of the magnetic field for different initial magnetic dipole strength $B_{R,0}$ with $B_{R,0} = 10^{14}$ G (top), $B_{R,0} = 10^{12}$ G (middle) and $B_{R,0} = 10^{10}$ G (bottom).</figcaption></figure>

It is interesting to note that the alternative model of Spruit (2002) predicts similar magnetic field strength at saturation:

$$B_{\phi}^{\text{S,sat}} = \sqrt{4\pi\rho R_{\text{TI}}^2} q \left(\frac{\Omega^2}{N}\right) \approx 3.2 \times 10^{17} \text{ G} \quad (34)$$

<table><caption>**Table 1.** Initial magnetic field strength $B_{R,0}$ and time to reach saturation in the one-zone model.</caption><thead><tr><th><i>B</i><sub>R,0</sub> [G]</th><th>Δ<i>t</i><sub>sat</sub> [s]</th></tr></thead><tbody><tr><td>10<sup>14</sup></td><td>≈0.11</td></tr><tr><td>10<sup>13</sup></td><td>≈0.48</td></tr><tr><td>10<sup>12</sup></td><td>≈2.4</td></tr><tr><td>10<sup>10</sup></td><td>≈86</td></tr></tbody></table>

$$B_r^{\text{S,sat}} = \sqrt{4\pi\rho R_{\text{TI}}^2} q^2 \left(\frac{\Omega^4}{N^3}\right) \approx 3.9 \times 10^{17} \text{ G.} \quad (35)$$

This can be interpreted by the fact that the rotation and Brunt-Väisälä frequencies are of the same order of magnitude. Overall, with an initial magnetic field of $10^{13}$ G, the TS-dynamo could generate an intense dipole and an even stronger toroidal magnetic field when the long-lived HMNS survives for longer than 0.5 s after the BNS mergers.

## 3.2. Dependence on the initial dipole magnetic field strength

For BNS mergers, the duration of the dynamo is quite essential as it may determine the time at which the remnant collapses into a black hole. This duration is expected to depend on the initial intensity of the large-scale poloidal magnetic field, which is uncertain as discussed above. We therefore vary the initial magnetic field strength to study its influence on the dynamo. Figure 3 shows the time evolution for different initial magnetic fields $B_{R,0}$. While the maximum magnetic field strength is similar for all cases, the time to reach saturation $\Delta t_{\text{sat}}$ decreases drastically with the increase of the initial magnetic field strength (see Table 1).

To realistically estimate the initial dipole field strength, the magnetic field in the core of the HMNS, amplified by the Kelvin-Helmholtz (KH) instability at the onset of the merger, should be explored with dissipation (Rasio & Shapiro 1999; Kiuchi et al. 2014, 2018, 2024). Due to this instability, the large-scale initial magnetic field could be amplified and potentially be stronger than $B_{R,0} = 10^{13}$ G. Therefore, the TS dynamo could saturate in less than 0.5 s.

To complete this picture, we can estimate the duration of each phase of the dynamo process. The first phase is the winding phase, during which the toroidal field is amplified. It lasts until either the toroidal field becomes strong enough for the Tayler instability to grow as fast as $B_{\phi}$ (Barrère et al. 2022), or that it becomes strong enough to reach the critical magnetic field for the instability to start with the high viscosity. The formula for the winding time is therefore

$$t_{\text{winding}} = \frac{\max(\omega_{\text{A,TI}}, \omega_{\text{A,crit}})}{q\Omega\omega_{\text{R,0}}} = \max\left((q^2\omega_{\text{R,0}}^2\Omega)^{-1/3}, \frac{\omega_{\text{A,crit}}}{q\Omega\omega_{\text{R,0}}}\right), \quad (36)$$

where $\omega_{\text{A,crit}}$ and $\omega_{\text{R,0}}$ are the Alfvén frequency of the critical magnetic field due to viscosity and the initial Alfvén frequency, respectively. $\omega_{\text{A,TI}}$ is the Alfvén frequency when the Tayler instability becomes dynamically relevant, that is when $\sigma_{\text{shear}} = \sigma_{\text{TI}}$, which is given by

$$\omega_{\text{A,TI}} = (q\Omega^2\omega_{\text{R,0}})^{1/3}. \quad (37)$$

At the end of the winding phase, the toroidal magnetic field is equal to

$$B_{\phi}(t_{\text{winding}}) = q\Omega t_{\text{winding}} B_{\text{R,0}}. \quad (38)$$
<figure><img src="image_7.png" /><figcaption>Fig. 4. Timescale of the different phases post-merger of the TS dynamo as a function of the initial radial magnetic field strength. The points are measured from the one-zone model, while the lines represent analytical estimates. The vertical dashed line indicates the volume-averaged poloidal magnetic field in Aguilera-Miret et al. (2024).</figcaption></figure>

For the Tayler instability phase, we adapt Eq. (51) of Barrère et al. (2022) considering that $B_\phi(t_{winding})$ is given by Eq. (38). Introducing the corresponding Alfvén frequency $\omega_{A,winding}$, we find

$$t_{TI} \sim \frac{\Omega}{(\omega_{A,winding})^2} \ln \left( \frac{\omega_{A,winding}^2 r}{\Omega \delta v_{A,0}} \right), \quad (39)$$

where $\delta v_{A,0} = \delta v_A(t=0)$.

Finally, we use the same formula as Barrère et al. (2022) for the dynamo phase

$$t_{dyn} = \left( \frac{B_{\phi}}{\partial_{t}^{2} B_{\phi}} \right)^{1/2} = \frac{\Omega}{\omega_{A,dyn}^{2}} \left( \frac{N}{q \omega_{A,dyn}} \right)^{1/2}, \quad (40)$$

where $\omega_{A,dyn} = (qN\Omega^4\omega_{R,0}^2)^{1/7}$ is the Alfvén frequency for which the winding rate is equal to the non-linear growth rate $\sigma_{NL}$ of the radial magnetic field $B_r$.

Figure 4 compares these timescale estimates to the ones measured in the one-zone model. The analytical predictions are overall in good agreement with only a slight overestimation compared to the model results. We also give a simple fit of the total dynamo time

$$t_{total} \approx 65.7 \left( \frac{B_{R,0}}{10^{10}} \right)^{-0.71} \text{ s.} \quad (41)$$

This shows that the TS dynamo can grow in a few seconds as long as the initial magnetic field strength is higher than $10^{12}$ G. This figure shows that the shear time dominates the total time for initial magnetic field strengths lower than $10^{12}$ G. In this case, a higher rotational angular velocity or shear rate would linearly decrease the total dynamo time. For higher initial magnetic field strength, the Tayler instability phase dominates the total time, and therefore, the Tayler-instability growth rate becomes a key parameter. In the first few milliseconds after the mergers, the Kelvin-Helmholtz instability amplifies the magnetic field, which places a lower bound on the radial component of the initial large-scale magnetic field for our one-zone model. However, the strength of this lower bound remains unknown. As

an initial guess, we take the volume-averaged poloidal magnetic field $B_{R,0} \sim 5 \times 10^{14}$ G found in the simulations of Kelvin-Helmholtz instability with a realistic initial field (Fig. 4 of Aguilera-Miret et al. 2024). An initial magnetic field $B_{R,0} \sim 5 \times 10^{14}$ G means that it is plausible that the TS dynamo saturates after $\le 50$ ms.

## 4. Comparison with a 3D GRMHD simulation

In the previous section, we showed that the TS dynamo is able to grow and reach saturation on a timescale of a few tens of milliseconds for an initial dipole higher than $10^{15}$ G, which is the typical dipole value used in GRMHD simulations of BNS mergers (e.g., Giacomazzo & Perna 2013; Kiuchi et al. 2014, 2024; Mösta et al. 2020). As a consequence, the magnetic fields would be expected to grow from the TS dynamo in such GRMHD simulations. In this section, we analyze the data of a 3D ideal GRMHD simulation in which the massive neutron star remnant is long-lived (Kiuchi et al. 2024). Since the remnant is differentially rotating during the duration of the simulation (see Figure 1), we call it an HMNS in the following. Note that it may become a supramassive neutron star once in solid body rotation or collapse to a black hole in the case of a strong spin-down.

Since we look for the TS dynamo, we focus our analysis on the core of the BNS merger remnant at radii lower than 10 km. In this region, the resolution is $\Delta x = 12.5$ m in the range [0, 4.5] km and $\Delta x = 25$ m in the range [4.5, 9] km, which is roughly the region where the positive differential rotation is found. In this region, the resolution becomes $\Delta x = 100$ m after $t - t_{merger} = 0.050$ s. The numerical techniques and the study of the other regions of these simulations are detailed in Kiuchi et al. (2024). Note that the outer surface of the merger remnant core for which $q$ is negative is subject to the $\alpha\Omega$ dynamo as reported in Kiuchi et al. (2024).

### 4.1. Dynamics of the GRMHD simulation

The 3D neutrino radiation GRMHD simulation for a BNS merger of a symmetric binary case, with masses of $1.35-1.35 M_\odot$ and a stiff equation of state, DD2 (Hempel & Schaffner-Bielich 2010), shows that the merger results in a HMNS that survives for timescales $>0.15$ s. The rotation profile in the equatorial plane of the HMNS shows a differential rotation with $q > 0$ until 10 km (see Fig. 1). This leads to the winding of the poloidal magnetic field in the core of the remnant. Figure 5 compares the time evolution of the magnetic field strength in the 3D GRMHD simulation and the one-zone model. The radial magnetic field in the one-zone model is amplified to its saturation values within 0.02 s, while in the GRMHD simulation, the radial magnetic field strength varies slightly and overall remains constant, staying close to its initial value around $10^{15}$ G. This shows that no TS dynamo processes occur in the simulation's HMNS core even though it would be expected according to the theoretical prediction. In addition, the non-axisymmetric vertical field $\delta B_z$ is only decreasing in this region, while it should be growing to a similar amplitude as the $\delta B_\phi$ of the one-zone model. In this simulation, the toroidal magnetic field quickly reaches $10^{17}$ G due to the winding around 7 km. This gives an Alfvén frequency of $\omega_A = 2100 s^{-1}$, which is smaller but of the same order of magnitude as $N \approx 5000 s^{-1} \sim \Omega \approx 5500 s^{-1}$. The local Alfvén frequency may even be closer to the rotation frequency, as the local field can be stronger than $10^{17}$ G.
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Time evolution of the axisymmetric magnetic field and $\delta B_z$ for the GRMHD simulation (solid lines) in the unstable region (see Section 2.1 and Figure 7) and for the one-zone model (dashed lines) that has been shifted temporally.</figcaption></figure>

To understand why there is no TS dynamo, we look at a snapshot of the simulation in a plane at $z = 2$ km and $t - t_{merger} = 0.03$ s (see Fig. 6). We see that the toroidal magnetic field is amplified to a maximum value of $4 \times 10^{17}$ G (top panel of Figure 6). Since this is larger than the critical magnetic field for a strong viscosity, the simulation should be unstable to the Tayler instability but the non-axisymmetric component of the radial magnetic field does not show any evidence for the Tayler instability. Indeed, its main mode is an $m = 2$ mode, and the geometry corresponds well to the velocity field in the equatorial plane (bottom panel of Fig. 6). The origin of this $m = 2$ mode is well-known and comes from the initial geometry and mass ratio of the binary. Figure 6 shows that the velocity field is not a purely $m = 2$ mode as the initial $m = 2$ mode can also degrade as a $m = 1$ mode, but it is different from the $m = 1$ mode from the Tayler instability. Indeed, the $m = 2$ and degraded $m = 1$ modes give the radial magnetic field exactly the same geometry as the radial velocity field and it has the ring shape of the toroidal field. In addition, according to simulations of the Tayler instability (Ji et al. 2023; Barrère et al. 2025), we would expect to see some opposing polarities of the radial field on a meridional slice, contrary to the geometry observed in Figure 7, which is the same as the toroidal field. Note that, in the case of an unequal mass binary neutron star merger, the $m = 1$ mode may be more easily triggered due to the asymmetry of the initial system.

Positive differential rotation can stabilize the Tayler instability, as shown by Kiuchi et al. (2011), but the toroidal magnetic field is decreasing vertically in the north hemisphere in the unstable region (see Fig. 7), which makes it unstable (see Eq. (27) of Kiuchi et al. 2011). From this criterion, the core of the remnant at $t - t_{merger} = 0.06$ s should be unstable to the Tayler instability with a fast-growing rate of $\sigma_{TI} \sim \omega_A$, because we are in the “slow” rotating regime.

To be unstable to the Tayler instability, the geometry of the magnetic field is also important, as the size of the magnetic field can limit the wavelength. To estimate what is the maximal wavelength available, we used the geometrical instability criterion for a wave from Goossens & Tayler (1980)

$$ \frac{B_{\phi}^{2}}{4\pi R_{\text{TI}}^{2} \sin \theta} (m^{2} - 2 \cos^{2} \theta - 2 \sin \theta \cos \theta \partial_{\theta} \log B_{\phi}) < 0. \quad (42) $$

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

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

<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Snapshot of the toroidal magnetic field (top), radial magnetic field (middle) and radial velocity (bottom) in the core of the HMNS, taken in the plane $z = 2$ km at $t - t_{merger} = 0.035$ s.</figcaption></figure>

Figure 7 shows the toroidal magnetic field with the gray contours showing where the Tayler instability is stable according to this criterion with $m = 1$. The magnetic field should indeed be unstable to the Tayler instability, but the unstable domain is relatively small, with a typical size of around 1 km. Contrary to the one-zone model, the hypothesis $\omega_A \ll \Omega$ is not verified in the unstable domain.

Lastly, another instability criterion is that $B_R/B_\phi < k_\phi/k_R$ which implies that the magnetic tension due to the radial component should not exceed the magnetic pressure due to the toroidal com-
ponent (Eq. (8)) to operate the Tayler instability. This criterion is usually given for an axisymmetric $B_R$ but it is unclear whether the magnetic tension from a non-axisymmetric $B_{R,m \neq 0}$ would stabilize the instability by overcoming the magnetic tension from the background toroidal field.

Looking at the middle panel of Fig. 6, this criterion may explain why the Tayler instability is temporarily stabilized as the radial length of $B_R$ is quite small. However, after some time, $B_\phi$ grows in strength and size with the winding (Figure 7). At the same time, the amplitude of the non-axisymmetric modes decreases and the instability criterion becomes verified.

In any case, these non-axisymmetric modes amplify the magnetic field to a lower strength than the ones in the one-zone model and are expected to dissipate less than a hundred milliseconds after the merger as the gravitational wave luminosity decreases (Radice & Bernuzzi 2023). These modes might cause a delay in the appearance of the TS dynamo, but this does not explain why it does not appear later in the simulation.

## 4.2. Numerical dissipation analysis

The absence of the Tayler instability in this 3D GRMHD simulation could be due to the wavelength of the Tayler instability being lower than the grid resolution or to the fact that the numerical resistivity and viscosity reduce the growth rate of the instability. Thus, we estimate the expected wavelength of the Tayler instability if it were growing in the simulation. In ideal MHD, the wavelength would be limited by the buoyancy and would give $\lambda_{TI} = \frac{\omega_A}{N}R \approx R$ in the simulation. This wavelength is well resolved so we need to look at the impact of diffusion on the wavelength. From the instability criterion, the wavelength $\lambda_{crit} = 2\pi/k_{TI,crit}$ below which the Tayler instability should occur is

$$ \lambda_{crit} = \frac{2\pi R_{TI} \omega_{A,crit}}{N} = \begin{cases} 2\pi R_{TI} \left(\frac{\Omega}{N}\right)^{1/2} \left(\frac{v}{R_{TI}^2 \Omega}\right)^{1/4} & \text{if } \omega_A \ll \Omega, \\ 2\pi R_{TI} \left(\frac{v}{R_{TI}^2 N}\right)^{1/3} & \text{if } \Omega \ll \omega_A. \end{cases} \quad (43) $$

As this wavelength depends on the numerical resistivity and viscosity inside the simulation, we estimate the numerical resistivity and viscosity. We use the fitting formula derived for Eulerian MHD codes in Rembiasz et al. (2017)

$$ v_{\star} = \mathcal{R}_{\nu}^{\Delta x} \mathcal{V} \mathcal{L} \left(\frac{\Delta x}{\mathcal{L}}\right)^{r}, \quad (44) $$

where $\mathcal{V}$ and $\mathcal{L}$ are the characteristic speed and length of the system, $\mathcal{R}_{\nu}^{\Delta x}$ and $r$ are fitting parameters. We use the fitting parameters $r = 4.95$ and $\mathcal{R}_{\nu}^{\Delta x} = 42$ for the HLLD solver with the fourth-order Runge Kutta method, as it is the same Riemann solver used in our GRMHD simulations. We use the fast magnetosonic speed for the characteristic velocity as found in the study $\mathcal{V} = \max(v_A, c_s) \approx c_s \approx 1.2 \times 10^{10} \text{ cm s}^{-1}$ for the characteristic speed. We take the characteristic length as a variable $\mathcal{L} = \lambda_{crit}$. Since the numerical resistivity depends on the resolution and on the scale we are considering, by combining Eqs. (43) and (44), we have to solve the equation

$$ \left\{ \begin{aligned} \lambda_{crit}^{1.9875} &= 2\pi R_{TI} \left(\frac{\Omega}{N}\right)^{1/2} \left(\frac{\mathcal{R}_{\nu}^{\Delta x} \mathcal{V} \Delta x^{4.95}}{R_{TI}^2 \Omega}\right)^{1/4}, & \text{if } \omega_A \ll \Omega, \\ \lambda_{crit}^{2.317} &= 2\pi \left(\frac{\mathcal{R}_{\nu}^{\Delta x} \mathcal{V} \Delta x^{4.95} R_{TI}}{N}\right)^{1/3}, & \text{if } \omega_A \gg \Omega, \end{aligned} \right. \quad (45) $$

<figure><img src="image_12.png" /><figcaption>Fig. 7. Snapshot of the axisymmetric toroidal magnetic field (top) and a slice of the radial magnetic field $B_R$ (bottom) in the core of the remnant neutron star at $t - t_{merger} = 0.060 \text{ s}$ and $y = 0 \text{ km}$. The gray contours correspond to the regions stable to the Tayler instability according to the geometric criterion (Eq. (42) with $m = 1$). Density contours for $\rho \in [10^{12}, 10^{13}, 10^{14}] \text{ g cm}^{-3}$ are respectively the dotted, dashed-dotted, solid black lines. The dashed green line corresponds to the peak of the rotation rate in the HMNS (see Figure 1).</figcaption></figure>

which, with $\Delta x = 100 \text{ m}$, gives $\lambda_{crit} = 1.55 \text{ km}$. Note that, in the fast-rotating limit ($\omega_A \ll \Omega$), we get $\lambda_{crit} = 2.4 \text{ km}$ at a $R_{TI} = 7 \text{ km}$. In the case where the geometry does not constrain the wavelength, the Tayler instability would be resolved by a hundred-meter resolution.

However, the unstable region of the magnetic field is limited by its vertical length $\approx 1 \text{ km}$ (Figure 7), and we therefore have to limit the wavelength of the Tayler instability. To check the impact of dissipation on a mode with a 1 km wavelength, we solve the local dispersion relation of Acheson & Gibbons (1978) with equal diffusivities, $v = \eta = \kappa$. We also assume the magnetic field to be increasing with radius and decreasing vertically, and we take the other quantities to be the same as in the simulation (see Appendix C for more details). Figure 8 shows that the growth rate in the limit of low dissipation is in fair agreement with $\sigma_{TI} \approx \omega_A$, which is consistent with the regime $\omega_A \gg \Omega$. The stability limit due to dissipation is found to be in good agreement with the criterion $\omega_A = k_{TI}^2 \eta$ (dashed line in Fig. 8). We
<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Growth rate of the Tayler instability as a function of dissipation, with $v = \eta = \kappa$. The red line shows the theoretical growth rate in the regime of $\Omega \ll \omega_A$. The dashed line is the stabilizing limit given by $k_{TI}^2\eta = \omega_A$.</figcaption></figure>

also note that even for a high dissipation like the neutrino viscosity $v_{HMNS} = 7.36 \times 10^7 \text{ cm}^2 \text{ s}^{-1}$, the growth rate is unchanged compared to the ideal MHD case. We would therefore expect the instability to grow in the absence of numerical limitations.

We thus use the limit given by Eq. (20) with the numerical dissipation as the criterion to see whether the Tayler instability can develop in the numerical simulation. Figure 9 shows the critical magnetic field for the numerical viscosity and resistivity in the code depending on the length scale below 1 km for different resolutions. The magnetic field found in the simulation is close to the critical magnetic field with the estimated dissipation of a 100-meter resolution, which may explain why the simulation is stable to Tayler instability at this resolution. Overall, in our understanding, the Tayler instability may be stabilized by the non-axisymmetric modes initially and later on by the resolution that is not enough to resolve the small region where it should be unstable.

# 5. Discussion

## 5.1. Comparison with Boussinesq simulations

Like the one-zone model of Barrère et al. (2022), the model of this article relies on strong assumptions to overcome the non-linearity of the dynamo mechanism, especially for the saturation phase (Spruit 2002; Zahn et al. 2007; Fuller et al. 2019). A comparison with the recent 3D numerical simulations of the TS dynamo with a positive differential rotation (Barrère et al. 2023, 2025) is therefore useful to test the validity of these necessary assumptions. Their simulations show the existence of several TS dynamo branches, but we only consider the strong, dipolar branch, which has been thoroughly studied by Barrère et al. (2025) and harbours the strongest magnetic fields. Furthermore, this solution is in good agreement with the scaling laws in Eqs. (31)–(33), derived by Fuller et al. (2019), which correspond to the formalism used in our one-zone model, with, however, a normalization factor significantly lower than 1.

In the region of the HMNS core with a positive shear (5–10 km), the stratification is in the range $N/\Omega \approx 0.5–1$, which is covered by the investigation of Barrère et al. (2025). By using the parameters of a typical HMNS ($R = 7$ km, $\rho =$

<figure><img src="image_14.png" /><figcaption>**Fig. 9.** Critical magnetic field $B_{\phi,crit}$ (Eq. (20)) in the numerical simulation depending on the wavelength of the Tayler instability for several resolutions with the numerical viscosity and resistivity estimated with the simulation grid size (Eq. (44)). The Tayler instability should be unstable in the simulation when the critical field is weaker than the magnetic field in the GRMHD simulation (gray region for $\Delta x = 12.5$ m).</figcaption></figure>

$3.7 \times 10^{14} \text{ g cm}^{-3}$, and $\Omega = 5468 \text{ s}^{-1}$) and extrapolating to $q \sim 1$, we infer from the numerical models of Barrère et al. (2025) with $N/\Omega = 0.5$ and $N/\Omega = 1$ that the axisymmetric toroidal, poloidal magnetic fields and the magnetic dipole should reach

$$B_{\text{tor}}^{m=0} \approx 0.63-1.4 \times 10^{16} \text{ G}, \qquad (46)$$

$$B_{\text{pol}}^{m=0} \approx 0.15-1.5 \times 10^{15} \text{ G}, \qquad (47)$$

$$B_{\text{dip}} \approx 0.79-6.3 \times 10^{14} \text{ G}, \qquad (48)$$

respectively. The axisymmetric toroidal magnetic field $B_{\text{tor}}^{m=0}$ is therefore weaker by a factor of ~20–45 than $B_{\phi}^{\text{F,sat}}$ (Eq. (31)). Likewise, the magnetic fields $B_{\text{tor}}^{m=0}$ and $B_{\text{dip}}$ are ~220–4200 times weaker than $B_{\text{r}}^{\text{F,sat}}$ (Eq. (32)), which is consistent with the fact that the TS dynamo in the simulations of Barrère et al. (2025) produces weaker large-scale magnetic fields than predicted by the formalism of Fuller et al. (2019). These intensities predicted by the scaling of MHD simulations are much lower than the one-zone model predicts.

These lower intensities are already reached in the GRMHD simulation we present in this study due to the high initial poloidal magnetic field. This can be interpreted in two different ways. First, we assume that the MHD simulation saturation gives the upper limit of the magnetic field in the GRMHD simulation. This means that, if the Tayler-Spruit dynamo were active in the GRMHD simulation, it would lead to a decrease in the magnetic field until it reaches the MHD saturation values, an evolution we do not see. In this case, the Tayler-Spruit dynamo would dissipate the axisymmetric toroidal magnetic field and prevent its growth to values $\ge 2 \times 10^{17} \text{ G}$. The second interpretation could be that the saturation intensities are different due to setup differences, as the background profiles, the differential rotation, and how it is generated in the MHD simulations could be quite important for the saturation values.

Moreover, these MHD simulations are likely to underestimate the magnetic field strength in a realistic HMNS. Indeed, for computational reasons, the simulations were run with much higher resistivities than realistic estimates, which hinders the dynamo process. Moreover, the quantities $B_{\text{tor}}^{m=0}$, $B_{\text{pol}}^{m=0}$, and $B_{\text{dip}}$
are volume-averaged. This suggests that the magnetic fields must be stronger locally when they are concentrated in some regions of the integration domain, which happens when the stratification increases.

## 5.2. MRI-driven dynamo vs TS dynamo

As it is argued in Margalit et al. (2022), the TS dynamo could also be relevant for the transport of angular momentum where the differential rotation is decreasing with radius. However, depending on the Brunt-Väisälä frequency, this region would be unstable to the MRI, which has a much faster growth rate than the TS dynamo in the fast-rotating case. For this reason, the MRI is expected to saturate first, and the Tayler instability would not impact this region. To confirm this picture, we compute the growth rate of both instabilities in the equatorial plane at r = 20 km with a magnetic field B₀ = 10¹⁵ G as it is found after the Kelvin-Helmholtz instability. We obtain

$$ 
\begin{aligned}
\sigma_{\text{MRI}} &= \frac{q\Omega}{2} \approx 2600 \text{ s}^{-1} \\
\sigma_{\text{TS}} &= \frac{\omega_A^2}{\Omega} \approx 0.072 \left(\frac{B_0}{10^{15} \text{ G}}\right)^2 \text{ s}^{-1} \ll \sigma_{\text{MRI}},
\end{aligned}
\quad (49) $$

due to $\omega_A \ll \Omega$. Interactions from the turbulence developed by the MRI and the TS dynamo would occur at the radius where the differential rotation changes sign, making this transition region difficult to study.

## 5.3. Angular momentum transport and spin-down time

To have the complete picture of the angular momentum transport in the HMNS, we first describe the transport in the outer region, R ≥ 10 km, which is unstable to the MRI. The angular momentum transport timescale due to MRI turbulence is

$$ t_{\text{MRI}} \approx \frac{R^2}{\alpha_{\text{turb}} H \Omega} \sim 5 \times 10^{-3} \left(\frac{\alpha_{\text{turb}}}{0.1}\right)^{-1} \left(\frac{\Omega}{\Omega_K}\right)^{-1} \left(\frac{H}{R/3}\right)^{-2} \quad [\text{s}]. \quad (50) $$

The spin-down in the outer region would be from the magnetic torque of the large-scale magnetic field generated by the $\alpha\Omega$ dynamo. On the other hand, our study shows that in the inner core, the angular momentum transport timescale due to the TS dynamo corresponds to the time it takes to reach saturation. Indeed, the very large saturated magnetic fields lead to very fast transport of angular momentum on a much shorter timescale than the dynamo growth timescale. This growth lies in the range $10^{-2} - 10^2$ s and scales with the initial magnetic field as $t_{\text{total}} = 65.7 (B_{R,0}/10^{10})^{-0.71}$ s. One interesting feature is that due to the positive differential rotation, the angular momentum would be transported inwards rather than outwards.

The spin-down time depends on how we consider the generated magnetic field inside the remnant. We first consider that the magnetic field stays buried in the core of the HMNS, then it would act as a torque and slow the core down following Eq. (5). By taking the saturated magnetic field values of the one-zone model (Eqs. 28 and 29), the core would slow down at a rate of

$$ \gamma_{\text{spin-down,core}} \equiv \frac{\dot{\Omega}}{\Omega} = -\frac{R_{\text{TI}}^3 B_{\text{R}} B_{\phi}}{I\Omega} \approx 3.5 \times 10^2 \text{ s}^{-1}. \quad (51) $$

We can also assume that the obtained dipole becomes the dipole at the surface of a remnant neutron star without matter outside.

Then, the dipole formula gives the following spin-down rate of

$$ \gamma_{\text{spin-down,dipole}} \equiv \frac{\dot{\Omega}}{\Omega} = -\frac{B_r^2 R_{\text{TI}}^6 \Omega^2}{6c^3 I} \approx 4.5 \times 10^{-2} \text{ s}^{-1} \quad (52) $$

assuming that the saturated value of the dipole is $B_r \sim B_R = 4 \times 10^{16}$ G. The spin-down due to the magnetic torque would therefore be faster than the dipole spin-down and slow the remnant down on a timescale of milliseconds.

Hence, the remnant would collapse in a few milliseconds when the dynamo saturates and the rotational kinetic energy would not be emitted in electromagnetic waves. If the magnetic braking follows the torque used in this study, this could therefore shorten the lifetime of a remnant neutron star after the merger depending on the initial radial field $B_{R,0}$. However, if we consider the saturated values from Barrère et al. (2025), the spin-down value due to the magnetic torque would be a thousand times slower $\gamma_{\text{spin-down,core}} \approx 0.33 \text{ s}^{-1}$. The collapse time would then be around $O(1)$ s after the dynamo saturation rather than $O(0.010)$ s.

# 6. Conclusions

We have investigated the TS dynamo in the context of BNS mergers. Following Barrère et al. (2022), we developed a one-zone model of the TS dynamo for a BNS merger remnant with realistic parameters estimated from a 3D GRMHD simulation.

We found that the TS dynamo could develop in the core of remnant neutron stars due to its positive differential rotation. Due to the impact of neutrino viscosity, the toroidal magnetic field must be amplified to higher values than $3.2 \times 10^{15}$ G by the winding in order to be Tayler-unstable.

The magnetic field's evolution can be divided into three phases: the winding phase, the Tayler instability phase and the non-linear growth of the dynamo. Saturation occurs when the differential rotation is effectively quenched by the Maxwell stress. The magnetic field generated by the dynamo saturates at a very high intensity of $B_\phi = 1.6 \times 10^{17}$ G and $B_R = 5.9 \times 10^{16}$ G for the magnetic dipole, according to the one-zone model. These strong intensities must be taken with caution as current simulations of the TS dynamo applied to PNSs predict a weaker toroidal field of $10^{16}$ G and a poloidal field of $10^{15}$ G.

The saturated magnetic field strength does not depend on the initial magnetic field, unlike the time required to reach saturation, which can range from 0.1 s for a magnetic dipole of $10^{14}$ G to 2.4 s for a magnetic dipole of $10^{12}$ G, which becomes comparable to the $O(1)$ s, which corresponds to the HMNS lifetime in the long-lived case. In the first case, the transport of angular momentum by the TS dynamo would lead to a faster collapse to a black hole or a faster spin-down of the remnant. This shows that the TS dynamo would be important in the case of a long-lived remnant neutron star as long as the initial magnetic field dipole is higher than $10^{12}$ G.

However, these results depend on whether the TS dynamo grows fast enough, and consequently on the resulting magnetic dipole after the Kelvin-Helmholtz instability. Using the study of Aguilera-Miret et al. (2024), we estimate the large-scale poloidal magnetic field to be of the order of the volume-averaged magnetic field $5 \times 10^{14}$ G, and the dynamo would saturate in less than 50 milliseconds. This would transport the angular momentum very efficiently in long-lived remnants and potentially reduce their lifetime, depending on the magnetic spin-down mechanism. Nonetheless, the volume-averaged magnetic field is different from the large-scale magnetic field. A more
detailed study is therefore needed, and it would require to use
realistic initial magnetic fields in the neutron stars and preferen-
tially have a converged growth of the Kelvin-Helmholtz instabil-
ity.

As we discussed in Section 5.1, the simulations by Barrère et al. (2025) give weaker saturated magnetic field strengths than our results but the differences could be due to a cylindrical rotation profile and the different background stratification. The value of the magnetic dipole needs to be confirmed by 3D numerical simulation with a rotation profile and a background stratification corresponding to a remnant neutron star. This is left to a further study, which should also address the question of the interaction between the TS dynamo and the MRI instability that operates at cylindrical radii larger than 10 km.

The TS dynamo operating in the remnant neutron star could
have significant consequences for astrophysical observations.
First, the angular momentum transport and the magnetic brak-
ing would happen when the dynamo saturates, which would lead
to the collapse of the remnant neutron star if it is more massive
than the maximum mass of a non-rotating neutron star M<sub>TOV</sub>.
Otherwise, it would lead to a stable, slow-rotating neutron star
with a strong magnetic field, a proto-magnetar. If the rotation
is damped as fast as described by the magnetic torque of the
saturated field, the lifetime of the remnant neutron star would
be reduced and could prevent the late-time emission of a fast-
rotating proto-magnetar. These results, in particular the rate of
magnetic braking, need to be tested with numerical simulations.

In addition, having motions due to the TS dynamo in the neu-
tron star remnant is expected to lead to some emission in gravita-
tional waves, which might be detectable for a close event. Using
Fuller's formalism, we estimate the amplitude of the velocity
perturbations due to TS dynamo with the formula

$$
\delta v_{r,sat} = \delta v_{A,sat} \frac{\omega_{A,sat}^2}{N\Omega} \approx 4.5 \times 10^9 \text{ cm s}^{-1}. \quad (53)
$$

For the sake of simplicity, we estimate the strain by using
the order of magnitude derived from the quadrupole formula
(Kokkotas 2002)

$$
h \approx \frac{G \varepsilon E_{\text{turb}}}{c^4 D_{\text{GW}}}, \qquad (54)
$$

where ε is the efficiency to convert the turbulent energy to grav-
itational waves that we assume to be low, with ε ≈ 5%. The
kinetic energy is computed by assuming that this velocity is con-
stant in the region from 5 to 10 km and oscillates at the Brunt-
Väisälä frequency N. A distance of 100 Mpc would then give a
strain of h ≈ 2 × 10⁻²⁵, which is of the same order but slightly
lower than the sensitivity of the future generation of gravita-
tional detectors. Knowing that the kinetic energy predicted here
can be viewed as optimistic, a direct detection of the signal will
be difficult but the TS dynamo would still impact the collapse
time of the remnant. Further investigations of the TS dynamo in
BNS mergers are therefore important to better understand future
multi-messenger observations.

Acknowledgements. This work was supported by the European Research Council (MagBURST grant 715368) and by the “Action thématique Phénomènes Extrêmes et Multi-messagers” (AT-PEM) of CNRS/INSU PN Astro with INP and IN2P3, co-funded by CEA and CNES. This work was in part supported by the Grant-in-Aid for Scientific Research (grant Nos. 23H04900, 23H01172, 23K25869) of Japan MEXT/JSPS.

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L13

Acheson, D. J., & Gibbons, M. P. 1978, Phil. Trans. R. Soc. London Ser. A, 289, 459

Aguilera-Miret, R., Palenzuela, C., Carrasco, F., Rosswog, S., & Viganò, D. 2024, Phys. Rev. D, 110, 083014

Barrère, P., Guilet, J., Reboul-Salze, A., Raynaud, R., & Janka, H. T. 2022, A&A, 668, A79

Barrère, P., Guilet, J., Raynaud, R., & Reboul-Salze, A. 2023, MNRAS, 526, L88

Barrère, P., Guilet, J., Raynaud, R., & Reboul-Salze, A. 2025, A&A, 695, A183

Baumgarte, T. W., Shapiro, S. L., & Shibata, M. 2000, ApJ, 528, L29

Bonanno, A., Guerrero, G., & Del Sordo, F. 2020, Mem. Soc. Astron. It., 91, 249

Breu, C., & Rezzolla, L. 2016, MNRAS, 459, 646

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

Combi, L., & Siegel, D. M. 2023, Phys. Rev. Lett., 131, 231402

Cook, G. B., Shapiro, S. L., & Teukolsky, S. A. 1994a, ApJ, 424, 823

Cook, G. B., Shapiro, S. L., & Teukolsky, S. A. 1994b, ApJ, 422, 227

Daniel, F., Petitdemange, L., & Gissinger, C. 2023, Phys. Rev. Fluids, 8, 123701

Denissenkov, P. A., & Pinsonneault, M. 2007, ApJ, 655, 1157

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

Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661

Gao, H., Ding, X., Wu, X.-F., Zhang, B., & Dai, Z.-G. 2013, ApJ, 771, 86

Giacomazzo, B., & Perna, R. 2013, ApJ, 771, L26

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

Goossens, M., & Tayler, R. J. 1980, MNRAS, 193, 833

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

Hempel, M., & Schaffner-Bielich, J. 2010, Nucl. Phys. A, 837, 210

Ji, S., Fuller, J., & Lecoanet, D. 2023, MNRAS, 521, 5372

Jouve, L., Lignières, F., & Gaurat, M. 2020, A&A, 641, A13

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

Kastaun, W., Ciolfi, R., Endrizzi, A., & Giacomazzo, B. 2017, Phys. Rev. D, 96, 043019

Keil, W., Janka, H. T., & Mueller, E. 1996, ApJ, 473, L111

Kiuchi, K., Yoshida, S., & Shibata, M. 2011, A&A, 532, A30

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

Kiuchi, K., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2018, Phys. Rev. D, 97, 124039

Kiuchi, K., Reboul-Salze, A., Shibata, M., & Sekiguchi, Y. 2024, Nat. Astron., 8, 298

Kokkotas, K. 2002, Gravitational wave Physics, https://phas.ubc.ca/~halpern/403/references/GW_Physics.pdf

Lasota, J.-P., Haensel, P., & Abramowicz, M. A. 1996, ApJ, 456, 300

Margalit, B., Jermyn, A. S., Metzger, B. D., Roberts, L. F., & Quataert, E. 2022, ApJ, 939, 51

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

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

Mösta, P., Radice, D., Haas, R., Schnetter, E., & Bernuzzi, S. 2020, ApJ, 901, L37

Müller, B., Janka, H.-T., & Marek, A. 2013, ApJ, 766, 43

Musolino, C., Ecker, C., & Rezzolla, L. 2024, ApJ, 962, 61

Petitdemange, L., Marcotte, F., & Gissinger, C. 2023, Science, 379, 300

Petitdemange, L., Marcotte, F., Gissinger, C., & Daniel, F. 2024, A&A, 681, A75

Pitts, E., & Tayler, R. J. 1985, MNRAS, 216, 139

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

Radice, D., & Bernuzzi, S. 2023, ApJ, 959, 46

Rasio, F. A., & Shapiro, S. L. 1999, Class. Quant. Grav., 16, R1

Raynaud, R., Guilet, J., Janka, H.-T., & Gastine, T. 2020, Sci. Adv., 6, eaay2732

Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2022, A&A, 667, A94

Rembiasz, T., Obergaulinger, M., Cerdá-Durán, P., Aloy, M.-A., & Müller, E. 2017, ApJS, 230, 18

Rüdiger, G., & Kitchatinov, L. L. 2010, Geophys. Astrophys. Fluid Dyn., 104, 273

Sarin, N., Omand, C. M. B., Margalit, B., & Jones, D. I. 2022, MNRAS, 516, 4949

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

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

Shibata, M., Fujibayashi, S., & Sekiguchi, Y. 2021, Phys. Rev. D, 104, 063026

Skoutnev, V. A., & Beloborodov, A. M. 2024, ApJ, 974, 290

Spruit, H. C. 2002, A&A, 381, 923

Szklarski, J., & Arlt, R. 2013, A&A, 550, A94

Tayler, R. J. 1973, MNRAS, 161, 365

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

Zahn, J. P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145
# Appendix A: Effects from general relativity

In the context of BNS merger remnant neutron stars, effects from general relativity can become important and the formalism needs to be adapted to a 3+1 decomposition. The three-component velocity $v^i$, which corresponds to the Newtonian velocity in the non-relativistic limit, is defined by

$$v^i = \frac{u^i}{u^t}, \qquad (A.1)$$

where $u^i$ is the four-component velocity, $u^t$, $u^j$ being its time and space components. In our simulation, the shift vector can be neglected at radii larger than 5 km. We can also assume that the movement of the flow is non-relativistic meaning that $u^t = W/\alpha \approx 1/\alpha$, where W is the Lorentz factor of the flow.

In ideal GRMHD, the induction equation is the same as the Newtonian case but for the magnetic field in the inertial frame $\mathcal{B}^i$, written as $B_i$ in the following, and the three-component velocity $v^i$. The magnetic field that enters the momentum equation is the magnetic field comoving with the fluid flow $b^i$ and the spatial components are defined as

$$b_i = \frac{1}{\alpha u^t} (B_i + u_i (B^j u_j)). \qquad (A.2)$$

Knowing that $u_\phi \sim \frac{1}{\alpha} v_\phi$ would be dominant the component of the velocity and that $v_\phi$ is of the order of magnitude 0.1c, the order of magnitude of the corrections would be 4% for $b_\phi$ as it would scale as $u_\phi^2$ and lower for the other components. Another correction would be in the expression of the angular momentum transport by the Maxwell stress tensor which is then given by

$$T_{R\phi}^{MAX} = b^2 u_R u_\phi - b_R b_\phi, \qquad (A.3)$$

where $b^2 = b^i b_\mu = \frac{1}{\alpha^2(u^i)^2} (B^2 + (B^i u_i)^2)$. The corrections would be of the order of $\frac{b_\phi}{b_R} u_R u_\phi$ and therefore around ~8% as $b_\phi/b_R \sim 2$.

In GRMHD, the definition of the Alfvén velocity also changes to

$$v_A^2 = \frac{b^2}{4\pi\rho h + b^2}, \qquad (A.4)$$

where h is the specific enthalpy. In the case of $b^2 \ll 4\pi\rho h$, the formalism can be simply adapted by replacing $\rho$ by $\rho h$. This assumption corresponds to $v_A \ll c$, which holds already in the case of $\omega_A \ll \Omega$ as the rotational velocity $R_{TI}\Omega$ is lower than c. This correction would reduce the Alfvén frequency and thus the growth rate $\sigma_{TI}$ by $h(R_{TI}) \sim 1.126c^2$ in the case of fast rotation.

# Appendix B: Force balance at high viscosity

Fuller et al. (2019) make the assumption of a Coriolis-Lorentz force balance for the perturbed velocity. However, another balance of the viscous forces and Lorentz force could be possible. To compare the Coriolis force and the viscous force, a local Ekman number can be defined as

$$E = \frac{\nu k_{TI}^2}{\Omega}, \qquad (B.1)$$

where $k_{TI}$ is the wavenumber of the displacements due to the Tayler instability. When $E < 1$, the Coriolis force is stronger than the viscous forces and a magnetostrophic balance can be safely assumed.

The Ekman number based on the conservative critical magnetic field $B_{\phi,crit}$ (Eq. (20)) and its wavelength (Eq. (22)) is $E \approx 1.6 \times 10^{-6} \ll 1$, so we can safely assume the magnetostrophic balance. From Eq. (18), we can estimate the range of magnetic fields where viscosity dominates over Coriolis force, corresponding to the regime $E > 1$, which is given by

$$B_{\phi} < B_{\phi,\text{min,visc}} = \sqrt{4\pi\rho N} \sqrt{\frac{\nu}{\Omega}} \approx 2.9 \times 10^{13} \text{G}. \qquad (B.2)$$

Therefore, the magnetostrophic balance can be assumed for both the Tayler instability and the dynamo phases, but not necessarily for the initial stages of the time evolution. However, the magnetic field strength is amplified to values greater than $B_{\text{min,visc}}$ by the winding, for which the magnetostrophic balance can also be assumed. Since we have $B_{\phi,\text{min,visc}} < B_{\phi,crit}$, the magnetostrophic balance can be assumed during the Tayler instability phase and dynamo phase. However, in the winding phase, $B_\phi$ could be smaller than $B_{\phi,\text{min,visc}}$ depending on the initial value of $B_R$. In such a case, a viscous-Lorentz force balance should be assumed. However, as this phase is not very important for the growth of the perturbed magnetic field and, therefore, we can assume a magnetostrophic balance for the whole duration. The viscous-Lorentz balance can be written as $\delta v_\perp \sim \delta v_A \omega_A / (\nu k_{TI}^2)$, which is derived from the magnetostrophic balance $\delta v_\perp \sim \delta v_A \omega_A / \Omega$ divided by the Ekman number. The magnetostrophic balance is used to derive the non-linear growth rate of the TS dynamo, and the balance needs to be divided by the Ekman number $E$ with the viscous-Lorentz. In case we wanted to adapt Eq. (27), the growth rate of $B_R$ would then be divided by $E$ if $E > 1$, which gives

$$\partial_t B_R = (\sigma_{\text{NL}} - \gamma_{\text{diss}}) B_R = \mathcal{B}_{\text{TI}} \frac{\omega_A^2}{N\Omega} \frac{\delta B_{\perp}^2}{\max(1, E) \sqrt{4\pi\rho R_{\text{TI}}^2}} - \frac{\omega_A^2}{\Omega} \left(\frac{\delta B_{\perp}}{B_{\phi}}\right)^2 B_R. \qquad (B.3)$$

# Appendix C: Linear study of Tayler instability

We solve the dispersion relation of Acheson & Gibbons (1978) with Python. We first use symPy to calculate the polynomial coefficients from the relation and then solve the roots of the polynomial numerically. We tested that we recover the dispersion relation in both the non-rotating ideal case and the rotating case as in Kiuchi et al. (2011). We can recover the theoretical growth rate in the two opposite limits $\omega_A \ll \Omega$ and $\Omega \ll \omega_A$. The code used to solve the linear dispersion relation is publicly available at the address in the footnote 2.

<footnote n="2">2 https://github.com/alexisreboulsalze/AchesonLinSolve/tree/main</footnote>