# A new scenario for magnetar formation: Tayler-Spruit dynamo in a proto-neutron star spun up by fallback

P. Barrère, Jérôme Guilet, Alexis Reboul-Salze, Raphaël Raynaud, Hans-Thomas Janka

► To cite this version:

P. Barrère, Jérôme Guilet, Alexis Reboul-Salze, Raphaël Raynaud, Hans-Thomas Janka. A new scenario for magnetar formation: Tayler-Spruit dynamo in a proto-neutron star spun up by fallback. Astronomy and Astrophysics - A&A, 2022, 668 (A79), pp.A79. 10.1051/0004-6361/202244172. hal-03699876

HAL Id: hal-03699876

https://hal.science/hal-03699876

Submitted on 7 Dec 2022

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
# A new scenario for magnetar formation: Tayler-Spruit dynamo in a proto-neutron star spun up by fallback

P. Barrère<fnref n="1" />, J. Guilet<fnref n="1" />, A. Reboul-Salze<fnref n="2" />, R. Raynaud<fnref n="1" />, and H.-T. Janka<fnref n="3" />

Received 2 June 2022 / Accepted 11 October 2022

## ABSTRACT

Magnetars are isolated young neutron stars characterised by the most intense magnetic fields known in the Universe, which power a wide variety of high-energy emissions from giant flares to fast radio bursts. The origin of their magnetic field is still a challenging question. In situ magnetic field amplification by dynamo action could potentially generate ultra-strong magnetic fields in fast-rotating progenitors. However, it is unclear whether the fraction of progenitors harbouring fast core rotation is sufficient to explain the entire magnetar population. To address this point, we propose a new scenario for magnetar formation involving a slowly rotating progenitor, in which a slow-rotating proto-neutron star is spun up by the supernova fallback. We argue that this can trigger the development of the Tayler-Spruit dynamo while other dynamo processes are disfavoured. Using the findings of previous studies of this dynamo and simulation results characterising the supernova fallback, we derive equations modelling the coupled evolution of the proto-neutron star rotation and magnetic field. Their time integration for different accreted masses is successfully compared with analytical estimates of the amplification timescales and saturation value of the magnetic field. We find that the magnetic field is amplified within 20–40 s after the core bounce, and that the radial magnetic field saturates at intensities between ~10¹<fnref n="3" /> and 10¹⁵ G, therefore spanning the full range of a magnetar’s dipolar magnetic fields. The toroidal magnetic field is predicted to be a factor of 10–100 times stronger, lying between ~10¹⁵ and 3 × 10¹⁶ G. We also compare the saturation mechanisms proposed respectively by H.C. Spruit and J. Fuller, showing that magnetar-like magnetic fields can be generated for a neutron star spun up to rotation periods of ≤8 ms and ≤28 ms, corresponding to accreted masses of ≥4 × 10⁻<fnref n="2" /> M☉ and ≥1.1 × 10⁻<fnref n="2" /> M☉, respectively. Therefore, our results suggest that magnetars can be formed from slow-rotating progenitors for accreted masses compatible with recent supernova simulations and leading to plausible initial rotation periods of the proto-neutron star.

**Key words.** stars: magnetars – supernovae: general – magnetohydrodynamics (MHD) – dynamo

## 1. Introduction

Magnetars represent two classes of isolated young neutron stars whose emission is powered by their ultrastrong magnetic field: anomalous X-ray pulsars and soft gamma repeaters. They feature a large spectrum of activity from short bursts (Gotz et al. 2006; Coti Zelati et al. 2018, 2021) to giant flares (Evans et al. 1980; Hurley et al. 1999, 2005; Svinkin et al. 2021), whose signal contains quasi-periodic oscillations (Israel et al. 2005; Strohmayer & Watts 2005; Gabler et al. 2018; Roberts et al. 2021). Moreover, a Galactic magnetar has recently been associated with a fast radio burst (FRB; Bochenek et al. 2020; CHIME/FRB Collaboration et al. 2020), which validates the capability of magnetar scenarios to explain at least a fraction of FRBs.

The pulsed X-ray activity of magnetars shows that they are characterised by a slow rotation period of 2–12 s and a fast spin-down. Under the assumption of a magnetic dipole spin-down, magnetars are therefore constrained to exhibit strong dipolar surface magnetic fields ranging from 10¹⁴ to 10¹⁵ G (Kouveliotou 1999; Kaspi & Beloborodov 2017), which are two orders of magnitude larger than in regular neutron stars. Furthermore, several lines of evidence suggest the presence of a non-dipolar magnetic field stronger than the dipolar component. Indeed, absorp-

tion lines have been detected in the X-ray spectra of two magnetars: SGR 0418+5729 (Tiendo et al. 2013) and SWIFT J1882.3-1606 (Rodríguez Castillo et al. 2016). If these lines are interpreted as proton cyclotron lines, they are respectively the signature of non-dipolar magnetic fields of ~2 × 10¹⁴–10¹⁵ G and ~6 × 10¹⁴–2.5 × 10¹⁵ G, which are stronger than their respective dipolar components by a factor of ~30–170 (Rea et al. 2010, 2012). Another sign of strong non-dipolar magnetic fields is the detection of a phase modulation in the hard-X-ray emission of a few magnetars. This may be explained by precession movements due to an internal toroidal magnetic field reaching a strength of ~10¹⁶ G (Makishima et al. 2014, 2016, 2019, 2021).

Proto-magnetars may be the central engine of extreme events if they are born rotating with a period of a few milliseconds. Indeed, their large-scale magnetic field can extract a large amount of rotational energy, which may create jets and lead to magnetorotational explosions (Burrows et al. 2007; Dessart et al. 2008; Takiwaki et al. 2009; Kuroda et al. 2020; Bugli et al. 2020, 2021; Obergaulinger & Aloy 2020, 2021, 2022). This process may explain hypernovae that are associated with long gamma-ray bursts (Duncan & Thompson 1992; Zhang & Mészáros 2001; Woosley & Bloom 2006; Drout et al. 2011; Nomoto et al. 2011; Gompertz & Fruchter 2017; Metzger et al. 2011, 2018). Moreover, their spin-down

<footnote n="1">¹ Université Paris-Saclay, Université Paris-Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: paul.barrere@cea.fr</footnote>

<footnote n="2">² Max Planck Institute for Gravitational Physics (Albert Einstein Institute), 14476 Potsdam, Germany</footnote>

<footnote n="3">³ Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany</footnote>
luminosity is invoked as a source of delayed energy injection to explain superluminous supernovae (SNe; Woosley 2010; Kasen & Bildsten 2010; Dessart et al. 2012; Inserra et al. 2013; Nicholl et al. 2013). Finally, millisecond magnetars, which may be formed in binary neutron star mergers, could also provide an explanation for the plateau phase in the X-ray emission of some short gamma-ray bursts (Metzger et al. 2008; Lü & Zhang 2014; Gompertz et al. 2014).

The central question to understand magnetar formation is the origin of their ultra-strong magnetic field. One type of scenario invokes magnetic flux conservation during the collapse of magnetised progenitors (Ferrario & Wickramasinghe 2006; Hu & Lou 2009). The magnetic field of these progenitors can originate from either a fossil field (Braithwaite & Spruit 2004, 2017) or dynamo action during main sequence star mergers (Schneider et al. 2019, 2020). While the surface magnetic field is constrained by observations (Petit et al. 2019), the magnetic field intensity in the iron core remains unknown, which makes this scenario uncertain. Another class of formation scenarios is the in situ amplification of the magnetic field by a dynamo process after the core collapse, especially at early stages of the proto-magnetar evolution. Two mechanisms have been studied so far: the convective dynamo (Thompson & Duncan 1993; Raynaud et al. 2020, 2022; Masada et al. 2022; White et al. 2022) and the magnetorotational instability (MRI)-driven dynamo (e.g. Obergaulinger et al. 2009; Mösta et al. 2014; Guilet & Müller 2015; Reboul-Salze et al. 2021a,b). The efficiency of these two dynamo mechanisms in the physical conditions relevant to a proto-neutron star (PNS) is still uncertain, in particular because the regime of very high magnetic Prandtl numbers (i.e. the ratio of viscosity to magnetic diffusivity) has not yet been thoroughly explored (Guilet et al. 2022; Lander 2021). Numerical simulations suggest that the efficiency of both dynamos increases for faster PNS rotation (Raynaud et al. 2020, 2022; Reboul-Salze et al. 2021a,b), which makes them good candidates to explain the central engine of extreme explosions. However, it may be more challenging for them to explain magnetar formation in standard SNe, which requires slower initial rotation of the PNS. Indeed, the observed SN remnants associated with Galactic magnetars have an ordinary kinetic explosion energy (Vink & Kuiper 2006; Martin et al. 2014; Zhou et al. 2019). This suggests that most Galactic magnetars are formed in standard SNe, which is consistent with the fact that extreme explosions represent about 1% of all SNe whereas magnetars constitute at least 10% of the whole Galactic young neutron-star population (Kouveliotou et al. 1994; Gill & Heyl 2007; Beniamini et al. 2019). Under the assumption that all the rotational energy of the PNS is injected into the kinetic energy of the explosion, the kinetic energy of the proto-magnetar must not exceed the standard kinetic energy of a SN explosion of $10^{51}$ erg, which translates into a constraint on its initial rotation period of $\ge 5$ ms (Vink & Kuiper 2006).

All things considered, the aforementioned scenarios require that the progenitor core be either strongly magnetised or fast rotating. It remains uncertain as to whether one of these conditions is met in a sufficient number of progenitors. This article presents our investigation of a new scenario wherein magnetars form from a slowly rotating, weakly magnetised progenitor. We consider the situation in which a newly formed PNS is spun up by the matter initially ejected by the SN explosion that remains gravitationally bound to the compact remnant and eventually falls back onto its surface. As the accretion is asymmetric, recent numerical simulations suggest that the fallback can bring a significant amount of angular momen-

tum to the PNS surface (Chan et al. 2020; Stockinger et al. 2020; Janka et al. 2022). We investigate the possibility that a magnetar may form due to the dynamo action triggered by the spin up from this fallback accretion. In this scenario, the MRI is expected to be stable because the PNS surface rotates faster than the core. The fallback starts roughly ~5–10 s after the core bounce (Stockinger et al. 2020; Janka et al. 2022), which may be too late for the development of a convective dynamo. Instead, we suggest that the magnetic field is amplified by another dynamo mechanism: the so-called Tayler-Spruit dynamo, which is driven by the Tayler instability. This instability feeds off a toroidal field in a stably stratified medium due to the presence of an electric current along the axis of symmetry (Tayler 1973; Pitts & Tayler 1985). Spruit (2002) proposed a first model of a dynamo driven by the Tayler instability in a differentially rotating stably stratified region. This model has received criticism from several authors (see Denissenkov & Pinsonneault 2007; Zahn et al. 2007), which has been addressed in the alternative description proposed by Fuller et al. (2019). The Tayler-Spruit dynamo has long been elusive in numerical simulations, but recent numerical simulations provide the first numerical evidence for its existence (Petitdemange et al. 2022). This dynamo is usually invoked for magnetic field amplification in the context of stellar interior physics, especially because of its suspected implications for angular momentum transport and the magnetic desert in Ap/Bp stars (e.g. Rüdiger & Kitchatinov 2010; Szklarski & Arlt 2013; Bonanno & Guarnieri 2017; Guerrero et al. 2019; Ma & Fuller 2019; Bonanno et al. 2020; Jouve et al. 2020). However, this dynamo process has never been studied in the framework of magnetar formation.

In the following, Sect. 2 presents the scenario in more detail and the formalism used by our model. We describe our results in Sect. 3 and discuss them in Sect. 4. Finally, we draw conclusions in Sect. 5.

## 2. Mathematical modelling of the scenario

To study our scenario, we built a one-zone model consisting in 'average' time evolution equations that capture the main stages sketched in Fig. 1. We start by describing the impact of the SN fallback on the PNS rotation (the differential rotation) and the magnetic field (the shearing of the radial field and the exponential growth of the Tayler instability). Finally, we present the mathematical formalism for the non-linear stages, that is, the saturation mechanism of the dynamo as modelled by Spruit (2002) and Fuller et al. (2019), which we complete by a description of the generation of the radial magnetic field through non-linear induction. For the computation of the time evolution, we only implement the description based on the work of Fuller et al. (2019); but in Sect. 3 we compare both models regarding the predictions of the saturated magnetic field.

### 2.1. Fallback accretion

Our scenario starts a few seconds after the core bounce when a fraction of the fallback matter gets accreted onto the PNS surface. This matter is initially ejected during the explosion but stays gravitationally bound to the PNS, and so begins to be asymmetrically accreted (Chan et al. 2020). This fallback matter is thought to have a large angular momentum, which can even reach the magnitude of the Keplerian angular momentum (Janka et al. 2022). Therefore, the spin of the PNS is strongly affected and the surface rotation can be accelerated up to millisecond periods. In our scenario, the core of the progenitor is
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of the different stages of our magnetar formation scenario. The dashed line encloses the region of the fallback (orange arrows). Red and white lines represent the magnetic field lines and fluid motions, respectively. Ω and E$_{φ}$ stand for the angular rotation frequency and the azimuthal component of the electromotive force, respectively. B$_{φ}$ and B$_{r}$ are the axisymmetric azimuthal and radial magnetic fields, and δB$_{⊥}$ is the non-axisymmetric perpendicular magnetic field.</figcaption></figure>

assumed to be slowly rotating. Thus, the PNS surface spins faster than the PNS interior, which creates differential rotation.

To model the accretion onto the PNS surface, we use the asymptotic scaling for the mass accretion rate $\dot{M}_{acc} \propto t^{-5/3}$ from Chevalier (1989). As the accretion mass rate must be finite at the beginning of this accretion regime, we define a start time $t_0$ such that

$$ \dot{M}_{acc} = \frac{A}{(t + t_0)^{5/3}}, \qquad (1) $$

where A is a constant. Then, the accreted mass during this regime is

$$ M_{acc} = \int_{0}^{\infty} \frac{A}{(t + t_0)^{5/3}} dt. \qquad (2) $$

As $M_{acc}$ is constant, we have $A = \frac{2}{3}t_0^{2/3}M_{acc}$ and so the accretion mass rate is

$$ \dot{M}_{acc} = \frac{2}{3} M_{acc} \frac{t_0^{2/3}}{(t + t_0)^{5/3}}. \qquad (3) $$

From the fallback matter, only a fraction with angular momentum as large as the Keplerian limit at most will be accreted by the PNS, as discussed by Janka et al. (2022). Therefore, the relation between the average angular rotation frequency of the PNS and the mass accretion rate is

$$ \dot{\Omega} = \frac{\dot{j}_{\text{kep}}}{I} \dot{M}_{\text{acc}}, \qquad (4) $$

where I stands for the PNS moment of inertia and $\dot{j}_{\text{kep}} = \sqrt{GM_{\text{PNS}}r}$ is the specific Keplerian angular momentum at the PNS surface. As the PNS mass changes little and the contraction of the PNS is almost over at the times considered for the fallback accretion, we assume I to be constant. As supposed in

Fuller et al. (2019), the angular momentum is transported faster latitudinally than radially due to stratification, meaning that the differential rotation is shellular, that is Ω is constant on spherical shells.

As the accretion process spins up only the outer part of the PNS but not its inner core, the shear rate $q \equiv r\partial_r \ln \Omega$ is also expected to evolve. To describe this effect, we use the approximate expression

$$ r\partial_r\Omega \sim \Omega - \Omega(r = 0), \qquad (5) $$

where Ω and Ω(r = 0) are the average and central angular rotation frequency, respectively. Assuming that the rotation frequency at the centre of the PNS is unchanged by the accretion process (it will change only due to angular momentum transport processes described in Sect. 4.1), we infer the time derivative of the shear rate as

$$ \dot{q} \sim \frac{\dot{\Omega}}{\Omega}(1 - q). \qquad (6) $$

## 2.2. Shearing and Tayler instability growth

The differential rotation generated by the fallback will shear the radial component of the large-scale radial magnetic field $B_r$ into the azimuthal field $B_φ$ as follows

$$ \partial_t B_\phi = q \Omega B_r. \qquad (7) $$

Therefore, we can define a growth rate<fnref n="1" /> for $B_φ$:

$$ \sigma_{\text{shear}} \equiv q \frac{B_r}{B_\phi}. \qquad (8) $$

<footnote n="1">¹ We note that this growth rate and the others defined below are rather instantaneous growth rates because they depend on the magnetic field.</footnote>
As $B_\phi$ grows, it becomes Tayler unstable. To depict the linear growth of the instability, we make the following assumptions. First, the stratified medium of the PNS interior is characterised by the Brunt-Väisälä frequency (Hüdepohl 2014):

$$N \equiv \sqrt{-\frac{g}{\rho} \left( \left. \frac{\partial \rho}{\partial S} \right|_{P,Y_e} \frac{dS}{dr} + \left. \frac{\partial \rho}{\partial Y_e} \right|_{P,S} \frac{dY_e}{dr} \right)} \sim 4 \times 10^3 \text{s}^{-1}, \quad (9)$$

where g, ρ, Y_e, and S are the gravitational acceleration, the PNS mean density, the electron fraction, and the entropy, respectively. In the remainder of this paper, we use the fiducial value N = 4 × 10³ s⁻¹ based on the results of the 1D core-collapse supernova (CCSN) simulations from Hüdepohl (2014, Chap. 5). The Brunt-Väisälä frequency is almost uniform in most of the PNS except near the surface where it peaks at ~10⁴ s⁻¹. Hüdepohl (2014) made a comparison between two different equations of state (EOS): Shen (Shen et al. 1998a,b, 2011) and LS220 (Lattimer & Swesty 1991), and found that the choice of the EOS mainly affects the localisation and duration of the convection but not the value of the Brunt-Väisälä frequency in the stably stratified region. Second, the main background azimuthal field is B_φ, which is associated with the Alfvén frequency:

$$\omega_A \equiv \frac{B_\phi}{\sqrt{4\pi\rho r^2}} \approx 11.6 \left( \frac{B_\phi}{10^{15} \text{ G}} \right) \text{s}^{-1}, \quad (10)$$

for r = 12 km and ρ = 4.1 × 10¹⁴ g cm⁻³. Finally, the frequencies characterising the PNS are ordered such that

$$N \gg \Omega \gg \omega_A. \quad (11)$$

The development of the Tayler instability is triggered after reaching the critical strength (Spruit 1999, 2002; Zahn et al. 2007)

$$\begin{aligned} B_\phi > B_{\phi,c} &\sim \Omega \left(\frac{N}{\Omega}\right)^{1/2} \left(\frac{\eta}{r^2 \Omega}\right)^{1/4} \sqrt{4\pi\rho r^2} \\ &\approx 2.5 \times 10^{12} \left(\frac{\Omega}{200\pi \text{ rad s}^{-1}}\right)^{1/4} \text{G}, \end{aligned} \quad (12)$$

where η ~ 10⁻⁴ cm² s⁻¹ (Thompson & Duncan 1993) is the magnetic diffusivity. The fastest-growing perturbations are the m = 1 modes with an associated rate of (Ma & Fuller 2019)

$$\sigma_{\text{TI}} \sim \frac{\omega_A^2}{\Omega} \approx 0.21 \left( \frac{B_\phi}{10^{15} \text{ G}} \right)^2 \left( \frac{\Omega}{200\pi \text{ rad s}^{-1}} \right)^{-1} \text{s}^{-1}. \quad (13)$$

As the PNS interior is strongly stratified, we can determine a maximum radial length scale for the instability

$$l_r \sim \frac{\omega_A}{N} l_\perp, \quad (14)$$

where the horizontal length scale is approximated by $l_\perp \sim r$.

## 2.3. Spruit's picture of the dynamo

Spruit (2002) proposes that the energy in the azimuthal large-scale field B_φ cascades to small scale, that is the form of the non-linear magnetic energy dissipation is

$$\dot{E}_{\text{mag}} \sim \gamma_{\text{turb}} |B_{\phi}|^2, \quad (15)$$

where γ_turb is the turbulent damping rate. To determine this rate, Spruit (2002) argues that the saturation of the instability occurs when the turbulent velocity field generates a sufficiently large

effective turbulent diffusivity to balance the growth rate of the instability, that is

$$\gamma_{\text{turb}} \sim \frac{\eta_e}{l_r^2} \sim \sigma_{\text{TI}}, \quad (16)$$

where η_e is an effective turbulent diffusivity.

The solenoidal character of the perturbed magnetic field implies $B_r/l_r \sim B_φ/l_\perp$, which leads to

$$B_r \sim \frac{\omega_A}{N} B_\phi, \quad (17)$$

using the relation between length scales of the instability given by Eq. (14). As the azimuthal magnetic field B_φ is generated via the shear of the radial magnetic field B_r, the dynamo is expected to saturate when the shear (Eq. (8)) balances the turbulent damping (Eq. (16)). Thus, the amplitudes of the magnetic field components saturate at

$$B_{\phi,S}^{\text{sat}} \sim \sqrt{4\pi\rho r^2} q \frac{\Omega^2}{N}, \quad (18)$$

$$B_{r,S}^{\text{sat}} \sim \sqrt{4\pi\rho r^2} q^2 \frac{\Omega^4}{N^3}. \quad (19)$$

This description of the dynamo mechanism has been criticised for two reasons: First of all, if the large-scale component of B_φ remains constant on larger length scales than l_r, the displacements produced by the instability are not expected to mix the large-scale field lines to damp B_φ through reconnection. Therefore, the damping rate estimated in Eq. (16) is overestimated for the large-scale components of the azimuthal field B_φ (see Fuller et al. 2019). Secondly, as m = 1 modes are dominant, the radial magnetic field B_r produced by the instability is non-axisymmetric, and therefore its shear generates a mostly non-axisymmetric azimuthal field B_φ. Hence, the axisymmetric component of the fields B_r and B_φ may not be related by Eq. (17) (see Zahn et al. 2007).

## 2.4. A revised model of the dynamo

This section presents a description of the dynamo that completes the model proposed by Fuller et al. (2019) in the sense that we consider the time evolution of the magnetic field. A clear distinction is now made between the 'axisymmetric' components B_φ, B_r, and the 'non-axisymmetric' perturbed components δB_⊥, δB_r, which become the ones connected by the solenoidal condition

$$\delta B_r \sim \frac{\omega_A}{N} \delta B_\perp. \quad (20)$$

To overcome the previously raised difficulties, Fuller et al. (2019) argue that the energy in the perturbed field δB dissipates to small scales and find that the damping rate is

$$\gamma_{\text{cas}} \sim \frac{\delta v_A}{r}, \quad (21)$$

where δv_A ≡ δB_⊥ / √(4πρ) is the perturbed Alfvén velocity. Thus, equating the instability growth rate (Eq. (13)) and the damping rate (Eq. (21)) gives the saturation strength of the perturbed field δB_⊥ for a given strength of azimuthal field B_φ,

$$\delta B_\perp \sim \frac{\omega_A}{\Omega} B_\phi. \quad (22)$$
When the instability is saturated, the non-linear magnetic energy dissipation is then

$$ \dot{E}_{\text{mag}} \sim \gamma_{\text{cas}} |\delta \mathbf{B}|^2 \sim \frac{\delta v_A}{r} |\delta B_{\perp}|^2. \qquad (23) $$

As the azimuthal field $B_\phi$ is the dominant magnetic component, $\dot{E}_{\text{mag}} \sim B_\phi \partial_t B_\phi$. Hence, a damping rate can be defined for the axisymmetric components $B_\phi$ and $B_r$:

$$ \gamma_{\text{diss}} \equiv \frac{\dot{E}_{\text{mag}}}{B_{\phi}^2}. \qquad (24) $$

As the previous expression of the magnetic energy (Eq. (23)) is only valid when the instability saturates, we use the expression

$$ \dot{E}_{\text{mag}} \sim \frac{\omega_A^2}{\Omega} |\delta B_{\perp}|^2, \qquad (25) $$

which is valid in both the saturated and non-saturated states. Therefore, the damping rate defined in Eq. (24) becomes

$$ \gamma_{\text{diss}} \sim \frac{\omega_A^2}{\Omega} \left( \frac{\delta B_{\perp}}{B_{\phi}} \right)^2. \qquad (26) $$

To close the dynamo loop, the Tayler instability must generate the axisymmetric radial magnetic field $B_r$ (α-effect), which will be sheared again (Ω-effect). In the framework of the mean field theory, the induction equation reads

$$ \partial_t \langle \mathbf{B} \rangle = \nabla \times (\langle \mathbf{v} \rangle \times \langle \mathbf{B} \rangle + \mathbf{E}) - \eta \Delta \langle \mathbf{B} \rangle, \qquad (27) $$

in which we ignore the resistive term. Considering the average symbol $\langle \cdot \rangle$ as an azimuthal average, we note $\langle \mathbf{B} \rangle = \mathbf{B}$ in order to remain consistent with our notation of the axisymmetric magnetic field. The electromotive force $\mathbf{E} \equiv \langle \delta \mathbf{v} \wedge \delta \mathbf{B} \rangle$ is the important non-linear quantity responsible for the generation of the axisymmetric radial field $B_r$. In spherical coordinates, the radial component of Eq. (27) is

$$ \partial_t B_r = \frac{1}{r \sin \theta} \left[ \partial_{\theta}(\sin \theta E_{\phi}) - \partial_{\phi} E_{\theta} \right]. \qquad (28) $$

As $B_r$ is axisymmetric, $E_\theta$ can be ignored. By definition, the azimuthal component of electromotive force is

$$ E_{\phi} = \delta v_r \delta B_{\theta} - \delta v_{\theta} \delta B_r. \qquad (29) $$

Supposing an incompressible perturbed velocity field and using Eq. (20), we write

$$ \frac{\delta v_r}{\delta v_{\perp}} \sim \frac{\delta B_r}{\delta B_{\perp}} \sim \frac{\omega_A}{N}, \qquad (30) $$

and so the azimuthal electromotive force reads

$$ E_{\phi} \sim \delta v_{\theta} \delta B_{r} \sim \delta v_{r} \delta B_{\theta} \sim \delta v_{r} \delta B_{\perp}, \qquad (31) $$

where we assume that the two terms on the right-hand side of Eq. (29) do not cancel. The production of the radial field $B_r$ can be approximated by

$$ \partial_t B_r \sim \frac{E_{\phi}}{r} \sim \frac{\delta v_r \delta B_{\perp}}{r}. \qquad (32) $$

We must note that this expression differs from Eq. (29) in Fuller et al. (2019), which appears to contain a typo. As in

Fuller et al. (2019) and Ma & Fuller (2019), we expect magnetostrophic balance $\delta v_\perp \sim \delta v_A \omega_A / \Omega$, which leads to

$$ \delta v_r \sim \frac{\omega_A^2}{N\Omega} \delta v_A. \qquad (33) $$

Thus,

$$ \partial_t B_r \sim \frac{E_\phi}{r} \sim \frac{\omega_A^2}{N\Omega} \frac{\delta B_{\perp}^2}{\sqrt{4\pi\rho r^2}}, \qquad (34) $$

and we can define a growth rate for $B_r$

$$ \sigma_{\text{NL}} \equiv \frac{1}{\sqrt{4\pi\rho r^2}} \frac{\omega_A^2}{N\Omega} \frac{\delta B_{\perp}^2}{B_r}. \qquad (35) $$

The radial field $B_r$ will saturate when its non-linear growth rate (Eq. (35)) is balanced by the turbulent dissipation (Eq. (26)). This way, we find the relation between the axisymmetric fields

$$ B_r \sim \frac{\omega_A}{N} B_\phi, \qquad (36) $$

using Eq. (22). We note that this relation is similar to Eq. (17) from Spruit (2002), which was derived for the non-axisymmetric components only. Fuller et al. (2019) also established the same relation arguing that the Tayler instability cannot operate when the magnetic tension forces become larger than the magnetic pressure forces leading to the instability.

The azimuthal magnetic field saturates when the shear rate (Eq. (8)) balances the dissipation rate (Eq. (26)). Thus, using the relations between the magnetic field components (Eqs. (22) and (36)), we find the magnetic field strengths in the saturated regime derived in Fuller et al. (2019):

$$ B_{\phi,F}^{\text{sat}} \sim \sqrt{4\pi\rho r^2 \Omega} \left(\frac{q\Omega}{N}\right)^{1/3}, \qquad (37) $$

$$ \delta B_{\perp,F}^{\text{sat}} \sim \sqrt{4\pi\rho r^2 \Omega} \left(\frac{q\Omega}{N}\right)^{2/3}, \qquad (38) $$

$$ B_{r,F}^{\text{sat}} \sim \sqrt{4\pi\rho r^2 \Omega} \left(\frac{q^2\Omega^5}{N^5}\right)^{1/3}. \qquad (39) $$

Finally, the angular momentum is redistributed in the PNS through Maxwell stresses associated with an effective angular momentum diffusivity (Spruit 2002; Fuller et al. 2019):

$$ v_{\text{AM}} = \frac{B_r B_\phi}{4\pi\rho q \Omega}, \qquad (40) $$

which affects the shear parameter at the rate

$$ \gamma_{\text{AM}} \equiv \frac{v_{\text{AM}}}{r^2}. \qquad (41) $$

## 2.5. Governing evolution equations

Now that the main equations involved in our scenario have been brought out, we can write the evolution equations for the rotation properties and the magnetic field. The evolution of PNS angular rotation frequency is driven by the fallback accretion rate (Eq. (3)) as described by Eq. (4). Hence,

$$ \dot{\Omega} = \frac{2}{3} \Delta\Omega \frac{t_0^{2/3}}{(t+t_0)^{5/3}}, \qquad (42) $$
where $\Delta\Omega = \Omega_{fin} - \Omega_{init} = M_{acc}j_{kep}/I$. As previously mentioned, the shear rate is also expected to decrease due to angular momentum transport (Eq. (41)) such that

$$ \dot{q} = \frac{\dot{\Omega}}{\Omega}(1-q) - \gamma_{AM}q = \frac{2}{3} \frac{\Delta\Omega}{\Omega} \frac{t_0^{2/3}}{(t+t_0)^{5/3}} - \frac{B_r B_\phi}{4\pi\rho\Omega r^2}. \quad (43) $$

Combining the different growth and damping rates given by Eqs. (8), (13), (21), (26), and (35), we find that the magnetic field evolution is governed by the following equations:

$$ \partial_t B_\phi = (\sigma_{\text{shear}} - \gamma_{\text{diss}}) B_\phi = q\Omega B_r - \frac{\omega_A^2}{\Omega} \frac{\delta B_\perp^2}{B_\phi}, \quad (44) $$

$$ \partial_t \delta B_\perp = (\sigma_{\text{TI}} - \gamma_{\text{eas}}) \delta B_\perp = \frac{\omega_A^2}{\Omega} \delta B_\perp - \frac{\delta v_A}{r} \delta B_\perp, \quad (45) $$

$$ \partial_t B_r = (\sigma_{\text{NL}} - \gamma_{\text{diss}}) B_r = \frac{\omega_A^2}{N\Omega} \frac{\delta B_\perp^2}{\sqrt{4\pi\rho r^2}} - \frac{\omega_A^2}{\Omega} \left( \frac{\delta B_\perp}{B_\phi} \right)^2 B_r. \quad (46) $$

Equations (42)–(46) are solved for a typical PNS of 5–10 s in age using the odeint function from the Python package SciPy. The PNS mass and radius are $M_{PNS} = 1.5 M_\odot$ and $r = 12$ km, and so the average density is $\rho = 4.1 \times 10^{14}$ g cm$^{-3}$. The moment of inertia is estimated using Eq. (12) from Lattimer & Schutz (2005):

$$ I = 0.237 M_{\text{PNS}} r^2 \left[ 1 + 4.2 \left( \frac{M_{\text{PNS}}}{M_{\odot}} \frac{1 \text{ km}}{r} \right) + 90 \left( \frac{M_{\text{PNS}}}{M_{\odot}} \frac{1 \text{ km}}{r} \right)^4 \right] \\ \approx 1.6 \times 10^{45} \text{ g cm}^2. \quad (47) $$

The PNS core is assumed to be initially in solid-body rotation (i.e. $q = 0$) and slowly rotating with an angular rotation frequency $\Omega_{init} = 2\pi \text{ rad s}^{-1}$ (i.e. an initial rotation period $P_{init} \equiv 2\pi/\Omega_{init} = 1$ s). We note that the results of the time integration weakly depend on these parameters as long as $P_{init} \ge 40$ ms. The parameters of the fallback are chosen to be consistent with the simulations of Janka et al. (2022), with a starting time at $t_0 = 7$ s. The initial magnetic field components $B_r$, $B_\phi$, and $\delta B_\perp$ are initialised at a strength of $10^{12}$ G, which is the typical magnetic field amplitude of regular neutron stars.

## 3. Results

We proceed with a twofold presentation of our model outputs: First we present the time evolution, in which we derive analytical predictions for the timescales of the different phases of the scenario and compare them to the integrated evolution. We then present the saturated regime where we focus on the magnetic field intensities reached via the Tayler-Spruit dynamo and provide an ‘upper’ limit of PNS rotation period (i.e. a ‘lower’ limit of fallback mass) to form magnetars.

### 3.1. Time evolution of the magnetic field

The time series for an asymptotic rotation period $P_{fin} \equiv 2\pi/\Omega_{fin} = 10$ ms displayed in Fig. 2 can be split into several phases, which are illustrated by the schematics at the top of the figure:

(i) $B_\phi$ is strongly amplified for ~4 s due to the winding of $B_r$ while the other components stay constant. As the mass-accretion rate is higher in this phase, strong increases in the shear rate and the rotation rate are also noted.

(ii) The Tayler instability develops and amplifies $\delta B_\perp$ for ~8 s.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Time evolution of the different components of the magnetic field (top), the dimensionless shear rate, and the angular rotation frequency (bottom) for an accreted mass $M_{acc} = 3.2 \times 10^{-2} M_\odot$ corresponding to an asymptotic rotation period $P_{fin} = 10$ ms. The different stages of the dynamo process are highlighted by the schematics at the top and their associated timescales by the double arrows. Their ends are illustrated by the dotted vertical lines: winding (black), linear development of the Tayler instability (dark blue), and the whole dynamo loop (red). The horizontal dashed lines (blue, orange, and green) show respectively the saturation intensities $B_{\phi,F}^{\text{sat}}$, $\delta B_{\perp,F}^{\text{sat}}$, and $B_{r,F}^{\text{sat}}$ (Eqs. 37–39) evaluated at the time of saturation. The dashed violet horizontal line represents the asymptotic angular rotation frequency $\Omega_{fin}$.</figcaption></figure>

(iii) The axisymmetric radial field $B_r$ is generated allowing the dynamo action to occur for ~5 s. The azimuthal magnetic field saturates at ~17 s, which corresponds to ~24 s after the core bounce.

(iv) The angular momentum transport by the magnetic field, which was negligible during the first ~16 s, becomes significant as the magnetic field saturates. This stage is discussed in Sect. 4.1.

The evolution of the magnetic field is shown until the whole angular momentum is transported, (i.e. when the shear rate reaches $q = 0$) at $t \sim 17.5$ s. Further evolution is not considered because our set of equations does not intend to describe either the relaxation phase of the magnetic field to a stable geometric configuration or the dynamics with very low shear where one would expect the Tayler-Spruit dynamo to stop or to act intermittently (Fuller & Lu 2022).

The different vertical lines in Fig. 3 show that the above phases occur at different times for different accreted masses. To better quantify the dependence on this parameter, we derive analytical estimates of the corresponding characteristic timescales $\tau_{shear}$, $\tau_{TI}$, and $\tau_{dyn}$ (see Fig. 2). First, the shearing phase begins
<figure><img src="image_4.png" /><figcaption>Fig. 3. Same as Fig. 2 but for total accreted masses of $M_{acc} = \{0.54, 0.8, 1.6, 6.4\} \times 10^{-2} M_{\odot}$ (corresponding to $P_{fin} = \{60, 40, 20, 5\}$ ms, respectively).</figcaption></figure>

when the fallback matter starts to be accreted on the PNS surface and finishes when the azimuthal magnetic field $B_{\phi}$ is strong enough to make the Tayler instability grow as fast as $B_{\phi}$, that is when the growth rate of the instability (Eq. (13)) is equal to the winding rate (Eq. (8)). Thus, the Alfvén frequency associated with the intensity of the azimuthal magnetic field at the end of this phase is

$$ \omega_A \sim \omega_{A,TI} \equiv (q\Omega^2\omega_{r,0})^{1/3}, \quad (48) $$

where $\omega_{r,0} \equiv B_r(t=0)/\sqrt{4\pi\rho r^2}$. Therefore, a characteristic timescale for the shearing stage can be defined as the inverse of the winding growth rate (Eq. (13)) evaluated at $\omega_A = \omega_{A,TI}$:

$$ \tau_{\text{shear}} \equiv \sigma_{\text{shear}}^{-1} \sim \frac{\omega_{A,TI}}{q\Omega\omega_{r,0}} = (q^2\omega_{r,0}^2\Omega)^{-1/3}. \quad (49) $$

Second, as the azimuthal field becomes unstable, the Tayler instability grows exponentially until the perturbed field reaches

the saturation intensity of Eq. (22). The perturbed field at saturation can be approximated by

$$ \delta B_{\perp}(t = t_{\text{sat}}) \sim \delta B_{\perp}(t = \tau_{\text{shear}}) \exp[\sigma_{\text{TI}}(t_{\text{sat}} - \tau_{\text{shear}})], \quad (50) $$

and so a characteristic timescale for this stage can be defined as

$$ \tau_{\text{TI}} \equiv t_{\text{sat}} - \tau_{\text{shear}} \sim \sigma_{\text{TI}}^{-1} \ln \left[ \frac{\delta B_{\perp}(t = t_{\text{sat}})}{\delta B_{\perp}(t = \tau_{\text{shear}})} \right]. \quad (51) $$

Using Eq. (22), we have

$$ \delta B_{\perp}(t = t_{\text{sat}}) \sim \frac{\omega_A(t = t_{\text{sat}})}{\Omega} B_{\phi}(t = t_{\text{sat}}). \quad (52) $$

In order to obtain a simple estimate, we make the rough approximations that $\delta B_{\perp}(t = \tau_{\text{shear}}) \sim \delta B_{\perp}(t = 0)$ and $B_{\phi}(t = \tau_{\text{shear}}) \sim B_{\phi}(t = t_{\text{sat}})$, which leads to

$$ \tau_{\text{TI}} \sim \frac{\Omega}{\omega_{A,TI}^2} \ln \left( \frac{\omega_{A,TI}^2 r}{\Omega \delta v_{A,0}} \right), \quad (53) $$
where $\delta v_{A,0} \equiv \delta v_A(t=0)$.

Third, when the perturbed field reaches a sufficient amplitude, the axisymmetric radial field is amplified through nonlinear induction, thus closing the dynamo loop. This phase ends when the magnetic field saturates at the intensities given by Eqs. (37)–(39). Likewise, we estimate the critical strength of the azimuthal field $B_{\phi,dyn}$ above which the dynamo loop is triggered by equating the growth rate of the radial field $B_r$ (Eq. (35)) and the winding rate (Eq. (8)). We obtain the Alfvén frequency associated with $B_{\phi,dyn}$

$$ \omega_{A,dyn} \equiv (qN\Omega^4\omega_{r,0}^2)^{1/7}, \qquad (54) $$

making use of Eq. (22). We define the dynamo characteristic timescale as

$$ \tau_{\text{dyn}} \equiv \left( \frac{B_{\phi}}{\partial_t^2 B_{\phi}} \right)^{1/2} . \qquad (55) $$

The time derivative of the radial magnetic field is

$$ \partial_t B_r \sim \frac{\omega_A^2 \delta v_A}{N \Omega r} \delta B_{\perp} \sim \frac{\omega_A^3 \delta v_A}{N \Omega^2 r} B_{\phi} \sim \frac{\omega_A^5}{N \Omega^3} B_{\phi}, \qquad (56) $$

using Eq. (22). Therefore,

$$ \partial_t^2 B_{\phi} \sim q\Omega\partial_t B_r \sim \frac{q\omega_A^5}{N\Omega^2} B_{\phi}, \qquad (57) $$

where $q$ and $\Omega$ are assumed to be constant during this phase. Thus, the dynamo characteristic timescale can be approximated as

$$ \tau_{\text{dyn}} \equiv \left( \frac{B_{\phi}}{\partial_t^2 B_{\phi}} \right)^{1/2} \sim \frac{\Omega}{\omega_{A,\text{dyn}}^2} \left( \frac{N}{q\omega_{A,\text{dyn}}} \right)^{1/2} . \qquad (58) $$

For the case $P_{\text{fin}} = 10$ ms, we have $\tau_{\text{shear}} \approx 2.3$ s, $\tau_{\text{TI}} \approx 8.5$ s, and $\tau_{\text{dyn}} \approx 5$ s, which are similar to the timescales illustrated by the dotted vertical lines in Fig. 2. The same observation can be made for $P_{\text{fin}} \le 30$ ms in Fig. 3. However, for $P_{\text{fin}} = 40$ ms ($M_{\text{acc}} = 0.008 M_{\odot}$) and $P_{\text{fin}} = 60$ ms ($M_{\text{acc}} = 0.0054 M_{\odot}$), the dynamo loop phase lasts respectively ~30 s and ~20 s (see Fig. 3), which is longer than the analytical predictions of $\tau_{\text{dyn}} \approx 11$ s and $\tau_{\text{dyn}} \approx 9.2$ s. This is due to the presence of a significant stage that is not included in the expression of $\tau_{\text{dyn}}$ where the growth of $B_{\phi}$ slows down before saturation.

The three characteristic timescales defined by Eqs. (49), (53), and (58) are plotted as a function of the fallback mass in Fig. 4 in addition to the characteristic timescale for the whole amplification process, which is defined as

$$ \tau_{\text{tot}} \equiv \tau_{\text{shear}} + \tau_{\text{TI}} + \tau_{\text{dyn}}. \qquad (59) $$

The vertical red dashed line represents the lower limit of fallback mass to form typical magnetars, which is estimated in Sect. 3.2 (corresponding to an asymptotic rotation period $P_{\text{fin}} \le 30$ ms). In the regime relevant to magnetar formation, the analytical and numerical estimates of the duration of the whole amplification process are in reasonable agreement, namely $\tau_{\text{tot}} \le 30$ s. In this regime, the phase which takes more time is the development of the Tayler instability. For $P_{\text{fin}} \ge 30$ ms, the comparison between the time at which $B_{\phi}$ saturates and $\tau_{\text{tot}}$ shows a significant difference, which is the consequence of the discrepancy noted above between the analytical estimate and numerical solution for $\tau_{\text{dyn}}$.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Different characteristic timescales as a function of the accreted mass: winding (black), Tayler instability (dark blue), and dynamo (red). The green line represents the sum of the three timescales. The shear rate is set at $q = 1$. The green crosses represent the entire amplification time obtained by integrating Eqs. (37)–(39). The red vertical line shows the lower limit of the accreted mass to form a magnetar with a radial field stronger than $B_Q \equiv m_e c^2 / e\hbar \approx 4.4 \times 10^{13}$ G using the predictions of Fuller et al. (2019).</figcaption></figure>

## 3.2. Magnetic field in the saturated regime

We now focus on the maximum magnetic field obtained at the end of the amplification phase. In the following discussion, this saturated magnetic field will be considered as a proxy for the magnetar’s magnetic field and its ‘radial’ component will be considered a proxy for the ‘dipolar’ component of the magnetic field. A more precise prediction would require a description of the relaxation towards a stable equilibrium, which is left for future studies.

In the top panel of Figs. 2 and 3, we see that the saturation intensities are close to their associated horizontal dashed lines, which illustrate the predictions of Eqs. (37)–(39) for values of the shear rate $q$ and the angular rotation frequency $\Omega$ reached at the time of magnetic field saturation. Therefore, these equations can be used to estimate the intensity of the saturated magnetic field. However, the angular frequency at $\tau_{\text{tot}}$ is still lower than its asymptotic value represented by the violet dashed line in the bottom panel. We estimate $\Omega(t = \tau_{\text{tot}})$ analytically by integrating Eq. (42)

$$ \Omega(t = \tau_{\text{tot}}) = \Omega_{\text{fin}} - \left( \frac{t_0}{\tau_{\text{tot}} + t_0} \right)^{2/3} (\Omega_{\text{fin}} - \Omega_{\text{init}}). \qquad (60) $$

Assuming that the timescales for the dynamo are roughly the same for the two models, that is, that of Fuller et al. (2019) and that of Spruit (2002), we also evaluate the expressions of the saturated magnetic field derived by Spruit (2002; Eqs. (18) and (19)) at $\Omega(t = \tau_{\text{tot}})$.

We see in Fig. 5 that our analytical estimates of the saturated fields (solid lines) are close to the numerical values at the peak of the solutions of Eqs. (44)–(46) (plus symbols). The small difference that appears for shorter rotation periods is due to the angular momentum transport, as discussed in see Sect. 4.1.

Using the maximum magnetic field as a proxy for the magnetar’s magnetic field may lead to an overestimation because a fraction of the magnetic energy can be dissipated during the relaxation to a stable magnetic configuration. Although our
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Predicted intensities for the saturated components of the magnetic field as a function of the accreted mass using the formalisms of Fuller et al. (2019; Eqs. (37)–(39); solid lines) and Spruit (2002; Eqs. (18)–(19); dash-dot lines). The shear rate is set at $q = 1$. These intensities are compared to the magnetic field reached at maximum intensity (blue and green plus signs) and at $q = 0$ (blue and green cross signs) by integrating Eqs. (42)–(46) for several fallback masses. Grey areas represent the estimated range of the dipolar magnetic field strength from regularly observed magnetars (dark grey) and from the three detected low-field Galactic magnetars (light grey). The vertical lines show the lower limit on the fallback mass (upper limit on the rotation period) needed to form a magnetar with a radial field stronger than $B_Q \approx 4.4 \times 10^{13}$ G for the prediction of Fuller et al. (2019; dashed red) and Spruit (2002; dotted red).</figcaption></figure>

model is unable to describe this relaxation process, we can get an idea of the robustness of our proxy by comparing to the magnetic field intensity at the time when $q = 0$. Figure 5 shows that this other proxy (cross signs) is smaller by between ~10% and ~50% but follows the same trends as the maximum magnetic field. Such a moderate difference would not change our main conclusions and suggests that the maximum magnetic field is a meaningful proxy for the final magnetic field.

In Fig. 5, the observed range of dipolar magnetic field for magnetars is fixed between the quantum electron critical field $B_Q \equiv m_e c^2 / e\hbar \approx 4.4 \times 10^{13}$ G and $B_{dip} \sim 2 \times 10^{15}$ G, the dipole field of the ‘most magnetised’ magnetar SGR 1806-20 (Olausen & Kaspi 2014). We find that the radial magnetic fields $B_{r,F}^{\text{sat}}$ and $B_{r,S}^{\text{sat}}$ fall in this range for accreted masses $M_{\text{acc}} \gtrsim 1.1 \times 10^{-2} M_{\odot}$ and $M_{\text{acc}} \gtrsim 4 \times 10^{-2} M_{\odot}$ (i.e. asymptotic rotation periods $P_{\text{fin}} \lesssim 28$ ms and $P_{\text{fin}} \lesssim 8$ ms), respectively. This confirms that magnetar-like magnetic fields can be formed over a wide range of accreted masses. The analytical predictions also show that, in the regime relevant for magnetar formation, the azimuthal component $B_{\phi} \gtrsim 4 \times 10^{15}$ G for both saturation mechanisms, which is significantly stronger than the radial component.

For lower accreted masses spinning-up the PNS to periods ranging from 28 to 64 ms (between 8 and 14 ms for the predictions of Spruit 2002), our scenario may produce radial magnetic fields $B_r$ as strong as the dipolar fields diagnosed in low-field magnetars (Rea et al. 2010, 2012, 2014). Moreover, the strength of the associated azimuthal field is $B_{\phi} \sim (1-3) \times 10^{15}$ G, which can be related to the non-dipolar magnetic field needed to produce the outbursts and chaotic bursts observed in magnetars (Thompson & Duncan 1995). This azimuthal magnetic field may also be the source of the proton cyclotron absorption lines observed in two low-field magnetars by Tiengo et al. (2013) and Rodríguez Castillo et al. (2016). Thus, our model provides a possible explanation of low-field magnetar formation. It is an alternative to the initial interpretation proposed by Rea et al. (2010), which invokes $\ge 1$ Myr ‘old’ (or ‘worn-out’) magnetars whose initial strong dipolar field of ~few $\times 10^{14}$ G has decayed due to Ohmic and Hall processes. This diffusion could be enhanced by the presence of a strong initial toroidal field $\ge 10^{16}$ G (Turolla et al. 2011).

As Eqs. (37)–(39) give orders of magnitude, Fuller et al. (2019) parameterised them with a prefactor denoted $\alpha$. We used $\alpha \sim 1$ as obtained by Fuller et al. (2019) for evolved stars in the subgiant and red giant branches by calibrating $\alpha$ on asteroseismic measurements. However, Eggenberger et al. (2019) find $\alpha \sim 0.5$ for subgiant stars on the one hand, and $\alpha \sim 1.5$ for red giant stars on the other. Also, Fuller & Lu (2022) argue that $\alpha \sim 0.25$ if intermittent dynamo action is considered in radiative zones with insufficient shear to trigger a sustained dynamo. This smaller prefactor would imply a larger limit of accreted mass of $\sim 2 \times 10^{-2} M_{\odot}$ (i.e. a rotation period of $\sim 15$ ms).

# 4. Discussion

## 4.1. Angular momentum transport

In the previous section, we focused on the magnetic field amplification and did not discuss the angular momentum transport due to the Tayler-Spruit dynamo. Our analytical estimate of the saturated magnetic field is based on the assumption that the differential rotation is not erased before the end of the amplification. Indeed, Figs. 2 and 3 show that the angular momentum transport due to Maxwell stresses starts to be significant around the time of magnetic field saturation and that most of the angular momentum transport occurs afterwards. This can be explained by comparing the characteristic timescales of the dynamo loop phase $\tau_{\text{dyn}}$ (Eq. (58)) with those of angular momentum transport. These can be estimated at saturation using the expression of $B_{\phi,F}^{\text{sat}}$ (Eq. (37)) as

$$ \tau_{\text{AM}} \equiv \gamma_{\text{AM}}^{-1} = \frac{r^2}{v_{\text{AM}}} \sim \left(\frac{N}{\Omega}\right)^2 \Omega^{-1}, \qquad (61) $$

and

$$ \tau_{\text{dyn}} = \left( \frac{N}{q\omega} \right)^{4/3} \Omega^{-1}. \qquad (62) $$

The ratio of these two timescales is

$$ \frac{\tau_{\text{dyn}}}{\tau_{\text{AM}}} \sim \left(\frac{\Omega}{q^2 N}\right)^{2/3} \sim 0.3 \left(\frac{P}{10 \text{ ms}}\right)^{-2/3}. \qquad (63) $$

For all values of the accreted mass and corresponding rotation period considered in this paper, the angular momentum transport is therefore longer than the dynamo timescale. This explains why most of the angular momentum takes place after the dynamo saturation and justifies a posteriori our analytical estimate of the saturated magnetic field. At fast rotation periods of a few milliseconds, the two timescales are nonetheless close to each other; as a consequence, the angular momentum transport before saturation is not negligible, which explains the moderate discrepancy between our analytical estimate of the saturated magnetic field and the numerical results for short rotation periods (Fig. 5).
On the other hand, angular momentum transport due to the neutrino viscosity can be neglected because its typical timescale is much longer than the evolution timescales considered:

$$ \tau_n \equiv \frac{r^2}{\nu_n} \ge 3 \times 10^4 \text{ s}, \qquad (64) $$

where the neutrino viscosity $\nu_n$ is estimated with the approximate analytical expression of Keil et al. (1996) and Guilet & Müller (2015):

$$ \nu_n \sim 3 \times 10^8 \left( \frac{\rho}{10^{14} \text{ g cm}^{-3}} \right)^{-2} \left( \frac{T}{5 \text{ MeV}} \right)^2 \text{cm}^2 \text{s}^{-1}. \qquad (65) $$

## 4.2. Neutrinos

We demonstrate above that angular momentum transport by either the magnetic field or the neutrino viscosity does not significantly impact our results. However, the neutrino flux coming from the accretion is also expected to extract a fraction of the angular momentum of the PNS (Janka et al. 2004; Bollig et al. 2021). To investigate whether it does not jeopardise the model, we use the following reasoning. As fallback is assumed to start several seconds after bounce in our model, one may assume that most of the angular momentum extraction by neutrino emission at these late times originate from fallback accretion rather than PNS cooling. Most of the fallback mass is likely to have a specific angular momentum $j_0$ which exceeds the Keplerian value at the PNS surface. It will therefore assemble into an accretion disk, settling at a radius $r_k$ where $j_{\text{kep}}(r_k) = j_0$. Its gravitational binding energy $E_{\text{bind}}$ will be at most all converted into neutrino radiation, that is per baryon with a rest mass $m_B$:

$$ \Delta E_v \le E_{\text{bind}}(r_k) \sim \frac{GM_{\text{PNS}}m_B}{r_k}, \qquad (66) $$

where we assume that the disk mass is small compared to the PNS mass. The corresponding specific angular momentum loss is

$$ \Delta j_v \le \frac{\Delta E_v}{m_B c^2} j_{\text{kep}}(r_k) \sim \frac{R_S}{2r_k} j_{\text{kep}}(r_k), \qquad (67) $$

where $R_S \equiv 2GM_{\text{PNS}}/c^2$ is the PNS Schwarzchild radius. $\Delta j_v$ is maximal at the PNS surface (i.e. when $r_k = r$), which implies

$$ \Delta j_v \le 0.185 j_{\text{kep}}(r), \qquad (68) $$

for the same parameters of a typical PNS introduced in Sect. 2.5. Therefore, the extraction of angular momentum by neutrino radiation is very inefficient.

## 4.3. Impact of the viscosity on the Tayler instability

In the reasoning developed above, we do not take into account the effects due to viscous processes, which might be important because they could be much larger than the effects of the resistivity in PNSs. Therefore, here we aim to address the question of their impact on the development of the Tayler instability. To the best of our knowledge, no analytical study of the Tayler instability has included the impact of viscosity. Hence, we use an approximate reasoning similar to that of Spruit (2002), which is based on a comparison of the instability growth timescale with the viscous damping timescale. This provides the following instability criterion:

$$ \sigma_{\text{TI}}^{-1} \sim \frac{\Omega}{\omega_A^2} \le \frac{l_r}{\nu}, \qquad (69) $$

with $\nu$ being the kinematic viscosity. Using the constraint on the radial length scale $l_r$ due to the stratification (Eq. (14)), we infer an instability criterion on the azimuthal magnetic field $B_\phi$ as a function of the viscosity:

$$ B_{\phi} > B_{\phi,c} \sim \Omega \left(\frac{N}{\Omega}\right)^{1/2} \left(\frac{\nu}{r^2 \Omega}\right)^{1/4}. \qquad (70) $$

This equation is similar to Eq. (12) but with the magnetic diffusivity substituted by the viscosity. To obtain an order of magnitude of this critical value, we must determine a value of the viscosity which is relevant for our scenario. As the fallback accretion occurs seconds to minutes after the PNS formation, the PNS has cooled down to temperatures $\le 1.1 \times 10^{11}$ K in the core and $\le 5 \times 10^{10}$ K in the outer region ($\le 10$ MeV and $\le 5$ MeV, respectively; Hüdepohl 2014). The neutrino mean free path can be approximated by (Thompson & Duncan 1993; Eq. (11)):

$$ l_n \sim 4 \times 10^4 \left( \frac{\rho}{10^{14} \text{ g cm}^{-3}} \right)^{-1/3} \left( \frac{T}{5 \text{ MeV}} \right)^{-3} \left( \frac{f(Y_p)}{1} \right) \text{cm}, \qquad (71) $$

where $f(Y_p)$ is function of the proton fraction close to unity. This length is larger than the maximum radial length scale (Eq. (14)):

$$ l_r \sim 4 \times 10^3 \left( \frac{B_\phi}{10^{15} \text{ G}} \right) \text{cm}. \qquad (72) $$

Therefore, neutrinos do not provide any relevant viscosity at the Tayler instability length scales and we must consider instead a microscopic viscosity such as the shear viscosity due to neutron-neutron scattering (Cutler & Lindblom 1987, Eq. (14)):

$$ v_s \sim 3 \times 10^{-2} \left( \frac{\rho}{10^{14} \text{ g cm}^{-3}} \right)^{5/4} \left( \frac{T}{5 \text{ MeV}} \right)^{-2} \text{cm}^2 \text{s}^{-1}. \qquad (73) $$

The associated critical magnetic field is therefore

$$ B_{\phi,c} \approx 10^{13} \left( \frac{v_s}{3 \times 10^{-2} \text{cm}^2 \text{s}^{-1}} \right)^{1/4} \left( \frac{\Omega}{200\pi \text{rad s}^{-1}} \right)^{1/4} \text{G}, \qquad (74) $$

which is four times stronger than the critical magnetic field inferred from the criterion of Spruit (2002) (Eq. (12)). However, this new critical magnetic field is still much weaker than the characteristic azimuthal magnetic field separating the winding phase from the phase in which the Tayler instability develops (when the growth rate of the Tayler instability reaches the winding rate; Eq. (48))

$$ B_{\phi} = \omega_{\text{A,TI}} \sqrt{4\pi\rho r^2} \sim 1.3 \times 10^{15} \left( \frac{\Omega}{200\pi \text{rad s}^{-1}} \right)^{2/3} \left( \frac{B_r}{10^{12} \text{G}} \right)^{1/3}. \qquad (75) $$

As a consequence, the viscosity is not expected to prevent the Tayler instability from growing and should not have a significant impact on our results. However, we note that our argument is approximate and would need to be upgraded through a linear analysis of the Tayler instability including the viscous processes.

## 4.4. Superfluidity and superconductivity

A last potential obstacle for our model may emerge from the crust formation and the superfluidity and superconductivity in the core, which occur during the cooling of the PNS. The outer crust is expected to start freezing a few minutes after the PNS
formation and the inner crust forms far later, between 1 and 100 yr after formation (Aguilera et al. 2008). Therefore, no part of the crust is formed during the time interval involved in our scenario.

The potential early appearance of superfluid neutrons or even superconductive protons in the PNS core at temperatures below $10^8-10^{10}$ K is worth discussing because the MHD approximation would not be sufficiently realistic and a multi-fluid approach would be more relevant (Glampedakis et al. 2011; Sinha & Sedrakian 2015). However, the 1D models of PNS cooling show higher temperatures than $10^{10}$ K in the PNS even after 15 s (e.g. Pons et al. 1999; Roberts 2012; Hüdepohl 2014; Roberts et al. 2017). Moreover, Gusakov & Kantor (2013) and Glampedakis & Jones (2014) brought forward a critical perturbed magnetic field strength above which superfluidity of neutrons dies out. Therefore, the MHD approximation is still valid for describing the PNS internal dynamics during the first 40 s following the core bounce.

## 5. Conclusions

In this paper, we propose a new scenario for magnetar formation, in which the Tayler-Spruit dynamo amplifies the large-scale magnetic field of a PNS spun up by SN fallback accretion. We develop a one-zone model describing the evolution of the magnetic field averaged over a PNS subject to fallback accretion. The equations describing the time evolution are solved numerically and compared successfully with analytical estimates of the final magnetic field and of the duration of each stage of the amplification process. Predictions for the different components of the magnetic field are therefore obtained as a function of the accreted mass for the two proposed saturation models of the Tayler-Spruit dynamo (Spruit 2002; Fuller et al. 2019). Our main conclusions can be summarised as follows:

- Radial magnetic fields spanning the full range of the magnetar dipole intensity can be formed for accreted masses compatible with the results of recent SN simulations. Our model predicts the formation of magnetar-like magnetic fields for accreted masses $M_{\text{acc}} \ge 1.1 \times 10^{-2} M_{\odot}$ for the saturation model of Fuller et al. (2019) and $M_{\text{acc}} \ge 4 \times 10^{-2} M_{\odot}$ for the saturation model of Spruit (2002). This corresponds to neutron star final rotation periods $P_{\text{fin}} \le 28$ ms and $P_{\text{fin}} \le 8$ ms, respectively.

- The azimuthal component of the magnetic field is predicted to be in the range $10^{15}-10^{16}$ G, which is stronger than the radial component by a factor of 10–100.

- In the regime relevant for magnetar formation, the magnetic field amplification lasts between 15 and 30 s. On such a timescale, the MHD equations assumed in the description of the Tayler-Spruit dynamo are expected to be valid. Furthermore, we have not identified any other process capable of interfering with the Tayler-Spruit dynamo by transporting angular momentum on a comparable or shorter timescale.

Our results therefore predict that magnetars can indeed be formed in our new scenario. Magnetar formation is possible at sufficiently long rotation periods to be compatible with the lower limit of 5 ms inferred from regular SN remnants associated with magnetars. With the saturation model of Fuller et al. (2019), the full range of magnetar fields can be obtained within this constraint, even those that exhibit a strong dipolar magnetic field $B_{\text{dip}} \sim 10^{15}$ G. On the other hand, with the saturation model proposed by Spruit (2002), only the lower end of the magnetar fields can be obtained with $P_{\text{fin}} < 5$ ms, while dipolar magnetic fields $\ge 2 \times 10^{14}$ G need faster rotation periods.

An important prediction of our scenario is the very intense toroidal magnetic field, which lies between $3 \times 10^{15}$ and $3 \times 10^{16}$ G for parameters corresponding to radial magnetic fields in the magnetar range. These values are compatible with the interpretation of the X-ray flux modulations observed in three magnetars as free precession driven by an intense toroidal magnetic field (Makishima et al. 2014, 2016, 2019, 2021).

The intense toroidal magnetic field predicted in our scenario also provides interesting perspectives from which to explain the formation of low-field magnetars. For radial magnetic fields in the range of the dipolar magnetic field deduced for these objects (Rea et al. 2010, 2012, 2013, 2014), our model predicts a toroidal magnetic field intensity of $\sim 1-3 \times 10^{15}$ G. Such non-dipolar magnetic fields are strong enough to be the energy source of the magnetar-like emission from these objects and to explain the variable absorption lines interpreted as proton cyclotron lines (Tiendo et al. 2013; Rodríguez Castillo et al. 2016). We therefore suggest that some of the low-field magnetars may be born with low dipolar magnetic fields, rather than evolve to this state as assumed in the ‘worn-out’ magnetar scenario.

A question arising from our study is the location of the magnetic field in the PNS, which cannot be captured by our one-zone model. As the shear due to fallback accretion is expected to be strongest in the outer region of the PNS, one may expect the magnetic field to be preferentially located in these outer layers. Such a concentration of the magnetic field near the surface would have interesting consequences for its long-term evolution because the magnetic field may be confined in the crust without significant magnetic flux threading the superconductive core. The long-term evolution of such a crust-confined magnetic field configuration has been thoroughly investigated by numerical simulations (e.g. Viganò 2013; Gourgouliatos et al. 2016; Pons & Viganò 2019). By contrast, if the magnetic field is also present in deeper regions, its evolution in the superconductive core and the transition layer with the crust must be taken into account. A few papers studied this evolution in numerical simulations (e.g. Henriksson & Wasserman 2013; Lander 2013; Ciolfi & Rezzolla 2013) but these lead to a slower magnetic field evolution that is incompatible with magnetar observations (Elfritz et al. 2016). These results would favour initial crust-confined magnetic fields but need to be confirmed in more realistic 3D simulations of the magneto-thermal evolution in the whole neutron star. We also note that the localisation of the magnetic field in our magnetar formation scenario should be studied in more detail. On the one hand, the stratification increases by a factor of $\sim 2.5$ close to the PNS surface, which might weaken the magnetic field and confine it closer to the surface. On the other hand, the shear can be expected to become significant in the bulk of the PNS after angular momentum has been partly redistributed by the Tayler-Spruit dynamo. Some of the magnetic field may also be transported to deeper regions via the Tayler instability or during the relaxation to a stable equilibrium.

Another relevant question is the geometry of the magnetic field amplified by the Tayler-Spruit dynamo. One should keep in mind that our comparison to magnetars relies on the assumption that the generated radial field $B_r$ is mostly dipolar. Although the real geometry of the poloidal magnetic field generated by the Tayler-Spruit dynamo is not known, it is likely to be partly non-dipolar, meaning that the large-scale dipolar magnetic field is a fraction of the radial field $B_r$. Therefore, corresponding predictions should be refined by studying dedicated multi-dimensional models. In Petitdemange et al. (2022), a dynamo similar to the Tayler-Spruit dynamo has been found through numerical simulations in a configuration where the
surface rotates slower than the core, which is therefore differ-
ent to the case of spun-up PNS. Moreover, the observed magne-
tars are cooled-down neutron stars with a stable configuration of
magnetic field. Hence, the study of the magnetic field relaxation
from a turbulent saturated state to a stable configuration is impor-
tant to estimate a more realistic intensity of the dipolar poloidal
field. Thus, numerical simulations will be essential to further
study of the evolution of the magnetic field geometry in our
framework.

A salient feature of our fallback scenario is that it decou-
ples magnetar formation from rapid progenitor rotation and from
strong magnetisation of the pre-collapse stars. Rapid progeni-
tor rotation is necessary for magnetar formation by the convec-
tive dynamo, which requires initial NS spin periods of ≤10 ms
(Raynaud et al. 2022), and by the magnetorotational instability
(Reboul-Salze et al. 2021a,b). Strong magnetisation of the pre-
collapse star on the other hand is a crucial aspect in the fos-
sil field scenario or the stellar merger scenario (Schneider et al.
2019). Instead of requiring fast rotation or strong magnetic
field in the progenitor core, our scenario predicts magnetar for-
mation when fallback deposits a sufficient amount of angular
momentum on the PNS surface. With the angular momentum
of the mass accreted by the NS being limited by the Keplerian
value, magnetars are formed for accreted masses of more than
~1.1 × 10⁻² M☉ (case of Fuller et al. 2019) and ~4 × 10⁻² M☉
(case of Spruit 2002) in our scenario. The fallback mass should
be several times larger than the accreted mass, because angular
momentum loss must be expected to lead to mass loss during
the accretion process. Therefore, fallback masses of more than a
few 10⁻² M☉ to 10⁻¹ M☉ seem to be needed. Based on 1D mod-
els of neutrino-driven core-collapse SN explosions, this indi-
cates a preference for single stars with zero-age-main-sequence
(ZAMS) masses above about ~18 M☉ (Sukhbold et al. 2016) and
helium stars (hydrogen-stripped stars in binaries) with ZAMS
masses above 30–40 M☉ (depending on details of the mass-loss
evolution); although the compactness differences between the
single-star models of Sukhbold & Woosley (2014) compared to
those of Sukhbold et al. (2018) as well as 3D explosion effects
(which increase the fallback mass; Janka et al. 2022) may shift
these ZAMS masses to lower values. This would be consis-
tent with the observations constraining magnetar progenitors
to masses higher than 30 M☉ (Gaensler et al. 2005; Bibby et al.
2008; Clark et al. 2008) and also with the case of the magne-
tar SGR1900+14, whose progenitor mass was estimated to be
17 ± 2 M☉ (Davies et al. 2009).

While our model avoids the uncertainty of the progenitor
core rotation and magnetic field, it implies coping with the
uncertainties on the fallback process. A precise modelling of the
fallback depends on such challenging questions as how a long-
lasting post-explosion phase where downflows to the PNS coex-
ist with outflows of neutrino-heated matter transitions into the
fallback accretion as discussed by Janka et al. (2022); the com-
plex dynamical processes that determine the fraction of fallback
matter that gets accreted by the PNS from a fallback disk; and
the efficiency of the accretion to spin-up the PNS. Our scenario
should therefore be explored in more depth by more realistic fall-
back models.

Following its saturation, the PNS magnetic field may inter-
act with the newly formed disk of fallback matter and is strong
enough to influence the fallback accretion mechanism. We did
not model this interaction because our study was focused on
the phase of magnetic field amplification. Nevertheless, this
could strongly influence the rotation of the newly born magne-
tar. The evolution of the PNS-fallback disk system depends on

three characteristic radii (Metzger et al. 2018; Beniamini et al. 2020; Lin et al. 2020; Ronchi et al. 2022): (i) the magnetospheric radius $r_m$, which is the radius at which the matter is blocked by the magnetic barrier, (ii) the corotation radius $r_c$ where the matter has the same rotation frequency as the PNS, and (iii) the light cylinder radius $r_{lc}$, which is the ratio of light speed to the PNS rotation frequency. The strong magnetic field repels the magnetosphere behind the corotation radius (i.e. $r_c < r_m$) which stops the accretion and so the PNS spin-up. If the fallback accretion rate of the disk is large enough, the inner part of the disk penetrates the light cylinder (i.e. $r_{lc} > r_m$) and opens up a part of the magnetic field lines. The PNS-fallback disk system enters the so-called propeller regime and the PNS angular momentum is transported towards the disk via the magnetic dipole torque. This mechanism is thought to extract the PNS angular momentum very efficiently; for instance Beniamini et al. (2020) even predict magnetars spun down to rotation periods of ~$10^6$ s after ~$10^3$ yr. For this reason, this scenario is often invoked to explain the ultra-long-period magnetars such as 1E 1613 (e.g. De Luca et al. 2006; Li 2007; Rea et al. 2016) or the recently observed GLEAM-X J162759.5-523504.3 (Ronchi et al. 2022), which have respective rotation periods of ~$2.4 \times 10^4$ s and ~$1.1 \times 10^3$ s. It would be interesting to include such a spin-down model in our magnetar formation scenario in order to obtain a prediction of the rotation period at later times.

Finally, the PNS-fallback disk system has also been invoked
to explain the light curve of luminous and extreme SNe of types
Ib/c (e.g. Dexter & Kasen 2013; Metzger et al. 2018; Lin et al.
2021). We may also expect our scenario to produce these
types of explosions depending on the amount of accreted mass
during the dynamo process. First, PNSs that have accreted
~2–3 × 10⁻² M☉ of fallback matter before the magnetic field
saturation have rotation periods of around 10–20 ms, which are
too slow to produce extreme explosions. According to our sce-
nario, their typical magnetic field is of 1–5 × 10¹⁴ G, which
would lead to regular luminous SNe Ib/c. Their light curve
would be dominated by the PNS spin-down luminosity instead
of the ⁵⁶Ni decay luminosity (Ertl et al. 2020; Afsariardchi et
2021). Second, for fallback masses spinning up the PNSs
to millisecond rotation periods, the magnetic field saturates
a few 10 s after the core bounce at $B_r \ge 5 \times 10^{14}$ G. The
rotational energy can be kept for later times and be slowly
extracted to irradiate its environment, which might lead to
superluminous SNe I (Woosley 2010; Kasen & Bildsten 2010;
Bersten et al. 2016; Margalit et al. 2018; Lin et al. 2020, 2021).
Finally, to produce extreme explosions such as hypernovae,
which have approximately ten times larger kinetic energies and
much higher ⁵⁶Ni yields than the vast majority of CCSNe,
an energy injection within a timescale of ≤1 s is required
(Barnes et al. 2018) to explain the large masses of ⁵⁶Ni ≥
0.2 M☉ inferred from their light curves (e.g. Woosley & Bloom
2006; Drout et al. 2011; Nomoto et al. 2011). For rotation peri-
ods of ≤1 ms, which correspond to rotational energies of
≥3 × 10⁵² erg, our model provides a radial magnetic field of $B_r \ge
2.6 \times 10^{16}$ G, which may be enough to inject the energy quickly
through magnetic dipole spin-down only (Suwa & Tominaga
2015). The presence of a propeller regime would enhance the
PNS spin-down such that a weaker dipolar magnetic field of
~2 × 10¹⁵ G would also produce a hypernova (Metzger et al.
2018) but through a propeller-powered explosion. Therefore,
our magnetic-field-amplification scenario by PNS accretion
or fallback accretion may be of relevance to a wide vari-
ety of magnetar-powered phenomena in different types of SN
events.
**Acknowledgements.** We thank the referee Jim Fuller for his thorough reading and relevant comments, which helped improve this paper. This research was supported by the European Research Council through the ERC starting grant MagBURST No. 715368. Hans-Thomas Janka acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through Sonderforschungsbereich (Collaborative Research Centre) SFB-1258 'Neutrinos and Dark Matter in Astro- and Particle Physics (NDM)' and under Germany's Excellence Strategy through Cluster of Excellence ORIGINS (EXC-2094)-390783311. We are grateful to Thierry Foglizzo and Ludovic Petitdemange for useful discussions.

## References

Afsariardchi, N., Drout, M. R., Khatami, D. K., et al. 2021, ApJ, 918, 89

Aguilera, D. N., Pons, J. A., & Miralles, J. A. 2008, A&A, 486, 255

Barnes, J., Duffell, P. C., Liu, Y., et al. 2018, ApJ, 860, 38

Beniamini, P., Hotokezaka, K., van der Horst, A., & Kouveliotou, C. 2019, MNRAS, 487, 1426

Beniamini, P., Wadiasingh, Z., & Metzger, B. D. 2020, MNRAS, 496, 3390

Bersten, M. C., Benvenuto, O. G., Orellana, M., & Nomoto, K. 2016, ApJ, 817, L8

Bibby, J. L., Crowther, P. A., Furness, J. P., & Clark, J. S. 2008, MNRAS, 386, L23

Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59

Bollig, R., Yadav, N., Kresse, D., et al. 2021, ApJ, 915, 28

Bonanno, A., & Guarrieri, F. 2017, Astron. Nachr., 338, 516

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

Braithwaite, J., & Spruit, H. C. 2004, Nature, 431, 819

Braithwaite, J., & Spruit, H. C. 2017, R. Soc. Open Sci., 4, 160271

Bugli, M., Guilet, J., Obergaulinger, M., Cerdá-Durán, P., & Aloy, M. A. 2020, MNRAS, 492, 58

Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443

Burrows, A., Dessart, L., Livne, E., Ott, C. D., & Murphy, J. 2007, ApJ, 664, 416

Chan, C., Müller, B., & Heger, A. 2020, MNRAS, 495, 3751

Chevalier, R. A. 1989, ApJ, 346, 847

CHIME/FRB Collaboration, Andersen, B. C., Bandura, K. M., et al. 2020, Nature, 587, 54

Ciolfi, D., & Rezzolla, L. 2013, MNRAS, 435, L43

Clark, J. S., Muno, M. P., Negueruela, I., et al. 2008, A&A, 477, 147

Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961

Coti Zelati, F., Borghese, A., Israel, G. L., et al. 2021, ApJ, 907, L34

Cutler, C., & Lindblom, L. 1987, ApJ, 314, 234

Davies, B., Figer, D. F., Kudritzki, R.-P., et al. 2009, ApJ, 707, 844

De Luca, A., Caraveo, P. A., Mereghetti, S., Tiengo, A., & Bignami, G. F. 2006, Science, 313, 814

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

Dessart, L., Burrows, A., Livne, E., & Ott, C. D. 2008, ApJ, 673, L43

Dessart, L., Hillier, D. J., Waldman, R., Livne, E., & Blondin, S. 2012, MNRAS, 426, L76

Dexter, J., & Kasen, D. 2013, ApJ, 772, 30

Drout, M. R., Soderberg, A. M., Gal-Yam, A., et al. 2011, ApJ, 741, 97

Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9

Eggenberger, P., den Hartogh, J. W., Buldgen, G., et al. 2019, A&A, 631, L6

Elfritz, J. G., Pons, J. A., Rea, N., Glampedakis, K., & Viganò, D. 2016, MNRAS, 456, 4461

Ertl, T., Woosley, S. E., Sukhbold, T., & Janka, H. T. 2020, ApJ, 890, 51

Evans, W. D., Klebesadel, R. W., Laros, J. G., et al. 1980, ApJ, 237, L7

Ferrario, L., & Wickramasinghe, D. 2006, MNRAS, 367, 1323

Fuller, J., & Lu, W. 2022, MNRAS, 511, 3951

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

Gabler, M., Cerdá-Durán, P., Stergioulas, N., Font, J. A., & Müller, E. 2018, MNRAS, 476, 4199

Gaensler, B. M., McClure-Griffiths, N. M., Oey, M. S., et al. 2005, ApJ, 620, L95

Gill, R., & Heyl, J. 2007, MNRAS, 381, 52

Glampedakis, K., & Jones, D. I. 2014, MNRAS, 439, 1522

Glampedakis, K., Andersson, N., & Samuelsson, L. 2011, MNRAS, 410, 805

Gompertz, B., & Fruchter, A. 2017, ApJ, 839, 49

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

Gotz, D., Israel, G. L., Mereghetti, S., et al. 2006, ATel, 953, 1

Gourgouliatos, K. N., Wood, T. S., & Hollerbach, R. 2016, Proc. Natl. Acad. Sci., 113, 3944

Guerrero, G., Del Sordo, F., Bonanno, A., & Smolarkiewicz, P. K. 2019, MNRAS, 490, 4281

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

Guilet, J., Reboul-Salze, A., Raynaud, R., Bugli, M., & Gallet, B. 2022, MNRAS, 516, 4346

Gusakov, M. E., & Kantor, E. M. 2013, MNRAS, 428, L26

Henriksson, K. T., & Wasserman, I. 2013, MNRAS, 431, 2986

Hu, R.-Y., & Lou, Y.-Q. 2009, MNRAS, 396, 878

Hüdepohl, L. 2014, Ph.D. Thesis, Technical University of Munich, Germany

Hurley, K., Cline, T., Mazets, E., et al. 1999, Nature, 397, 41

Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098

Inserra, C., Smartt, S. J., Jerkstrand, A., et al. 2013, ApJ, 770, 128

Israel, G. L., Belloni, T., Stella, L., et al. 2005, ApJ, 628, L53

Janka, H. T. 2004, in IAU Symp., 218, 3

Janka, H.-T., Wongwathanarat, A., & Kramer, M. 2022, ApJ, 926, 9

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

Kasen, D., & Bildsten, L. 2010, ApJ, 717, 245

Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261

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

Kouveliotou, C. 1999, Proc. Natl. Acad. Sci., 96, 5351

Kouveliotou, C., Fishman, G. J., Meegan, C. A., et al. 1994, Nature, 368, 125

Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, ApJ, 896, 102

Lander, S. K. 2013, Phys. Rev. Lett., 110, 071101

Lander, S. K. 2021, MNRAS, 507, L36

Lattimer, J. M., & Schutz, B. F. 2005, ApJ, 629, 979

Lattimer, J. M., & Swesty, D. F. 1991, Nucl. Phys. A, 535, 331

Li, X.-D. 2007, ApJ, 666, L81

Lin, W. L., Wang, X. F., Wang, L. J., & Dai, Z. G. 2020, ApJ, 903, L24

Lin, W., Wang, X., Wang, L., & Dai, Z. 2021, ApJ, 914, L2

Lü, H.-J., & Zhang, B. 2014, ApJ, 785, 74

Ma, L., & Fuller, J. 2019, MNRAS, 488, 4338

Makishima, K., Enoto, T., Hiraga, J. S., et al. 2014, Phys. Rev. Lett., 112, 171102

Makishima, K., Enoto, T., Murakami, H., et al. 2016, PASJ, 68, S12

Makishima, K., Murakami, H., Enoto, T., & Nakazawa, K. 2019, PASJ, 71, 15

Makishima, K., Enoto, T., Yoneda, H., & Odaka, H. 2021, MNRAS, 502, 2266

Margalit, B., Metzger, B. D., Thompson, T. A., Nicholl, M., & Sukhbold, T. 2018, MNRAS, 475, 2659

Martin, J., Rea, N., Torres, D. F., & Papitto, A. 2014, MNRAS, 444, 2910

Masada, Y., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 75

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

Metzger, B. D., Beniamini, P., & Giannios, D. 2018, ApJ, 857, 95

Mösta, P., Richers, S., Ott, C. D., et al. 2014, ApJ, 785, L29

Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2013, Nature, 502, 346

Nomoto, K., Maeda, K., Tanaka, M., & Suzuki, T. 2011, Ap&SS, 336, 129

Obergaulinger, M., & Aloy, M. Å. 2020, MNRAS, 492, 461

Obergaulinger, M., & Aloy, M. Å. 2021, MNRAS, 503, 4942

Obergaulinger, M., & Aloy, M. Å. 2022, MNRAS, 512, 2489

Obergaulinger, M., Cerdá-Durán, P., Müller, E., & Aloy, M. A. 2009, A&A, 498, 241

Olausen, S. A., & Kaspi, V. M. 2014, ApJS, 212, 6

Petit, V., Wade, G. A., Schneider, F. R. N., et al. 2019, MNRAS, 489, 5669

Petitdemange, L., Marcotte, F., & Gissinger, C. 2022, arXiv e-prints [arXiv:2206.13819]

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

Pons, J. A., & Viganò, D. 2019, Liv. Rev. Comput. Astrophys., 5, 3

Pons, J. A., Reddy, S., Prakash, M., Lattimer, J. M., & Miralles, J. A. 1999, ApJ, 513, 780

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

Raynaud, R., Cerdá-Durán, P., & Guilet, J. 2022, MNRAS, 509, 3410

Rea, N., Esposito, P., Turolla, R., et al. 2010, Science, 330, 944

Rea, N., Israel, G. L., Esposito, P., et al. 2012, Nature, 524, 27

Rea, N., Israel, G. L., Pons, J. A., et al. 2013, ApJ, 770, 65

Rea, N., Viganò, D., Israel, G. L., Pons, J. A., & Torres, D. F. 2014, ApJ, 781, L17

Rea, N., Borghese, A., Esposito, P., et al. 2016, ApJ, 828, L13

Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2021a, A&A, 645, A109

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

Roberts, L. F. 2012, Ph.D. Thesis, University of California, Santa Cruz, USA

Roberts, L. F., & Reddy, S. 2017, in Handbook of Supernovae, eds. A. W. Alsabti, & P. Murdin (Springer International Publishing), 1605

Roberts, O. J., Veres, P., Baring, M. G., et al. 2021, Nature, 589, 207

Rodríguez Castillo, G. A., Israel, G. L., Tiengo, A., et al. 2016, MNRAS, 456, 4145

Ronchi, M., Rea, N., Graber, V., & Hurley-Walker, N. 2022, ApJ, 934, 184

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

Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211
Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2020, MNRAS, 495, 2796

Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 1998a, Nucl. Phys. A, 637, 435

Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 1998b, Prog. Theor. Phys., 100, 1013

Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 2011, ApJS, 197, 20

Sinha, M., & Sedrakian, A. 2015, Phys. Rev. C, 91, 035805

Spruit, H. C. 1999, A&A, 349, 189

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

Stockinger, G., Janka, H. T., Kresse, D., et al. 2020, MNRAS, 496, 2039

Strohmayer, T. E., & Watts, A. L. 2005, ApJ, 632, L111

Sukhbold, T., & Woosley, S. E. 2014, ApJ, 783, 10

Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38

Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93

Suwa, Y., & Tominaga, N. 2015, MNRAS, 451, 282

Svinkin, D., Frederiks, D., Hurley, K., et al. 2021, Nature, 589, 211

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

Takiwaki, T., Kotake, K., & Sato, K. 2009, ApJ, 691, 1360

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

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

Thompson, C., & Duncan, R. C. 1995, MNRAS, 275, 255

Tiengo, A., Esposito, P., Mereghetti, S., et al. 2013, Nature, 500, 312

Turolla, R., Zane, S., Pons, J. A., Esposito, P., & Rea, N. 2011, ApJ, 740, 105

Viganò, D. 2013, Ph.D. Thesis, University of Alacant, Spain

Vink, J., & Kuiper, L. 2006, MNRAS, 370, L14

White, C. J., Burrows, A., Coleman, M. S. B., & Vartanyan, D. 2022, ApJ, 926, 111

Woosley, S. E. 2010, ApJ, 719, L204

Woosley, S. E., & Bloom, J. S. 2006, ARA&A, 44, 507

Zahn, J. P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145

Zhang, B., & Mészáros, P. 2001, ApJ, 552, L35

Zhou, P., Vink, J., Safi-Harb, S., & Miceli, M. 2019, A&A, 629, A51