11institutetext: Université Paris-Saclay, Université Paris Cité, CEA, CNRS,
AIM, 91191, Gif-sur-Yvette, France
Université Paris Cité, Université Paris-Saclay, CNRS, CEA, AIM, F-91191 Gif-
sur-Yvette, France
Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
D-14476 Potsdam, Germany

# Tayler-Spruit dynamo in stably stratified rotating fluids: Application to
proto-magnetars

P. Barrère<EMAIL_ADDRESS>11 J. Guilet 11 R. Raynaud 22 A. Reboul-Salze
33112233

(Received ?????; accepted ?????)

The formation of highly magnetized young neutron stars, called magnetars, is
still a strongly debated question. A promising scenario invokes the
amplification of the magnetic field by the Tayler-Spruit dynamo in a proto-
neutron star (PNS) spun up by fallback. Barrère et al. (2023) supports this
scenario by demonstrating that this dynamo can generate magnetar-like magnetic
fields in stably stratified Boussinesq models of a PNS interior. To further
investigate the Tayler-Spruit dynamo, we perform 3D-MHD numerical simulations
with the MagIC code varying the ratio between the Brunt-Väisälä frequency and
the rotation rate. We first demonstrate that a self-sustained dynamo process
can be maintained for a Brunt-Väisälä frequency about 4 times higher than the
angular rotation frequency. The generated magnetic fields and angular momentum
transport follow the analytical scaling laws of Fuller et al. (2019), which
confirms our previous results. We also report for the first time the existence
of an intermittent Tayler-Spruit dynamo. For a typical PNS Brunt-Väisälä
frequency of ${10}^{3}\text{\,}\mathrm{s}^{-1}$, the axisymmetric toroidal and
dipolar magnetic fields range between
$$1.2\text{\times}{10}^{15}\text{\,}$-$2\text{\times}{10}^{16}\text{\,}\mathrm{G}$$
and
$$1.4\text{\times}{10}^{13}\text{\,}$-$3\text{\times}{10}^{15}\text{\,}\mathrm{G}$$,
for rotation periods of $1-$10\text{\,}\mathrm{m}\mathrm{s}$$. Thus, our
results provide numerical evidence that our scenario can explain the formation
of magnetars. As the Tayler-Spruit dynamo is often invoked for the angular
momentum transport in stellar radiative zones, our results are also of
particular importance in this field and we provide a calibration of the Fuller
et al.’s prescription based on our simulations, with a dimensionless
normalisation factor of the order of ${10}^{-2}\text{\,}$.

###### Key Words.:

stars: magnetic fields – stars: magnetars – supernovae: general –
magnetohydrodynamics (MHD) – dynamo – methods: numerical

##  1 Introduction

Soft gamma repeaters and anomalous X-ray pulsars are two classes of neutron
stars (NSs) that exhibit a wide variety of high-energy emissions from short
chaotic bursts during outbursts phases (e.g. Gotz et al., 2006; Younes et al.,
2017; Coti Zelati et al., 2018, 2021) to giant flares (Evans et al., 1980;
Hurley et al., 1999, 2005; Svinkin et al., 2021), which are the brightest
events observed in the Milky Way. These neutron stars are called magnetars
because their emissions have been shown to be powered by the dissipation of
their ultra-strong magnetic fields (Kouveliotou et al., 1994). Indeed, these
emissions show that they rotate with periods of
$$2\text{\,}$-$12\text{\,}\mathrm{s}$$ and have stronger rotation braking than
typical NSs (e.g. Rea et al., 2012; Olausen & Kaspi, 2014). If the spin-down
is due to the extraction of rotational energy by a magnetic dipole, we can
infer that most magnetars exhibit a surface magnetic dipole of
$${10}^{14}\text{\,}$-${10}^{15}\text{\,}\mathrm{G}$$, which are the strongest
known in the Universe. Three magnetars, however, display weaker magnetic
dipoles of
$$6\text{\times}{10}^{12}\text{\,}$-$4\text{\times}{10}^{13}\text{\,}\mathrm{G}$$
(Rea et al., 2010, 2012, 2013, 2014) but absorption lines detected in the
X-ray spectra of two of these magnetars suggest the presence of stronger non-
dipolar magnetic fields of
$$2\text{\times}{10}^{14}\text{\,}$-$2\text{\times}{10}^{15}\text{\,}\mathrm{G}$$
(Tiengo et al., 2013; Rodríguez Castillo et al., 2016). These ‘low-field’
magnetars therefore demonstrate that an ultra-strong surface magnetic dipole
is not necessary for a neutron star to produce magnetar-like emission.

Magnetars are also suspected to be the central engine of extreme events. In
combination with a millisecond rotation period, magnetars in their proto-
neutron star (PNS) stage may power magnetorotational supernova (SN) explosions
which are more energetic than typical neutrino-driven SNe (e.g. Burrows et
al., 2007; Dessart et al., 2008; Takiwaki et al., 2009; Kuroda et al., 2020;
Bugli et al., 2020, 2021, 2023; Obergaulinger & Aloy, 2020, 2021, 2022). The
formation of a millisecond magnetar is a popular scenario to explain super-
luminous SNe (Woosley, 2010; Kasen & Bildsten, 2010; Dessart et al., 2012;
Inserra et al., 2013; Nicholl et al., 2013) and hypernovae, of which the
latter are associated to long gamma-ray bursts (GRBs ; 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).
In the case of binary NS mergers, the NS remnant may be a magnetar whose
magnetic fields power the plateau phase observed in short GRBs afterglows
(e.g. Lü & Zhang, 2014; Gompertz et al., 2014; Kiuchi et al., 2024).

Recently, the observation of the fast radio burst FRB 200428 was associated to
X-ray bursts of the magnetar SGR 1935+2154 (Bochenek et al., 2020; CHIME/FRB
Collaboration et al., 2020; Mereghetti et al., 2020; Zhu et al., 2023; Tsuzuki
et al., 2024), which supports magnetar-powered emission scenarios to explain
at least a fraction of FRBs.

In order to better understand these phenomena, it is thus essential to
investigate the question of magnetar formation and especially the origin of
their ultra-strong magnetic fields. The association of a few magnetars with SN
remnants suggests that they are born during core-collapse SNe (Vink & Kuiper,
2006; Martin et al., 2014; Zhou et al., 2019). The magnetic fields could be
amplified during the core-collapse due to the magnetic flux conservation
(Ferrario & Wickramasinghe, 2006; Hu & Lou, 2009; Schneider et al., 2019;
Shenar et al., 2023). However, the magnetic field of the iron core of the
progenitor star is not constrained by observations and it is uncertain whether
this scenario can explain the whole magnetar population (Makarenko et al.,
2021). A second type of scenario invokes a dynamo action in the newly formed
PNS to generate strong large-scale magnetic fields. Three mechanisms have been
studied: the convective dynamo (Thompson & Duncan, 1993; Raynaud et al., 2020,
2022; Masada et al., 2022; White et al., 2022), the magnetorotational
instability (MRI)-driven dynamo (Obergaulinger et al., 2009; Mösta et al.,
2014; Rembiasz et al., 2017; Reboul-Salze et al., 2021, 2022; Guilet et al.,
2022), and the Tayler-Spruit dynamo (Barrère et al., 2022, 2023).

The two former dynamos have been shown to form magnetar-like magnetic fields
in the case of fast rotation. In the framework of the millisecond magnetar
model, these mechanisms are therefore promising to explain the formation of
the central engine extreme explosions. Nevertheless, two uncertainties remain.
First, SN remnants associated with magnetars show typical explosion energy of
$\sim${10}^{51}\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}$$. This implies that
most magnetars were born in standard core-collapse supernovae (CCSNe), which
require slower rotation periods of at least $5\text{\,}\mathrm{m}\mathrm{s}$
(Vink & Kuiper, 2006). Second, the rotation is assumed to stem from a fast-
rotating progenitor core. It is unclear whether there is a large enough
fraction of these progenitors to explain the entire magnetar population.

To address these points, we developed in Barrère et al. (2022) a new formation
scenario in which the PNS rotation is determined by the fallback, which is the
matter that is initially ejected by the SN explosion before eventually falling
back onto the PNS. 3D CCSNe numerical simulations show that the fallback
accretion can significantly spin up the PNS surface (Chan et al., 2020;
Stockinger et al., 2020; Janka et al., 2022). We argued that the differential
rotation caused by this spin-up triggers the development of the Tayler-Spruit
dynamo. This dynamo mechanism is driven by the combination of the differential
rotation and the Tayler instability, which is an instability triggered by
perturbations of a purely toroidal magnetic field (Tayler, 1973; Goossens &
Tayler, 1980; Goossens et al., 1981). Studies in stellar evolution (such as
Eggenberger et al., 2005; Cantiello et al., 2014; Eggenberger et al., 2019a,
b; den Hartogh et al., 2020; Griffiths et al., 2022) often rely on this
mechanism to explain the strong angular momentum transport (AMT) inferred via
asteroseismology in sub-giant/red giant stars (e.g. Mosser et al., 2012;
Deheuvels et al., 2014, 2015; Gehan et al., 2018). As our model in Barrère et
al. (2022), these works use analytical prescriptions to model the AMT produced
by the large-scale magnetic fields in stellar radiative zones. Two distinct
analytical models of the Tayler-Spruit dynamo are used: the original model of
Spruit (2002) and the revised one of Fuller et al. (2019), which tackles
previous criticism of the original model. The prescriptions, however, can not
take into account the strong non-linearity behind the dynamo mechanism, which
impels numerical investigations of its 3D complex dynamics to better
characterize its effects in both astrophysical contexts.

Petitdemange et al. (2023, 2024); Daniel et al. (2023) performed 3D direct
numerical simulations of dynamo action in a stably-stratified Couette flow in
the context of stellar radiative layers and so with a differential rotation in
which the inner core rotates faster than the outer layer. They argued that the
dynamo was driven by the Tayler instability and found scaling laws in
agreement with the prescriptions of Spruit (2002). In our numerical study
Barrère et al. (2023), the setup is different with an outer sphere rotating
faster than the core, which is relevant to our magnetar formation scenario. We
demonstrated the existence of more complex dynamics with two Tayler-Spruit
dynamo branches, which have distinct magnetic field strengths and geometries:
the weaker branch shows a hemispherical field while the strongest one displays
a dipolar symmetry, i.e. the magnetic field is equatorially symmetric.
Furthermore, the former follows the analytical scaling of Spruit (2002), while
the latter is in agreement with the predictions of Fuller et al. (2019).
Lastly, the dipolar dynamo could reach axisymmetric toroidal and dipole
magnetic fields up to $\sim$2\text{\times}{10}^{15}\text{\,}\mathrm{G}$$ and
$\sim$3\text{\times}{10}^{14}\text{\,}\mathrm{G}$$, respectively. Although
such intensities seem relevant to form magnetars, these models considered only
a fixed ratio of the Brunt-Väisälä to the surface angular frequency
$N/\Omega_{o}=0.1$, whereas it is expected to cover the range
$N/\Omega_{o}\in[0.1,10]$ in real PNSs.

Therefore, this article aims at investigating the impact of $N/\Omega_{o}$ on
the dipolar Tayler-Spruit dynamo discovered in Barrère et al. (2023). The
study of the hemispherical dynamo will lead to another paper more focused on
the complex physics behind the Tayler-Spruit dynamos. In the following, Sect.
2 describes the governing equations and the numerical setup of our
simulations. We present the results in Sect. 3, which will be applied to the
question of magnetar formation in Sect. 4. Finally, we discuss the results and
conclude in Sect. 5 and Sect. 6, respectively.

Figure 1: Viscous Elsasser number (and root mean square magnetic field) as a
function of the ratio of the Brunt-Väisälä frequency to the rotation rate at
the outer sphere. Filled and empty markers represent self-sustained and
transient dynamos, respectively. The black dashed vertical line and arrow
indicate the zone in which the fluid is hydrodynamically unstable. The inset
represents a 3D plot of the radial velocity (violet and green isosurfaces are
the positive and negative values, respectively) and the radial magnetic field
(red and blue isosurfaces are the positive and negative values, respectively)
in a run at $Pm=2$ and $N/\Omega_{o}=2$. The grey arrow points to the run
location in the diagram.

##  2 Numerical setup

###  2.1 Governing equations

As in Barrère et al. (2023), we model the PNS interior as a stably stratified
and electrically conducting fluid. We also adopt the Boussinesq approximation
and consider a fluid with a constant density
$\rho=$3.8\text{\times}{10}^{14}\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$,
which corresponds to a PNS with a radius of
$r_{o}=$12\text{\,}\mathrm{k}\mathrm{m}$$ and a mass of
$M=$1.4\text{\,}\mathrm{M}_{\odot}$$. The fluid evolves in a spherical Taylor-
Couette configuration, i.e. between two concentric spheres of radius
$r_{i}=$3\text{\,}\mathrm{k}\mathrm{m}$$ and $r_{o}$ which rotates with the
respective rates $\Omega_{i}=2\pi\times 25\,{\rm rad\,s^{-1}}$ and
$\Omega_{o}=2\pi\times 100\,{\rm rad\,s^{-1}}$. In the reference frame
rotating with the surface at the angular velocity
$\mathbf{\Omega}_{o}=\Omega_{o}\mathbf{e}_{z}$, the Boussinesq MHD equations
read

$\displaystyle\mathbf{\nabla\cdot v}$ $\displaystyle=0\,,$ (1) $\displaystyle
D_{t}\mathbf{v}$ $\displaystyle=-\frac{1}{\rho}{\nabla
p^{\prime}}-2\Omega_{o}\mathbf{e}_{z}\times\mathbf{v}-N^{2}\Theta\mathbf{e}_{r}+\frac{1}{4\pi\rho}(\nabla\times\mathbf{B})\times\mathbf{B}+\nu\Delta\mathbf{v}\,,$
(2) $\displaystyle D_{t}\Theta$ $\displaystyle=\kappa\Delta\Theta\,,$ (3)
$\displaystyle\partial_{t}\mathbf{B}$
$\displaystyle=\nabla\times(\mathbf{u}\times\mathbf{B})+\eta\Delta\mathbf{B}\,,$
(4) $\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=0\,,$ (5)

where $\mathbf{B}$ is the magnetic field, $\mathbf{v}$ is the velocity field,
$p^{\prime}$ is the non-hydrostatic pressure, $\rho$ is the mean density of
the PNS, $g=g_{o}r/r_{o}$ is the gravitation field, and
$\alpha\equiv\rho^{-1}(\partial_{T}\rho)_{p}$ is the thermal expansion
coefficient. $\mathbf{e}_{z}$ and $\mathbf{e}_{r}$ are the unit vectors of the
axial and the spherical radial directions, respectively. $\theta$ is the
buoyancy variable defined by

$\Theta\equiv-\frac{g}{N^{2}}\frac{\rho^{\prime}}{\rho}\,,$ (6)

where $\rho^{\prime}$ is the density perturbation due to the combined effect
of the electron fraction and entropy perturbations and

$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)}\,,$ (7)

is the Brunt-Väisälä frequency with the electron fraction $Y_{e}$, and the
entropy $S$, respectively.

In the above equations, we assume that the magnetic diffusivity $\eta$, the
kinematic viscosity $\nu$, and the “thermal” diffusivity $\kappa$ are
constant. We also assume that the thermal and lepton number diffusivities are
equal, which allows us to describe the buoyancy associated with both entropy
and lepton number gradients with the use of a single buoyancy variable
$\theta$ (Guilet et al., 2015).

Apart from the magnetic diffusivity which relates to the electrical
conductivity of electrons, the physical interpretation of the other transport
coefficients can lead to different estimates, depending on whether neutrinos
are considered or not to be the main source of diffusive processes (see Sect.
1.3 of the supplementary materials in Barrère et al., 2023).

Finally, we apply no-slip, electrically insulating, and fixed buoyancy
variable boundary conditions on both shells.

###  2.2 Numerical methods

We use the open source pseudo-spectral code MagIC111Commit 2266201a5 on
https://github.com/magic-sph/magic (Wicht, 2002; Gastine & Wicht, 2012;
Schaeffer, 2013) to integrate Eqs. (1)–(5) in 3D spherical geometry. To
satisfy the solenoidal conditions (1) and (5), the velocity and magnetic
fields are decomposed in poloidal and toroidal components (Mie
representation),

$\displaystyle\rho\mathbf{u}$ $\displaystyle=\nabla\times(\nabla\times
W\mathbf{e}_{r})+\nabla\times Z\mathbf{e}_{r}\,,$ (8)
$\displaystyle\mathbf{B}$ $\displaystyle=\nabla\times(\nabla\times
b\mathbf{e}_{r})+\nabla\times a_{j}\mathbf{e}_{r}\,,$ (9)

where $W$ and $Z$ ($b$ and $a_{j}$) are the poloidal and toroidal potentials
for the velocity (magnetic) field. The whole system of equations is then
solved in spherical coordinates by expanding the scalar potentials in
Chebyshev polynomials in the radial direction, and spherical harmonic
functions in the angular directions. The time-stepping scheme used is the
implicit/explicit Runge-Kutta BPR353 (Boscarino et al., 2013). We refer the
reader to the MagIC online documentation222https://magic-sph.github.io for an
exhaustive presentation of the numerical techniques.

###  2.3 Input parameters

The resistivity is controlled by the magnetic Prandtl number
$Pm\equiv\nu/\eta$. Though its realistic value in PNSs ($Pm\sim 10^{11}$,
Barrère et al., 2023) can not be reached by numerical simulations, we stay in
the regime $Pm\geqslant 1$ as we impose $Pm\in[1,4]$. We keep fixed the other
dimensionless control parameters: the shell aspect ratio $\chi\equiv
r_{i}/r_{o}=0.25$ and width $d\equiv r_{o}-r_{i}$, the Ekman number
$E\equiv\nu/(d^{2}\Omega_{o})=10^{-5}$, the thermal Prandtl numbers
$Pr\equiv\nu/\kappa=0.1$, and the Rossby number $Ro\equiv
1-\Omega_{i}/\Omega_{o}=0.75$, which controls the imposed differential
rotation.

The imposed stable stratification is characterized by the Brunt-Väisälä
frequency $N$ (Eq. (7)). In our parameter study, the ratio $N/\Omega_{o}$ is
varied between $0.1$ and $10$ and so covers the PNS regime. In practice, this
ratio is related to the Rayleigh number $Ra\equiv-(N/\Omega)^{2}Pr/E^{2}$,
which is negative in the regime of stable stratification.

The resolution is fixed at $(n_{r},n_{\theta},n_{\phi})=(257,256,512)$ for all
the runs. A few simulations were rerun with a higher resolution of
$(n_{r},n_{\theta},n_{\phi})=(481,512,1024)$ but showed no significant change
compared to runs with the usual resolution (see Appendix).

The simulations are initialized either from a nearby saturated state or a
strong ($B_{\phi}=$3.4\text{\times}{10}^{14}\text{\,}\mathrm{G}$$) toroidal
axisymmetric field with a dipolar equatorial symmetry, i.e. equatorially
symmetric333For the choice of these definitions, see Gubbins & Zhang (1993).
with $l=2,m=0$. We define a turbulent resistive time
$\bar{\tau}_{\eta}=\left(\pi r_{o}/\bar{\ell}\right)^{2}/\eta\sim
0.2d^{2}/\eta$, where $\bar{\ell}=10$ is the typical value of the average
harmonic degree of the time-averaged magnetic energy spectrum. In the
following, we will term a solution ‘transient’ when a steady state is
sustained for a time interval $\Delta t>0.3\bar{\tau}_{\eta}$, and ‘stable’
for $\Delta t\geqslant\bar{\tau}_{\eta}$.

We start with the run named ‘Ro0.75s’ from Barrère et al. (2023), where the
stratification is $N/\Omega_{o}=0.1$. The saturated state of this dynamo is
used to initialise the next simulation with a stronger stratification. The
whole set of simulations is initiated similarly using the nearby saturated
state of a less stratified run. With this procedure, $N/\Omega_{o}$ is
increased gradually in order to study the evolution of the dynamo branch.

###  2.4 Output parameters

We first characterize our models by computing the time average of the kinetic
and magnetic energy densities (after filtering out any initial transient). The
latter is expressed in terms of the viscous Elsasser number
$\Lambda_{\nu}\equiv\Lambda/Pm=B_{\rm rms}^{2}/(4\pi\rho\nu\Omega_{o})$ and
used to compute different root-mean-square (RMS) estimates of the magnetic
field. In addition to the total field, we distinguish the poloidal and
toroidal fields based on the Mie representation (Sect. 2.2), while the dipole
field refers to the $l=1$ poloidal component.

Figure 2: Meridional slices of the angular frequency and the non-axisymmetric
radial velocity (left and right slices respectively) for different values of
$N/\Omega_{o}$. $\Omega$ and $v_{r}^{m\neq 0}$ are scaled by $\Omega_{o}$ and
$d\Omega_{o}$, respectively.

##  3 Results

The following sections gather the different results we obtain from the set of
numerical simulations listed in appendix B. We first describe the global
dynamics of the dipolar Tayler-Spruit dynamo in the parameter space in Sect.
3.1. Then, we analyse the influence of stratification on the modes of Tayler
instability and on the generated axisymmetric magnetic fields in their
saturated state in Sect. 3.3 and Sect. 3.4, respectively. We also present the
angular momentum transport by both Reynolds and Maxwell stresses due to the
dynamo and compare the efficiencies of mixing and angular momentum transport
in Sect. 3.5. Finally, we examine a new intermittent behaviour of the Tayler-
Spruit dynamo at $N/\Omega_{o}\geq 2$, which is observed for the first time
(see Sect. 3.6).

###  3.1 Subcritical dynamo sustained at PNS-like stratifications

Fig. 1 shows that a self-sustained Tayler-Spruit dynamo can be maintained up
to $N/\Omega_{o}=1$ for $Pm=1$. For stronger stratifications, we have to
increase $Pm$ (i.e. decrease the resistivity) to maintain the dynamo. For
$Pm=4$, the stationary state is self-sustained up to $N/\Omega_{o}=4$ and we
obtained transient states up to $N/\Omega_{o}=10$. The self-sustained dynamo
is therefore present above the threshold for the fluid to be hydrodynamically
stable at $N/\Omega_{o}\sim 1.5$. This confirms the subcritical nature of the
Tayler-Spruit dynamo, which was already observed in previous studies
(Petitdemange et al., 2023; Barrère et al., 2023). We did not simulate fluids
at greater $Pm$ values for reasons of numerical costs. Given the trend with
$Pm$ observed in our simulations as well as theoretical expectations on the
Tayler instability threshold, we would expect the Tayler-Spruit dynamo to
exist at still higher values of $N/\Omega$ for the higher values of $Pm$
relevant to a PNS.

###  3.2 Impact on the differential rotation

Figure 3: Meridional slices of the axisymmetric azimuthal and the
$s=r\sin\theta$-component of the magnetic field (respective left and right
slices) for increasing values of $N/\Omega_{o}$. The magnetic field is scaled
by $\sqrt{4\pi\rho d^{2}\Omega_{o}^{2}}$. The hatched regions represent
Tayler-stable zones defined by the geometrical criterion of Goossens & Tayler
(1980) (see Supplementary Materials in Barrère et al. (2023)).

The meridional slices of the angular rotation frequency $\Omega$ illustrate
the impact of stable stratification on the rotation profile: we see that the
shear concentrates closer to the inner sphere and increases with
$N/\Omega_{o}$. At the same time, the rotation profile smoothly transits from
a quasi-cylindrical to a spherical geometry, which is an effect already
observed in stably stratified flows. Analytical and numerical studies of these
flows (e.g. Barcilon & Pedlosky, 1967a, b, c; Gaurat et al., 2015; Philidet et
al., 2020) indicate that this transition is controlled by the dimensionless
parameter $Q\equiv Pr(N/\Omega_{o})^{2}$, which varies between
${10}^{-3}\text{\,}$ and $10$ in our set of runs. The change in the flow
geometry is therefore explained by a transition from a case where neither the
rotation nor the buoyancy dominate ($E^{2/3}<Q<1$) to a buoyancy-dominated
flow ($Q\gg 1$).

###  3.3 Impact on the Tayler modes

As seen in Fig. 3, the unstable magnetic modes are located close to the poles
where the latitudinal gradient of $B_{\phi}$ is positive, which is a first
indication of the presence of Tayler modes. To confirm this statement, we use
the geometrical criterion of Goossens & Tayler (1980) for the stability of
$m=1$-modes,

$B_{\phi}^{2}\left(1-2\cos^{2}\theta\right)-\sin\theta\cos\theta\frac{\partial
B_{\phi}^{2}}{\partial\theta}>0\,.$ (10)

The stability regions displayed by the hatched zones in Fig. 3 match very well
regions where the unstable modes are absent. This confirms that the Tayler
instability is clearly identifiable, no matter the values of $N/\Omega_{o}$.

Moreover, the impact of stratification on the mode structure is striking. The
stable stratification tends to stabilise displacements in the radial
direction, as we can see looking at the non-axisymmetric radial velocity
$v_{r}^{m\neq 0}$ field in Fig. 2. As a consequence, the radial length scale
of the instability strongly decreases for increasing values of $N/\Omega_{o}$.
This feature is not surprising because Spruit (1999) already constrained the
mode maximum radial length scale

$l_{\rm TI}<l_{{\rm max},\,N}\equiv r\frac{\omega_{\rm A}}{N}\,,$ (11)

where $\omega_{\rm A}\equiv B_{\phi}^{m=0}/\sqrt{4\pi\rho r^{2}}$ is the
Alfvén frequency. Note that a lower limit due to resistivity is also predicted

$l_{\rm TI}^{2}>l_{{\rm min}}^{2}\equiv\frac{\eta\Omega_{o}}{\omega_{\rm
A}^{2}}\,.$ (12)

The length scales measured in our models are compared to these constraints in
Fig. 4. Since thermal diffusion can mitigate the effect of stratification, we
also define an effective Brunt-Väisälä frequency

$N_{\rm eff}\equiv N\sqrt{\eta/\kappa}=N\sqrt{Pr/Pm}$ (13)

and so

$l_{{\rm max},\,N_{\rm eff}}\equiv r\frac{\omega_{\rm A}}{N_{\rm eff}}\,$ (14)

to take this effect into account (Spruit, 2002). The Tayler modes in our
simulations have length scales ranging from
$r_{o}/4=$3\text{\,}\mathrm{k}\mathrm{m}$$ at $N/\Omega_{o}=0.1$ to
$r_{o}/80=$0.15\text{\,}\mathrm{k}\mathrm{m}$$ at $N/\Omega_{o}=10$. This
implies that the Tayler-Spruit dynamo requires higher and higher resolutions
at greater stratifications to be resolved. The measured $l_{\rm TI}$ follows
very well the upper limit $l_{{\rm max},\,N_{\rm eff}}$ (red points in Fig.
4), but is around one order of magnitude larger than $l_{{\rm max},\,N}$. This
demonstrates the importance of including the mitigation of the stratification
by diffusion. The minimum length scale $l_{\rm min}$ (Eq. (12)) is almost
equal to $l_{\rm TI}$ from $N/\Omega_{o}=0.5$ to $N/\Omega=4$, which indicates
that we are close to the instability threshold. For $N/\Omega_{o}\geqslant 6$,
however, $l_{\rm min}\sim 2l_{{\rm max},\,N_{\rm eff}}\sim 2-3l_{\rm TI}$. The
fluid is therefore stable, which is consistent with the transient state we
find in our simulations. Thus, the analytical limits for the Tayler modes to
develop are validated by our numerical simulations and suggest that the
Tayler-Spruit dynamo could be maintained for $N/\Omega_{o}\in[6,10]$ with
$Pm\gtrsim 16-36$.

Figure 4: Length scale of the Tayler instability mode measured in the
simulations (black stars) as a function of $N/\Omega_{o}$. The theoretical
lower ($l_{\rm min}$ in blue) and upper boundaries of the length scale are
also plotted using the classical ($l_{{\rm max},\,N}$ in orange) and the
effective ($l_{{\rm max},\,N_{\rm eff}}$ in red) Brunt-Väisälä frequencies.
Filled and empty markers represent self-sustained and transient dynamos,
respectively.

In addition to the decrease of $l_{\rm TI}$, the Tayler instability modes are
strongly affected by high values of $N/\Omega_{o}$. The time and volume
averaged spectrum of the magnetic energy in Fig. 5 show that the energy of the
large-scale ($l=1-10$) non-axisymmetric modes (solid lines) drop by two orders
of magnitude between $N/\Omega=0.25$ and at $N/\Omega=2$ compared to the
energy of the dominant axisymmetric toroidal component (blue dotted line).
This difference is represented more quantitatively by comparing the total non-
axisymmetric magnetic field $B^{m\neq 0}_{\rm tot}$ to $B^{m=0}_{\phi}$ in
Fig. 6. The ratio drops from $\sim 1$ to
$\sim$2\text{\times}{10}^{-3}\text{\,}$$ and follows a power law $B^{m\neq
0}_{\rm tot}/B^{m=0}_{\phi}\propto N_{\rm eff}/\Omega_{o}^{-1.8\pm 0.1}$.
Fuller et al. (2019) analytically derived that the ratio between the magnetic
field generated by the Tayler instability (noted $\delta B_{\perp}\sim
B^{m\neq 0}_{\rm tot}$) and $B^{m=0}_{\phi}$ follows $\omega_{\rm
A}/\Omega_{o}$. Since $\omega_{\rm A}\propto N_{\rm eff}/\Omega_{o}^{-1/3}$
(Fuller et al., 2019, and our Sect. 3.4), our simulations therefore do not
match the analytical prediction. Fuller et al. (2019) derived the ratio by
equating the Tayler instability growth rate and a turbulent damping rate
$\omega_{\rm A}^{2}/\Omega_{o}\sim\delta v_{\rm A}/r$, where $\delta v_{\rm
A}\equiv\delta B_{\perp}/\sqrt{4\pi\rho}$. As the growth rate of the Tayler
instability is robust (Zahn et al., 2007; Ma & Fuller, 2019) and well verified
in numerical simulations (Ji et al., 2023), our study then questions the
prediction of the turbulent damping rate.

Figure 5: Time and volume averaged spectra of the magnetic energy for the
parameters $Pm=1,N/\Omega_{o}=0.25$ (top) and $Pm=2,N/\Omega_{o}=2$ (bottom).
The magnetic energy is normalized by the energy of the dominant
$(\ell=2,m=0)$-mode of the toroidal component. Figure 6: Ratio of the RMS non-
axisymmetric magnetic field to the RMS axisymmetric toroidal magnetic field.
The dotted line shows the best fit for a power law of $N_{\rm
eff}/\Omega_{o}$. Filled and empty markers represent self-sustained and
transient dynamos, respectively.

###  3.4 Magnetic field saturation

As in Barrère et al. (2023), we confront the saturated large-scale magnetic
fields in our simulations to the analytical predictions. To this end, we first
measure the impact of the stratification on the local shear rate $q$, which
influences the magnetic field saturation. Indeed, the rotation profiles of
Fig. 2 show that the shear concentrates closer to the inner sphere and
increases with $N/\Omega_{o}$. The quantification of this effect is described
in Appendix A. These larger values of $q$ explain the increase of the magnetic
energy with $N/\Omega_{o}$ observed in Fig. 1.

In order to study the relation of the magnetic field components with $N_{\rm
eff}/\Omega_{o}$ while taking into account the variation of $q$, we use the
analytical prescriptions derived by Fuller et al. (2019):

$\displaystyle B^{m=0}_{\rm tor}$ $\displaystyle\sim\sqrt{4\pi\rho
r_{o}^{2}}\Omega_{o}\left(\frac{q\Omega_{o}}{N_{\rm eff}}\right)^{1/3},$ (15)
$\displaystyle B^{m=0}_{\rm pol},\,B_{\rm dip}$
$\displaystyle\sim\sqrt{4\pi\rho
r_{o}^{2}}\Omega_{o}\left(\frac{q^{2}\Omega_{o}^{5}}{N_{\rm
eff}^{5}}\right)^{1/3}\,.$ (16)

The exponents of $q$ are all the more robust as they are confirmed by
numerical simulations (Barrère et al., 2023). We define dimensionless magnetic
field components compensated for the effect of the shear in the following way:

$\displaystyle B^{m=0}_{\rm tor}$
$\displaystyle\longrightarrow\frac{B^{m=0}_{\rm tor}}{\sqrt{4\pi\rho
r_{o}^{2}\Omega_{o}^{2}}q^{1/3}}$ (17) $\displaystyle B^{m=0}_{\rm
pol},\,B_{\rm dip}$ $\displaystyle\longrightarrow\frac{B^{m=0}_{\rm
pol}}{\sqrt{4\pi\rho r_{o}^{2}\Omega_{o}^{2}}q^{2/3}},\,\frac{B_{\rm
dip}}{\sqrt{4\pi\rho r_{o}^{2}\Omega_{o}^{2}}q^{2/3}}\,.$ (18)

These compensated dimensionless components are plotted in Fig. 7as a function
of $N_{\rm eff}/\Omega_{o}$. The theoretical scaling laws (dotted black lines)
qualitatively match our data. Since the point at $N_{\rm
eff}/\Omega_{o}=$3\text{\times}{10}^{-2}\text{\,}$$ diverges from the scalings
due to the weaker effect of stable stratification, we exclude it while
calculating the best fits. We obtain the following power-laws $B_{\rm
tor}^{m=0}\propto({N_{\rm eff}/\Omega_{o}})^{-0.11\pm 0.05}$, $B^{m=0}_{\rm
pol}\propto({N_{\rm eff}/\Omega_{o}})^{-1.1\pm 0.2}$, and $B_{\rm
dip}\propto({N_{\rm eff}/\Omega_{o}})^{-1.5\pm 0.1}$. While $B_{\rm
tor}^{m=0}$ and $B_{\rm pol}^{m=0}$ follow power-laws slightly less steep than
predicted in Eqs. (15) and (16), $B_{\rm dip}$ is in good agreement with Eq.
(16).

This agreement with the theory is also found for the ratio
$B^{m=0}_{r}/B^{m=0}_{\phi}\sim\omega_{\rm A}/N_{\rm eff}$ (Spruit, 2002;
Fuller et al., 2019) as seen in Fig. 8. Our data is fitted by the power law
$B^{m=0}_{\rm pol}/B^{m=0}_{\rm tor}\propto(\omega_{\rm A}/N_{\rm
eff})^{0.93\pm 0.18}$, which is very close to the prediction. On the other
hand, the ratio of the magnetic dipole to the axisymmetric toroidal field
follows a somewhat steeper scaling law $B_{\rm dip}/B^{m=0}_{\rm
tor}\propto(\omega_{\rm A}/N_{\rm eff})^{1.3\pm 0.1}$.

Figure 7: RMS toroidal and poloidal axisymmetric magnetic fields (top), and
RMS magnetic dipole (bottom) compensated with the measured shear rate as a
function of the ratio between the effective Brunt-Väisälä frequency to the
rotation rate at the outer sphere $N_{\rm eff}/\Omega_{o}$. The magnetic field
is rendered dimensionless and compensated for the effect of the shear using
Eqs. (17) and (18). Dotted lines show the best fits of the data with Fuller’s
theoretical scaling laws (Eqs. (15) and (16)) within a multiplying factor.
Filled and empty markers represent self-sustained and transient dynamos,
respectively. Figure 8: Ratio between the RMS axisymmetric poloidal (top) and
the RMS dipolar (bottom) magnetic fields to the axisymmetric toroidal magnetic
field. Dotted lines show the best fits of the data with Fuller’s theoretical
scaling law $B_{r}/B_{\phi}\sim\omega_{\rm A}/N_{\rm eff}$ within a
multiplying factor. Filled and empty markers represent self-sustained and
transient dynamos, respectively.

###  3.5 Angular momentum transport and mixing

Figure 9: RMS Maxwell (top) and Reynolds (bottom) torques compensated with the
measured shear rate as a function of the ratio between the effective Brunt-
Väisälä frequency to the rotation rate at the outer sphere. Dotted lines shows
the best fits obtained with Fuller’s theoretical scaling laws. Filled and
empty markers represent self-sustained and transient dynamos, respectively.

The angular momentum transport due to the large-scale magnetic field and
turbulence in our simulations is also consistent with the theory of Fuller et
al. (2019), as shown in Fig. 9. For the Maxwell torque $T_{\rm M}$, we find
$B_{s}B_{\phi}\propto({N_{\rm eff}/\Omega_{o}})^{-1.8\pm 0.1}$ and
$B_{s}^{m=0}B_{\phi}^{m=0}\propto({N_{\rm eff}/\Omega_{o}})^{-1.6\pm 0.1}$
depending on whether we take the non-axisymmetric components into account in
$T_{\rm M}$. Note that the torque is more and more dominated by the
axisymmetric magnetic fields as $N_{\rm eff}/\Omega_{o}$ increases. This
dominance was assumed by Fuller et al. (2019) and can be expected given the
results of Sect. 3.3. The Reynolds torque values are more dispersed as a
function of the stratification, but fit the power laws $v^{m\neq
0}_{r}v^{m\neq 0}_{\phi}\propto({N_{\rm eff}/\Omega_{o}})^{-3.5\pm 0.2}$ and
$v^{m\neq 0}_{r}v^{m\neq 0}_{\phi}\propto({N_{\rm eff}/\Omega_{o}})^{-3.4\pm
0.2}$. Despite some scattering at high values of $N_{\rm eff}/\Omega_{o}$ in
the points corresponding to transient dynamos, our data therefore follows well
the analytical predictions $T_{\rm M}\propto({N_{\rm eff}/\Omega_{o}})^{-2}$
and $T_{\rm R}\propto({N_{\rm eff}/\Omega_{o}})^{-10/3}$ (dotted lines in
Figs. 9). Moreover, we find $T_{\rm
M}\sim${10}^{2}\text{\,}$-${10}^{3}\text{\,}$T_{\rm R}$, so the magnetic field
is much more efficient than turbulence at transporting angular momentum.

The mixing processes are also a crucial question in astrophysics, especially
in stars. The Tayler-Spruit dynamo is expected to produce a very limited
mixing efficiency compared to the angular momentum transport (Spruit, 2002;
Fuller et al., 2019). To measure this effect in our simulations, we define the
effective angular momentum transport diffusivity $\nu_{\rm AM}\equiv T_{\rm
M}/(\rho q\Omega_{o})$ and roughly approximate the effective mixing
diffusivity as $\nu_{\rm mix}\equiv q^{-5/3}v^{m\neq 0}l_{\rm TI}$, with the
rms turbulent velocity $v^{m\neq 0}\equiv\sqrt{E_{\rm kin}^{m\neq 0}/(2\rho)}$
calculated from the mean non-axisymmetric energy $E_{\rm kin}^{m\neq 0}$. We
divide by the power law $q^{5/3}$ in the expression of $\nu_{\rm mix}$ to take
into account the variation of $q$ like in Figs 7 and 9.

The ratio $\nu_{\rm mix}/\nu_{\rm AM}$ is plotted in Fig.10 and shows that our
data is in fair agreement with the scaling $\nu_{\rm mix}/\nu_{\rm
AM}\propto{N_{\rm eff}/\Omega_{o}}^{-5/3}$ of Fuller et al. (2019). The power
law $\nu_{\rm mix}/\nu_{\rm AM}\propto{N_{\rm eff}/\Omega_{o}}^{-1.2\pm 0.2}$
best fits our data, which is mildly less steep than predicted. Moreover, our
simulations also confirm that $\nu_{\rm mix}/\nu_{\rm
AM}\sim${10}^{-6}\text{\,}$-${10}^{-3}\text{\,}$\ll 1$ for Tayler-Spruit
dynamo. The use of passive scalars evolving in the velocity field in our
simulations could help measure more precisely $\nu_{\rm mix}$ even though the
approximation we used is satisfactory as a first analysis.

Figure 10: Ratio of the effective mixing diffusivity $\nu_{\rm mix}$ to the
effective angular momentum diffusivity $\nu_{\rm AM}$ as a function of $N_{\rm
eff}/\Omega_{o}$. Filled and empty markers represent self-sustained and
transient dynamos, respectively. Table 1: Table that sums up the theoretical
and measured scaling laws of the different quantities discussed in Sects. 3.4
and 3.5, and the dimensionless normalisation factor $\alpha$ defined by Fuller
et al. (2019) (see Eq. 22)

.  Quantity (dimensionless) Fuller et al. (2019)’s scaling law Best fit
exponent $\alpha$ $B^{m\neq 0}_{\rm tot}/B^{m\neq 0}_{\rm tor}$ $\omega_{\rm
A}/\Omega_{o}$ $(N_{\rm eff}/\Omega_{o})^{-0.18\pm 0.1}$ $B^{m=0}_{\rm
tor}/(\sqrt{4\pi\rho r_{o}^{2}\Omega_{o}^{2}}q^{1/3})$ $\alpha(N_{\rm
eff}/\Omega_{o})^{-1/3}$ $(N_{\rm eff}/\Omega_{o})^{-0.11\pm 0.05}$ $0.017\pm
0.001$ $B^{m=0}_{\rm pol}/(\sqrt{4\pi\rho r_{o}^{2}\Omega_{o}^{2}}q^{2/3})$
$\alpha^{2}(N_{\rm eff}/\Omega_{o})^{-5/3}$ $(N_{\rm eff}/\Omega_{o})^{-1.1\pm
0.2}$ $0.009\pm 0.002$ $B_{\rm dip}/(\sqrt{4\pi\rho
r_{o}^{2}\Omega_{o}^{2}}q^{2/3})$ $\alpha^{2}(N_{\rm eff}/\Omega_{o})^{-5/3}$
$(N_{\rm eff}/\Omega_{o})^{-1.5\pm 0.1}$ $0.007\pm 0.001$ $B^{m=0}_{\rm
pol}/B^{m=0}_{\rm tor}$ $\omega_{\rm A}/N_{\rm eff}$ $(\omega_{\rm A}/N_{\rm
eff})^{0.93\pm 0.2}$ — $B_{\rm dip}/B^{m=0}_{\rm tor}$ $\omega_{\rm A}/N_{\rm
eff}$ $(\omega_{\rm A}/N_{\rm eff})^{1.3\pm 0.1}$ — $B_{s}B_{\phi}/(4\pi\rho
r_{o}^{2}\Omega_{o}^{2}q)$ $\alpha^{3}(N_{\rm eff}/\Omega_{o})^{-2}$ $(N_{\rm
eff}/\Omega_{o})^{-1.8\pm 0.1}$ $0.016\pm 0.004$
$B_{s}^{m=0}B_{\phi}^{m=0}/(4\pi\rho r_{o}^{2}\Omega_{o}^{2}q)$
$\alpha^{3}(N_{\rm eff}/\Omega_{o})^{-2}$ $(N_{\rm eff}/\Omega_{o})^{-1.6\pm
0.1}$ $0.01\pm 0.004$ $v_{r}^{m\neq 0}v_{\phi}^{m\neq
0}/(r_{o}^{2}\Omega_{o}^{2}q^{5/3})$ $(N_{\rm eff}/\Omega_{o})^{-10/3}$
$(N_{\rm eff}/\Omega_{o})^{-3.5\pm 0.2}$ — $v_{s}^{m\neq 0}v_{\phi}^{m\neq
0}/(r_{o}^{2}\Omega_{o}^{2}q^{5/3})$ $(N_{\rm eff}/\Omega_{o})^{-10/3}$
$(N_{\rm eff}/\Omega_{o})^{-3.4\pm 0.2}$ — $\nu_{\rm mix}/\nu_{\rm AM}$
$(N_{\rm eff}/\Omega_{o})^{-5/3}$ $(N_{\rm eff}/\Omega_{o})^{-1.2\pm 0.2}$ —

Table 1 sums up the comparisons we have done between our data and the
different scalings derived by Fuller et al. (2019). Our results thus
consolidate the validity of Fuller et al. (2019)’s formalism for the
saturation of large-scale magnetic fields and angular momentum transport.
Besides, our simulations are not compatible with the analytical prescriptions
of Spruit (2002), which read

$\displaystyle B^{m=0}_{\rm tor}$ $\displaystyle\sim\sqrt{4\pi\rho
r_{o}^{2}}\Omega_{o}\left(\frac{q\Omega_{o}}{N_{\rm eff}}\right)$ (19)
$\displaystyle B^{m=0}_{\rm pol},\,B_{\rm dip}$
$\displaystyle\sim\sqrt{4\pi\rho
r_{o}^{2}}\Omega_{o}\left(\frac{q^{2}\Omega_{o}^{3}}{N_{\rm eff}^{3}}\right)$
(20) $\displaystyle T_{\rm M}$ $\displaystyle\sim
r_{o}^{2}\Omega_{o}^{2}q^{3}\left(\frac{\Omega_{o}}{N_{\rm
eff}}\right)^{4}\,.$ (21)

While our simulations support the scaling law of Fuller et al. (2019), we can
also constrain the dimensionless normalisation factor, (noted $\alpha$ in
Fuller et al. (2019)), that parametrises the saturated strength of the
axisymmetric toroidal magnetic field

$\frac{B^{m=0}_{\rm tor}}{\sqrt{4\pi\rho
r_{o}^{2}}}=\alpha\Omega_{o}\left(\frac{q\Omega_{o}}{N_{\rm
eff}}\right)^{1/3}\,.$ (22)

We infer the value of $\alpha$ by fitting our data by the theoretical scaling
law. The measures are listed in the last column of Table 1 and we find a mean
value of $\alpha\sim${10}^{-2}\text{\,}$$. This value is small compared to
those inferred by adjusting $\alpha$ in 1D stellar evolution models to the
asteroseismic observations of sub-/red giants, which is $\sim 0.25-1$ (Fuller
et al., 2019; Fuller & Lu, 2022; Eggenberger et al., 2019b). Either way, our
numerical simulations provide a more physically motivated value of $\alpha$
that could be implemented in 1D stellar evolution codes including the Tayler-
Spruit dynamo to transport angular momentum.

###  3.6 Intermittency

When $N/\Omega_{o}\geq 2$, we find that the Tayler-Spruit dynamo displays an
intermittent behaviour, which is clearly visible in the time series of Fig. 11
where the non-axisymmetric magnetic energy drops and increases cyclically by
two orders of magnitude. This corresponds to the loss and growth of the Tayler
instability. The same cycle also occurs for the axisymmetric $B_{r}$ and
$B_{\theta}$, which illustrates the loss of the dynamo. Those two cycles show
a very short lag of $\sim$2.4\text{\,}\mathrm{s}$$. We then notice that the
oscillations of the axisymmetric toroidal and poloidal magnetic energies are
in antiphase. This is also observed in the butterfly diagrams in which
$B_{\phi}$ decreases locally, and so in the volume average when $B_{r}$ is the
strongest. These cycles can be interpreted qualitatively as follows:

  * (i)

$B^{m=0}_{\phi}$ is close but above the critical strength for the Tayler
instability derived by combining Eqs. (14) and (12)

$B^{m=0}_{\phi,\rm c}\equiv\sqrt{4\pi\rho
r_{o}^{2}}\Omega_{o}\left(\frac{N_{\rm
eff}}{\Omega_{o}}\right)^{1/2}\left(\frac{\eta}{r_{o}^{2}\Omega_{o}}\right)^{1/4}$
(23)

and the dynamo is acting to generate $B^{m=0}_{r}$ ;

  * (ii)

$B^{m=0}_{\phi}$ decreases slightly below the critical strength due to
turbulent dissipation, which kills the Tayler instability and so the dynamo
loop ;

  * (iii)

the axisymmetric poloidal magnetic energy drops and the axisymmetric toroidal
component increases because of the winding and the lack of turbulent
dissipation ;

  * (iv)

$B^{m=0}_{\phi}$ exceeds the critical strength and the dynamo is active again.

An intermittent Tayler-Spruit dynamo was already proposed by Fuller & Lu
(2022) to explain the angular momentum transport in stellar stellar radiative
regions with a low shear.

Figure 11: Top: Time series of the magnetic energy. Bottom: Butterfly diagram
showing the latitudinal structure time evolution of different axisymmetric
magnetic field components averaged between the radii $r=$5\text{\,}$$ and
$r=$6\text{\,}\mathrm{k}\mathrm{m}$$. The magnetic energy was converted to
physical units by fixing $N=${10}^{-3}\text{\,}\mathrm{s}^{-1}$$ Figure 12:
Period of the cycle $P_{\rm cyc}$ (top) and the duty cycle $\alpha_{\rm cyc}$
(bottom) of the intermittent dynamo as a function of the input $N/\Omega_{o}$.
Filled and empty markers represent self-sustained and transient dynamos,
respectively.

Quantitatively, we find $B_{\phi,\rm
c}^{m=0}\sim$1.4\text{\,}$-$2.1\text{\times}{10}^{15}\text{\,}\mathrm{G}$$ for
the models with $N/\Omega_{o}\in[2,10]$, which is very close to the maximum
values
$B_{\phi}^{m=0}\sim$2.5\text{\,}$-$3\text{\times}{10}^{15}\text{\,}\mathrm{G}$$
measured in the same models. The proximity to the instability threshold
supports our interpretation. To characterise the time evolution of the
intermittency, we measure its duty cycle $\alpha_{\rm cyc}$, i.e. the ratio of
the time when the dynamo is active to the period of the cycle. We find that it
varies between $0.38$ and $0.66$, with a tendency to decrease with
$N/\Omega_{o}$ as seen in Fig. 12. The same trend is observed for the period
of these cycles $P_{\rm cyc}$, which ranges between $3\text{\,}\mathrm{s}$ and
$30\text{\,}\mathrm{s}$. This is consistent with the fact that we get closer
to the dynamo threshold.

##  4 Application to magnetar formation

In this section, we apply our numerical results to the magnetar formation
scenario proposed by Barrère et al. (2022), whose semi-analytical modelling
was based on the formalism by Fuller et al. (2019). To this end, the magnetic
field is converted into physical units by fixing the following parameters to
typical values in PNSs: the PNS radius
$r_{o}=$12\text{\,}\mathrm{k}\mathrm{m}$$, mass
$M=$1.4\text{\,}\mathrm{M}_{\odot}$$ that corresponds to a constant PNS
density of
$\rho\sim$4.1\text{\times}{10}^{14}\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$,
and Brunt-Väisälä frequency $N=$1\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$.
Fig. 13 shows the obtained magnetic field strength as a function of the
angular frequency of the PNS surface, for the axisymmetric toroidal and
poloidal components ($B^{m=0}_{\rm tor}$, $B^{m=0}_{\rm pol}$, upper panel)
and for the dipolar component $B_{\rm dip}$ (lower panel). The red markers
correspond to the magnetic field measured in the simulations, while the blue
markers correspond to the values extrapolated to $q=1$ (as was assumed in
Barrère et al. (2022)). This plot can be compared to Fig. 5 in Barrère et al.
(2022), the main difference being that we define, here, a low-field magnetar
as a magnetar with $B^{m=0}_{\rm tor}\geqslant${10}^{14}\text{\,}\mathrm{G}$$
but $B_{\rm dip}<$4.4\text{\times}{10}^{13}\text{\,}\mathrm{G}$$. The magnetic
field intensity show a similar trend with rotation frequency as was predicted
by Barrère et al. (2022): the axisymmetric toroidal magnetic field increases
as $\propto\Omega^{1.2}$ ($\propto\Omega^{4/3}$ in Barrère et al. (2022))
while the poloidal and dipolar components increase as $\propto\Omega^{2.4}$
($\propto\Omega^{8/3}$ in Barrère et al. (2022)). This agreement is linked to
the fact that the magnetic field in our numerical simulations follows well
Fuller et al. (2019)’s scaling law, as shown in the previous sections. The
main difference is that the saturated magnetic field in our simulations is
$\sim 17$ times weaker than in the model in Barrère et al. (2022). This
difference mainly comes from the assumption in Barrère et al. (2022) of a
dimensionless normalisation factor $\alpha=1$, while our simulations indicate
$\alpha\simeq 0.01$ (see Table 1). Note that this value is smaller than the
inverse of the measured dimensionless normalisation factor $\alpha^{-1}\sim
100$ because Barrère et al. (2022) used $N$ instead of the $N_{\rm eff}$.

The weaker magnetic field in our numerical simulations shifts the upper limit
of the rotation period to form magnetar-like magnetic fields to
$P\lesssim$6\text{\,}\mathrm{m}\mathrm{s}$$. This new limit corresponds to a
lower accreted fallback mass limit of
$\sim$5\text{\times}{10}^{-2}\text{\,}\mathrm{M}_{\odot}$$, which is still
consistent with recent supernova simulations (see the discussion in Barrère et
al., 2022). For rotation periods longer than $6\text{\,}\mathrm{m}\mathrm{s}$,
the magnetic dipole is too weak for a classical magnetar but the Tayler-Spruit
dynamo still produces strong total magnetic fields above
${10}^{14}\text{\,}\mathrm{G}$. We suggest that these could correspond to the
formation of low dipolar field magnetars. Indeed, the observations of
absorption lines in the X-ray spectra of low-field magnetars (Tiengo et al.,
2013; Rodríguez Castillo et al., 2016) and 3D numerical simulations of
magnetic field long-term evolution in NSs (Igoshev et al., 2021) suggest that
such strength is enough to produce magnetar-like luminous activity.

Figure 13: Magnetic strength of the axisymmetric toroidal $B^{m=0}_{\rm tor}$
(pentagons) and poloidal $B^{m=0}_{\rm pol}$ (squares) components (upper
panel), as well as that of the magnetic dipole $B_{\rm dip}$ (triangles, lower
panel) as a function of the angular frequency of the outer sphere, which
represents the PNS surface. The red markers correspond to the magnetic field
measured in the simulations, while the blue markers correspond to the values
extrapolated to $q=1$. The dotted lines are the best power-law fit of the
data. The dark and light grey regions represent the range of magnetic field
for classical magnetars ($B_{\rm
dip}\geqslant$4.4\text{\times}{10}^{13}\text{\,}\mathrm{G}$$) and low-field
magnetars ($B^{m=0}_{\rm tor}\geqslant${10}^{14}\text{\,}\mathrm{G}$$). The
black dashed line and arrow illustrate the rotation period below which the
dynamo can form classical magnetar-like magnetic fields. Filled and empty
markers represent self-sustained and transient dynamos, respectively.

##  5 Discussion

Here, we discuss the simplifications we used for the modelling of the PNS
interior evolution: the mechanism to force the differential rotation (Sect.
5.2) and the Boussinesq approximation (Sect. 5.3). In Sect. 5.4, we finally
compare our results on the Tayler-Spruit dynamo (Barrère et al., 2023, this
article) and the Tayler-Spruit dynamo obtained in other numerical simulations
(Petitdemange et al., 2023; Daniel et al., 2023; Petitdemange et al., 2024).

###  5.1 Dependence on diffusion processes

By assuming that all the rotational energy of the PNS is converted into
kinetic energy of the explosions, observations of SN remnants associated to
magnetars constrained a minimum initial PNS rotation period of
$P_{0}\sim$5\text{\,}\mathrm{m}\mathrm{s}$$ (Vink & Kuiper, 2006). This value
is close to our new constrain of the maximum $P$ to form magnetars through our
scenario, which would leave only a small parameter space for the formation of
magnetars by the Tayler-Spruit dynamo
($$5\text{\,}\mathrm{m}\mathrm{s}$<P<$6\text{\,}\mathrm{m}\mathrm{s}$$).
However, we argue that our simulations may underestimate the strength of
magnetic fields that would be generated in a realistic PNS. First, our
simulations are far from the realistic regime of $Pm$, which reaches
$Pm\sim${10}^{11}\text{\,}$$ assuming a neutrino viscosity
($Pm\sim${10}^{4}\text{\,}$$ assuming a shear viscosity) in PNSs
$\sim$10\text{\,}\mathrm{s}$$ after the core bounce (see Supplementary
Materials in Barrère et al., 2023). Second, the thermal diffusivity is much
larger than the resistivity in a PNS, which implies a small effective Brunt-
Väisälä frequency of $N_{\rm eff}=$2.2\text{\times}{10}^{-8}\text{\,}$N$
assuming a neutrino viscosity ($N_{\rm
eff}=$8\text{\times}{10}^{-4}\text{\,}$N$ assuming a shear viscosity). The
high-$Pm$ regime and the alleviation of the stable stratification favour most
likely the development of the dynamo and stronger magnetic fields. The
extrapolation of our results to the high Pm regime is an important open
question, as it is unclear whether the theoretical scalings will hold in this
regime.

###  5.2 Forcing of the differential rotation

To force the differential rotation, we chose to use a spherical Taylor-Couette
configuration, in which a constant rotation rate is imposed on both inner and
outer spheres. In this setup, the rotation profile is free to evolve as the
angular momentum is transported by turbulence and large-scale magnetic fields.
The imposed rotation of the outer sphere roughly mimics the maintenance of the
surface rotation due to fallback accretion, once the PNS surface is already
spun up significantly. However, the rotation profile evolution does not
describe the beginning of the accretion during which the surface is spun up
and the differential rotation, first concentrated close to the surface, is
transported in the PNS interior.

Maintaining the rotation on both spheres allows us to inject energy into the
flow and try to control the shear rate. As noticed in Sect. 3.4 and quantified
in App. A, the stable stratification however significantly changes the shear
rate. This complicates the measure of the respective scaling exponents with
$N/\Omega_{o}$ and $q$ independently. In addition, we observe in Fig. 2 that
most of the shear is concentrated closer and closer to the inner sphere. As
confirmed by our simulations, this restricts significantly the domain in which
the Tayler-Spruit dynamo can operate and participate to make the dynamo more
difficult to sustain. Thus, to investigate stronger stratification regimes, it
will be necessary to change the forcing method and perhaps opt for a
volumetric forcing as used for instance by Meduri et al. (2024).

###  5.3 Validity of the Boussinesq approximation

To model the PNS interior, we used the Boussinesq approximation, which reduces
the numerical cost and allows us to produce a few tens of models to better
understand the physics of the Tayler-Spruit dynamo. Despite the importance of
the density gradient, this approximation is reasonable in the case of PNS
interior:

  * (i)

The sound speed is close to the speed of light
$c_{s}\sim${10}^{10}\text{\,}\mathrm{c}\mathrm{m}\,\mathrm{s}^{-1}$$
(Hüdepohl, 2014; Pascal, 2021, private communication), so $v_{\rm
A}/c_{s}\lesssim v_{\phi}/c_{s}\lesssim${10}^{-2}\text{\,}$$, where
$v_{a}\equiv r_{o}\omega_{\rm A}$ and $v_{\phi}$ are the typical Alfvén and
azimuthal speeds.

  * (ii)

The density perturbation associated to the buoyancy term is small compared to
the PNS mean density: $\delta\rho/\rho=\theta
N^{2}/g\lesssim$9\text{\times}{10}^{-2}\text{\,}$$, with
$N=${10}^{3}\text{\,}\mathrm{s}^{-1}$$, $g\sim
GM/r_{o}\sim$1.3\text{\times}{10}^{13}\text{\,}\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$$,
and $\theta\lesssim r_{o}$ is the buoyancy variable (Eq. 6).

The impact of density gradient on the Tayler-Spruit dynamo has never been
investigated so far in numerical simulations. Therefore, future work should
consider more realistic PNS density profiles.

###  5.4 Comparison with other numerical models

In the literature, only a few other studies investigate numerically the
Tayler-Spruit dynamo (Petitdemange et al., 2023, 2024; Daniel et al., 2023).
The main difference between our setup and theirs is the opposite shear, i.e.
in their setup the inner boundary rotates faster than the outer one. As in our
studies, they find a subcritical bifurcation at the Tayler instability
threshold to a self-sustained state with a dominant axisymmetric toroidal
magnetic field. However, many differences can be noticed:

  * •

The generated magnetic structure in their simulations has a smaller scale and
is localized near the inner sphere in the equatorial plane. The impact of
stable stratification on the length scale of these modes may deserve a deeper
analysis. It is still unclear why this configuration is stable for one sign of
the shear and not the other.

  * •

As in Barrère et al. (2023), a hemispherical dynamo solution is also found by
Petitdemange et al. (2024) as they move from a laminar dynamo solution to the
strong Tayler-Spruit dynamo by increasing $N/\Omega_{o}$. However, they do not
find bistability between the hemispherical and the strong solutions as in
Barrère et al. (2023).

  * •

While the dipolar and hemispherical dynamos we found in Barrère et al. (2023)
are in good agreement with the predictions of Fuller et al. (2019) and Spruit
(2002), respectively, all their models, including those of the hemispherical
solution, are in agreement with the analytical model of Spruit (2002).

Therefore, the few numerical studies of the Tayler-Spruit dynamo indicate a
much more complex physics than anticipated analytically, with the existence of
a wide variety of dynamo solutions. So far, only Daniel et al. (2023) propose
a non-linear model of the subcritical transition to the Tayler-Spruit dynamo
of Petitdemange et al. (2023). In order to include the other solutions we
discovered, we must further investigate the dynamics of the dynamo using tools
from dynamic system theory.

##  6 Conclusions

###  6.1 Summary

Following our previous study Barrère et al. (2023) of the Tayler-Spruit dynamo
with a fixed ratio of rotation to Brunt-Väisälä frequency, we performed
numerical simulations of the dipolar Tayler-Spruit dynamo to investigate its
dependence on the level of stratification. The main results can be summarized
as follows:

  * •

We find self-sustained Tayler-Spruit dynamos for stratifications up to
$N/\Omega_{o}=4$ (and up to $N/\Omega_{o}=10$ for transient dynamos). The
dynamo also becomes intermittent as the saturated $B^{m=0}_{\phi}$ is close to
the Tayler instability threshold for $N/\Omega_{o}\geqslant 2$.

  * •

We observe a good agreement with the scaling laws of Fuller et al. (2019) for
large-scale magnetic fields and the angular momentum transport, the latter of
which is dominated by Maxwell torques.

  * •

With increasing $N/\Omega_{o}$, the Tayler modes have reduced radial length
scales as expected but their energy decreases faster than theoretically
predicted, which may indicate an underestimation of the turbulent dissipation
by the analytical models.

  * •

By measuring an approximate mixing diffusivity, we also determine the
efficiency of the mixing process due to the Tayler-Spruit dynamo. We find that
mixing is far less efficient than the angular momentum transport, as
analytically predicted.

  * •

Finally, as Fuller et al. (2019), we have defined a dimensionless
normalisation factor $\alpha$ parameterising the scaling law of $B_{\rm
tor}^{m=0}$ and numerically constrained its value, which is
$\alpha\sim${10}^{-2}\text{\,}$$. Therefore, the large-scale magnetic fields
in our simulations are weaker than theoretically foreseen.

Applying these results to the magnetar formation scenario of Barrère et al.
(2022), the lower limit of the angular frequency to form classical magnetar-
like dipoles is found to be larger than derived in Barrère et al. (2022) with
a rotation period of $P\sim$6\text{\,}\mathrm{m}\mathrm{s}$$. This value
corresponds to an accreted fallback mass of
$\sim$5\text{\times}{10}^{-2}\text{\,}\mathrm{M}_{\odot}$$, which is still
reasonable according to CCSN simulations (e.g. Sukhbold et al., 2016, 2018;
Chan et al., 2020; Janka et al., 2022). This new constraint is close to the
minimum initial PNS period derived from the kinetic energy of SN remnants
associated to magnetars (Vink & Kuiper, 2006). Nevertheless, as discussed in
Sect. 5.1, stronger magnetic fields are likely to form in the realistic PNS
regime, alleviating this tension.

###  6.2 Long-term evolution of the magnetic field

After $\sim$100\text{\,}\mathrm{s}$$, the fallback accretion becomes too weak
to maintain the differential rotation in the PNS. The newly formed strong
large-scale magnetic fields transport the angular momentum efficiently, which
damps the differential rotation and the dynamo will eventually stop. The
magnetic field is expected to enter a relaxation phase in which its structure
changes to reach a stable configuration. The exact shape of this magnetic
field is still an open question and, more generally, the magnetic relaxation
problem in astrophysics remains debated (e.g. Braithwaite, 2006; Duez &
Mathis, 2010; Duez et al., 2010; Akgün et al., 2013; Becerra et al., 2022a,
b). It is however well acknowledged that the magnetic configuration is
complex, mixing both large-scale poloidal and toroidal components. Thus, 3D
numerical simulations including rotation and thermal/density stratifications
are required to investigate this stage of the PNS evolution.

On longer timescales of $\sim 1-$100\text{\,}\mathrm{k}\mathrm{y}\mathrm{r}$$,
the remaining strong toroidal magnetic fields located in the NS crust are
prone to Hall diffusion and instability (Rheinhardt & Geppert, 2002), which
modifies their structures and so can influence the magnetar emission. The
strong magnetic field-induced stresses could also cause failures or plastic
deformations, which are suspected to explain the origin of magnetar bursts
(e.g. Thompson & Duncan, 1995; Perna & Pons, 2011; Lander et al., 2015; Lander
& Gourgouliatos, 2019). It is therefore crucial to run 3D numerical
simulations of magnetic field evolution in a NS structure using dynamo-
generated initial magnetic configuration to better constrain these properties.
Further investigations could also include the relaxation of the dynamo-
generated magnetic field to a stable configuration before the PNS becomes a
cooled stable NS.

###  6.3 Interaction with a remaining fallback disc

The magnetic dipole generated by the dipolar Tayler-Spruit dynamo may not be
strong enough to spin the magnetar down to the observed
$8-$12\text{\,}\mathrm{s}$$ via the magnetic spin-down mechanism. A good
alternative would be the propeller mechanism (e.g. Gompertz et al., 2014;
Beniamini et al., 2019; Lin et al., 2021; Ronchi et al., 2022). This operates
when the magnetosphere is large enough to interact with the remaining fallback
disc, i.e. when the Alfvén radius is larger than the corotation radius. In the
propeller regime, the inner disk matter is accelerated to super-Keplerian
velocity, which produces an outflow and so an angular momentum transfer from
the magnetar to the disc. If this mechanism operates in some magnetars, the
magnetic dipole which is inferred from the values of the NS rotation period
and its associated derivative will be overestimated. It thus fosters numerical
studies of the fallback matter in 3D simulations of core-collapse SNe and
investigations on the evolution of the potential remaining disc. This will
help constrain which progenitors are the best candidates to form magnetars via
our fallback scenario.

###  6.4 Implications for stellar physics

Our findings are also of importance for the study of stellar radiative zones.
Indeed, the scaling laws and the dimensionless normalisation factor $\alpha$
derived from our simulations could be implemented in 1D stellar evolution
codes. Evolution models using the prescriptions of Fuller et al. (2019) have
already been computed for sub-giant/red giant stars but with larger values of
$\alpha\sim 0.25-1$. These studies find a strong flattening of the rotation
profile and conclude that the prescribed Tayler-Spruit dynamo can not explain
both rotation profiles of sub-giant and red giant stars (Eggenberger et al.,
2019b), which suggests that different angular momentum transport mechanisms
occur during these two phases (Eggenberger et al., 2019a). The future
asteroseismic measurements of the magnetic fields in stellar interiors with
PLATO will be crucial to clarify the question of the transport mechanisms.
Though the first measurements of magnetic fields in some red giant cores
suggest a strong fossil field (Li et al., 2022, 2023; Deheuvels et al., 2023),
it will be essential to infer the asteroseismic signature of magnetic fields
generated by the simulated Tayler-Spruit dynamos for the future observations.
Evolution models including MHD instabilities effects were also performed in
the case of massive stars to constrain the rotation rate of the remaining PNS
or black hole (Griffiths et al., 2022; Fuller & Lu, 2022). They suggest that
the angular momentum transport by MHD instabilities is significant in every
stage of the massive star evolution. This stresses the importance of
performing 3D anelastic simulations with realistic background profiles of
radiative zones at different evolution stages to better constrain the angular
momentum transport and infer more robust rotation rates of stellar cores.

###### Acknowledgements.

This work was supported by the European Research Council (MagBURST grant
715368), and the PNPS and PNHE programs of CNRS/INSU, co-funded by CEA and
CNES. Numerical simulations have been carried out at the CINES on the Jean-Zay
supercomputer and at the TGCC on the supercomputer IRENE-ROME (DARI project
A130410317).

## References

  * Akgün et al. (2013) Akgün, T., Reisenegger, A., Mastrano, A., & Marchant, P. 2013, MNRAS, 433, 2445
  * Barcilon & Pedlosky (1967a) Barcilon, V. & Pedlosky, J. 1967a, Journal of Fluid Mechanics, 29, 1–16
  * Barcilon & Pedlosky (1967b) Barcilon, V. & Pedlosky, J. 1967b, Journal of Fluid Mechanics, 29, 673–690
  * Barcilon & Pedlosky (1967c) Barcilon, V. & Pedlosky, J. 1967c, Journal of Fluid Mechanics, 29, 609–621
  * Barrère et al. (2023) Barrère, P., Guilet, J., Raynaud, R., & Reboul-Salze, A. 2023, MNRAS, 526, L88
  * Barrère et al. (2022) Barrère, P., Guilet, J., Reboul-Salze, A., Raynaud, R., & Janka, H. T. 2022, A&A, 668, A79
  * Becerra et al. (2022a) Becerra, L., Reisenegger, A., Valdivia, J. A., & Gusakov, M. 2022a, MNRAS, 517, 560
  * Becerra et al. (2022b) Becerra, L., Reisenegger, A., Valdivia, J. A., & Gusakov, M. E. 2022b, MNRAS, 511, 732
  * Beniamini et al. (2019) Beniamini, P., Hotokezaka, K., van der Horst, A., & Kouveliotou, C. 2019, MNRAS, 487, 1426
  * Bochenek et al. (2020) Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59
  * Boscarino et al. (2013) Boscarino, S., Pareschi, L., & Russo, G. 2013, SIAM Journal on Scientific Computing, 35, A22
  * Braithwaite (2006) Braithwaite, J. 2006, A&A, 449, 451
  * Bugli et al. (2023) Bugli, M., Guilet, J., Foglizzo, T., & Obergaulinger, M. 2023, MNRAS, 520, 5622
  * Bugli et al. (2021) Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443
  * Bugli et al. (2020) Bugli, M., Guilet, J., Obergaulinger, M., Cerdá-Durán, P., & Aloy, M. A. 2020, MNRAS, 492, 58
  * Burrows et al. (2007) Burrows, A., Dessart, L., Livne, E., Ott, C. D., & Murphy, J. 2007, ApJ, 664, 416
  * Cantiello et al. (2014) Cantiello, M., Mankovich, C., Bildsten, L., Christensen-Dalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93
  * Chan et al. (2020) Chan, C., Müller, B., & Heger, A. 2020, MNRAS, 495, 3751
  * CHIME/FRB Collaboration et al. (2020) CHIME/FRB Collaboration, Andersen, B. C., Bandura, K. M., et al. 2020, Nature, 587, 54
  * Coti Zelati et al. (2021) Coti Zelati, F., Borghese, A., Rea, N., et al. 2021, The Astronomer’s Telegram, 14674, 1
  * Coti Zelati et al. (2018) Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961
  * Daniel et al. (2023) Daniel, F., Petitdemange, L., & Gissinger, C. 2023, Physical Review Fluids, 8, 123701
  * Deheuvels et al. (2015) Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96
  * Deheuvels et al. (2014) Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27
  * Deheuvels et al. (2023) Deheuvels, S., Li, G., Ballot, J., & Lignières, F. 2023, A&A, 670, L16
  * den Hartogh et al. (2020) den Hartogh, J. W., Eggenberger, P., & Deheuvels, S. 2020, A&A, 634, L16
  * Dessart et al. (2008) Dessart, L., Burrows, A., Livne, E., & Ott, C. D. 2008, ApJ, 673, L43
  * Dessart et al. (2012) Dessart, L., Hillier, D. J., Waldman, R., Livne, E., & Blondin, S. 2012, MNRAS, 426, L76
  * Drout et al. (2011) Drout, M. R., Soderberg, A. M., Gal-Yam, A., et al. 2011, ApJ, 741, 97
  * Duez et al. (2010) Duez, V., Braithwaite, J., & Mathis, S. 2010, ApJ, 724, L34
  * Duez & Mathis (2010) Duez, V. & Mathis, S. 2010, A&A, 517, A58
  * Duncan & Thompson (1992) Duncan, R. C. & Thompson, C. 1992, ApJ, 392, L9
  * Eggenberger et al. (2019a) Eggenberger, P., Deheuvels, S., Miglio, A., et al. 2019a, A&A, 621, A66
  * Eggenberger et al. (2019b) Eggenberger, P., den Hartogh, J. W., Buldgen, G., et al. 2019b, A&A, 631, L6
  * Eggenberger et al. (2005) Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9
  * Evans et al. (1980) Evans, W. D., Klebesadel, R. W., Laros, J. G., et al. 1980, ApJ, 237, L7
  * Ferrario & Wickramasinghe (2006) Ferrario, L. & Wickramasinghe, D. 2006, MNRAS, 367, 1323
  * Fuller & Lu (2022) Fuller, J. & Lu, W. 2022, MNRAS, 511, 3951
  * Fuller et al. (2019) Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661
  * Gastine & Wicht (2012) Gastine, T. & Wicht, J. 2012, Icarus, 219, 428
  * Gaurat et al. (2015) Gaurat, M., Jouve, L., Lignières, F., & Gastine, T. 2015, A&A, 580, A103
  * Gehan et al. (2018) Gehan, C., Mosser, B., Michel, E., Samadi, R., & Kallinger, T. 2018, A&A, 616, A24
  * Gompertz & Fruchter (2017) Gompertz, B. & Fruchter, A. 2017, ApJ, 839, 49
  * Gompertz et al. (2014) Gompertz, B. P., O’Brien, P. T., & Wynn, G. A. 2014, MNRAS, 438, 240
  * Goossens et al. (1981) Goossens, M., Biront, D., & Tayler, R. J. 1981, Astrophys. Space Sci., 75, 521
  * Goossens & Tayler (1980) Goossens, M. & Tayler, R. J. 1980, MNRAS, 193, 833
  * Gotz et al. (2006) Gotz, D., Israel, G. L., Mereghetti, S., et al. 2006, The Astronomer’s Telegram, 953, 1
  * Griffiths et al. (2022) Griffiths, A., Eggenberger, P., Meynet, G., Moyano, F., & Aloy, M.-Á. 2022, A&A, 665, A147
  * Gubbins & Zhang (1993) Gubbins, D. & Zhang, K. 1993, Physics of the Earth and Planetary Interiors, 75, 225
  * Guilet et al. (2015) Guilet, J., Müller, E., & Janka, H.-T. 2015, MNRAS, 447, 3992
  * Guilet et al. (2022) Guilet, J., Reboul-Salze, A., Raynaud, R., Bugli, M., & Gallet, B. 2022, MNRAS, 516, 4346
  * Hu & Lou (2009) Hu, R.-Y. & Lou, Y.-Q. 2009, MNRAS, 396, 878
  * Hüdepohl (2014) Hüdepohl, L. 2014, PhD thesis, Technical University of Munich, Germany
  * Hurley et al. (2005) Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098
  * Hurley et al. (1999) Hurley, K., Cline, T., Mazets, E., et al. 1999, Nature, 397, 41
  * Igoshev et al. (2021) Igoshev, A. P., Hollerbach, R., Wood, T., & Gourgouliatos, K. N. 2021, Nature Astronomy, 5, 145
  * Inserra et al. (2013) Inserra, C., Smartt, S. J., Jerkstrand, A., et al. 2013, ApJ, 770, 128
  * Janka et al. (2022) Janka, H.-T., Wongwathanarat, A., & Kramer, M. 2022, ApJ, 926, 9
  * Ji et al. (2023) Ji, S., Fuller, J., & Lecoanet, D. 2023, MNRAS, 521, 5372
  * Kasen & Bildsten (2010) Kasen, D. & Bildsten, L. 2010, ApJ, 717, 245
  * Kiuchi et al. (2024) Kiuchi, K., Reboul-Salze, A., Shibata, M., & Sekiguchi, Y. 2024, Nature Astronomy, 8, 298
  * Kouveliotou et al. (1994) Kouveliotou, C., Fishman, G. J., Meegan, C. A., et al. 1994, Nature, 368, 125
  * Kuroda et al. (2020) Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, ApJ, 896, 102
  * Lander et al. (2015) Lander, S. K., Andersson, N., Antonopoulou, D., & Watts, A. L. 2015, MNRAS, 449, 2047
  * Lander & Gourgouliatos (2019) Lander, S. K. & Gourgouliatos, K. N. 2019, MNRAS, 486, 4130
  * Li et al. (2022) Li, G., Deheuvels, S., Ballot, J., & Lignières, F. 2022, Nature, 610, 43
  * Li et al. (2023) Li, G., Deheuvels, S., Li, T., Ballot, J., & Lignières, F. 2023, A&A, 680, A26
  * Lin et al. (2021) Lin, W., Wang, X., Wang, L., & Dai, Z. 2021, ApJ, 914, L2
  * Lü & Zhang (2014) Lü, H.-J. & Zhang, B. 2014, ApJ, 785, 74
  * Ma & Fuller (2019) Ma, L. & Fuller, J. 2019, MNRAS, 488, 4338
  * Makarenko et al. (2021) Makarenko, E. I., Igoshev, A. P., & Kholtygin, A. F. 2021, MNRAS, 504, 5813
  * Martin et al. (2014) Martin, J., Rea, N., Torres, D. F., & Papitto, A. 2014, MNRAS, 444, 2910
  * Masada et al. (2022) Masada, Y., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 75
  * Meduri et al. (2024) Meduri, D. G., Jouve, L., & Lignières, F. 2024, A&A, 683, A12
  * Mereghetti et al. (2020) Mereghetti, S., Savchenko, V., Ferrigno, C., et al. 2020, ApJ, 898, L29
  * Metzger et al. (2018) Metzger, B. D., Beniamini, P., & Giannios, D. 2018, ApJ, 857, 95
  * Metzger et al. (2011) Metzger, B. D., Giannios, D., Thompson, T. A., Bucciantini, N., & Quataert, E. 2011, MNRAS, 413, 2031
  * Mosser et al. (2012) Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10
  * Mösta et al. (2014) Mösta, P., Richers, S., Ott, C. D., et al. 2014, ApJ, 785, L29
  * Nicholl et al. (2013) Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2013, Nature, 502, 346
  * Nomoto et al. (2011) Nomoto, K., Maeda, K., Tanaka, M., & Suzuki, T. 2011, Astrophys. Science Sci., 336, 129
  * Obergaulinger & Aloy (2020) Obergaulinger, M. & Aloy, M. Á. 2020, MNRAS, 492, 4613
  * Obergaulinger & Aloy (2021) Obergaulinger, M. & Aloy, M. Á. 2021, MNRAS, 503, 4942
  * Obergaulinger & Aloy (2022) Obergaulinger, M. & Aloy, M. Á. 2022, MNRAS, 512, 2489
  * Obergaulinger et al. (2009) Obergaulinger, M., Cerdá-Durán, P., Müller, E., & Aloy, M. A. 2009, A&A, 498, 241
  * Olausen & Kaspi (2014) Olausen, S. A. & Kaspi, V. M. 2014, ApJS, 212, 6
  * Pascal (2021) Pascal, A. 2021, Theses, Université Paris sciences et lettres
  * Perna & Pons (2011) Perna, R. & Pons, J. A. 2011, ApJ, 727, L51
  * Petitdemange et al. (2023) Petitdemange, L., Marcotte, F., & Gissinger, C. 2023, Science, 379, 300
  * Petitdemange et al. (2024) Petitdemange, L., Marcotte, F., Gissinger, C., & Daniel, F. 2024, A&A, 681, A75
  * Philidet et al. (2020) Philidet, J., Gissinger, C., Lignières, F., & Petitdemange, L. 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 336
  * Raynaud et al. (2022) Raynaud, R., Cerdá-Durán, P., & Guilet, J. 2022, MNRAS, 509, 3410
  * Raynaud et al. (2020) Raynaud, R., Guilet, J., Janka, H.-T., & Gastine, T. 2020, Sci. Adv., 6, eaay2732
  * Rea et al. (2010) Rea, N., Esposito, P., Turolla, R., et al. 2010, Science, 330, 944
  * Rea et al. (2012) Rea, N., Israel, G. L., Esposito, P., et al. 2012, ApJ, 754, 27
  * Rea et al. (2013) Rea, N., Israel, G. L., Pons, J. A., et al. 2013, ApJ, 770, 65
  * Rea et al. (2014) Rea, N., Viganò, D., Israel, G. L., Pons, J. A., & Torres, D. F. 2014, ApJ, 781, L17
  * Reboul-Salze et al. (2021) Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2021, A&A, 645, A109
  * Reboul-Salze et al. (2022) Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2022, A&A, 667, A94
  * Rembiasz et al. (2017) Rembiasz, T., Obergaulinger, M., Cerdá-Durán, P., Aloy, M.-Á., & Müller, E. 2017, ApJS, 230, 18
  * Rheinhardt & Geppert (2002) Rheinhardt, M. & Geppert, U. 2002, Phys. Rev. Lett., 88, 101103
  * Rodríguez Castillo et al. (2016) Rodríguez Castillo, G. A., Israel, G. L., Tiengo, A., et al. 2016, MNRAS, 456, 4145
  * Ronchi et al. (2022) Ronchi, M., Rea, N., Graber, V., & Hurley-Walker, N. 2022, ApJ, 934, 184
  * Schaeffer (2013) Schaeffer, N. 2013, Geochemistry, Geophysics, Geosystems, 14, 751
  * Schneider et al. (2019) Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211
  * Shenar et al. (2023) Shenar, T., Wade, G. A., Marchant, P., et al. 2023, Science, 381, 761
  * Spruit (1999) Spruit, H. C. 1999, A&A, 349, 189
  * Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923
  * Stockinger et al. (2020) Stockinger, G., Janka, H. T., Kresse, D., et al. 2020, MNRAS, 496, 2039
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38
  * Sukhbold et al. (2018) Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93
  * Svinkin et al. (2021) Svinkin, D., Frederiks, D., Hurley, K., et al. 2021, Nature, 589, 211
  * Takiwaki et al. (2009) Takiwaki, T., Kotake, K., & Sato, K. 2009, ApJ, 691, 1360
  * Tayler (1973) Tayler, R. J. 1973, MNRAS, 161, 365
  * Thompson & Duncan (1993) Thompson, C. & Duncan, R. C. 1993, ApJ, 408, 194
  * Thompson & Duncan (1995) Thompson, C. & Duncan, R. C. 1995, MNRAS, 275, 255
  * Tiengo et al. (2013) Tiengo, A., Esposito, P., Mereghetti, S., et al. 2013, Nature, 500, 312
  * Tsuzuki et al. (2024) Tsuzuki, Y., Totani, T., Hu, C.-P., & Enoto, T. 2024, MNRAS, 530, 1885
  * Vink & Kuiper (2006) Vink, J. & Kuiper, L. 2006, MNRAS, 370, L14
  * White et al. (2022) White, C. J., Burrows, A., Coleman, M. S. B., & Vartanyan, D. 2022, ApJ, 926, 111
  * Wicht (2002) Wicht, J. 2002, Physics of the Earth and Planetary Interiors, 132, 281
  * Woosley (2010) Woosley, S. E. 2010, ApJ, 719, L204
  * Woosley & Bloom (2006) Woosley, S. E. & Bloom, J. S. 2006, ARAA, 44, 507
  * Younes et al. (2017) Younes, G., Kouveliotou, C., Jaodand, A., et al. 2017, ApJ, 847, 85
  * Zahn et al. (2007) Zahn, J. P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145
  * Zhang & Mészáros (2001) Zhang, B. & Mészáros, P. 2001, ApJ, 552, L35
  * Zhou et al. (2019) Zhou, P., Vink, J., Safi-Harb, S., & Miceli, M. 2019, A&A, 629, A51
  * Zhu et al. (2023) Zhu, W., Xu, H., Zhou, D., et al. 2023, Science Advances, 9, eadf6198

##  Appendix A Measure of the shear rate

The differential rotation is characterized by a dimensionless shear rate
$q=r\partial_{r}{\rm ln}{\Omega}$. We define an effective shear rate based on
the time average of the radial rotation profile in the saturated state at the
colatitude of $\theta=\pi/8\,{\rm rad}$. We measure an average slope in the
range of radii where half of the Tayler mode energy (approximated by the
latitudinal magnetic energy $E_{B_{\theta}}$) is concentrated around its
maximum. We chose this particular method because this range of radii is the
region where the dynamo occurs. The measures are displayed in Fig. 14 (red
plot) along with other measures made with different methods. Whatever the
method used, we see that all the measures follow the same trend with an
increase of $q\propto N$ until $N/\Omega_{o}=4$ after which the values of $q$
stay almost constant.

Figure 14: Shear rates $q$ measured locally in the simulations as a function
of $N/\Omega_{o}$. The different colours represent distinct methods to measure
$q$: slope in the rotation profile between $3.2\text{\,}$ and
$4\text{\,}\mathrm{k}\mathrm{m}$ (pink), $q$ at the maximum of
$B^{m=0}_{\phi}$ and $B_{\theta}$ (green and orange, respectively), and slope
in the range of radii where half of the Tayler mode energy (approximated by
the latitudinal magnetic energy $E_{B_{\theta}}$) is concentrated around its
maximum (red).

##  Appendix B List of models

Tables 2–4 summarize the key parameters and output quantities of the
simulations carried out in this study.

Table 2: Overview of the stable (or failed) dynamo solutions. All the simulations have the same aspect ratio $\chi=0.25$, Ekman number $E=10^{-5}$, Rossby number $Ro=0.75$, thermal and magnetic Prandtl numbers $Pr=0.1$ and $Pm=1$, and the same resolution $(n_{r},n_{\theta},n_{\phi})=(256,256,512)$. Note that the run named Ro0.75s is the same as in Barrère et al. (2023). This table displays the input parameter of the runs. Name | $Pm$ | ${N/\Omega_{o}}$ | ${N_{\rm eff}/\Omega_{o}}$ | ${\Lambda_{i}}$
---|---|---|---|---
Ro0.75s | $1$ | $0.1$ | $0.03$ | $10$
Pm1Pr0.1NO0.25 | $1$ | $0.25$ | $0.08$ | $\Lambda(\mathrm{Ro0.75s})$
Pm1Pr0.1NO0.5 | $1$ | $0.5$ | $0.16$ | $\Lambda(\mathrm{Ro0.75s})$
Pm1Pr0.1NO1 | $1$ | $1$ | $0.32$ | $\Lambda(\mathrm{Pm1Pr0.1NO0.5})$
Pm2Pr0.1NO2 | $2$ | $2$ | $0.45$ | $\Lambda(\mathrm{Pm1Pr0.1NO1})$
Pm2Pr0.1NO4 | $2$ | $4$ | $0.89$ | $\Lambda(\mathrm{Pm2Pr0.1NO2})$
Pm2Pr0.1NO8 | $2$ | $8$ | $1.79$ | $\Lambda(\mathrm{Pm2Pr0.1NO4})$
Pm4Pr0.1NO4 | $4$ | $4$ | $0.63$ | $\Lambda(\mathrm{Pm2Pr0.1NO4})$
Pm4Pr0.1NO6 | $4$ | $6$ | $0.95$ | $\Lambda(\mathrm{Pm4Pr0.1NO4})$
Pm4Pr0.1NO8 | $4$ | $8$ | $1.26$ | $\Lambda(\mathrm{Pm2Pr0.1NO8})$
Pm4Pr0.1NO10 | $4$ | $10$ | $1.58$ | $\Lambda(\mathrm{Pm2Pr0.1NO10})$
Table 3: Same as Table 2 but this table displays measured values in the simulations used to produce the plots of the paper. Name | ${q}$ | $\Lambda$ | $B_{\mathrm{tor}}^{m=0}$ | $B_{\mathrm{pol}}^{m=0}$ | $B_{\mathrm{dip}}$ | $B_{\mathrm{tot}}^{m\neq 0}$ | $v_{\mathrm{tot}}^{m\neq 0}$
---|---|---|---|---|---|---|---
|  |  | $\left[10^{-3}\sqrt{4\pi\rho d^{2}\Omega_{o}}\right]$ | $\left[10^{-3}\sqrt{4\pi\rho d^{2}\Omega_{o}}\right]$ | $\left[10^{-3}\sqrt{4\pi\rho d^{2}\Omega_{o}}\right]$ | $\left[10^{-3}\sqrt{4\pi\rho d^{2}\Omega_{o}}\right]$ | $[10^{-4}d\Omega_{o}]$
Ro0.75s | $0.06$ | $51.85$ | $21$ | $1.7$ | $1.2$ | $9.3$ | $35$
Pm1Pr0.1NO0.25 | $0.17$ | $69.03$ | $24$ | $1.8$ | $1.4$ | $9.7$ | $34$
Pm1Pr0.1NO0.5 | $0.35$ | $74.71$ | $26$ | $2.0$ | $0.87$ | $6.6$ | $36$
Pm1Pr0.1NO1 | $0.69$ | $99.44$ | $31$ | $0.64$ | $0.36$ | $2.6$ | $19$
Pm2Pr0.1NO2 | $2.37$ | $251.25$ | $36$ | $0.82$ | $0.31$ | $1.1$ | $6.0$
Pm2Pr0.1NO4 | $2.57$ | $332.08$ | $41$ | $0.49$ | $0.18$ | $0.61$ | $3.5$
Pm2Pr0.1NO8 | $2.54$ | $531.21$ | $52$ | $0.18$ | $0.10$ | $0.59$ | $3.4$
Pm4Pr0.1NO4 | $2.57$ | $478.89$ | $35$ | $0.46$ | $0.077$ | $0.36$ | $1.6$
Pm4Pr0.1NO6 | $2.59$ | $495.22$ | $35$ | $0.31$ | $0.077$ | $0.31$ | $1.8$
Pm4Pr0.1NO8 | $2.54$ | $714.16$ | $42$ | $0.31$ | $0.077$ | $0.28$ | $1.6$
Pm4Pr0.1NO10 | $3.74$ | $862.86$ | $47$ | $0.15$ | $0.051$ | $0.2$ | $1.1$
Table 4: Following of Table 3. Name | $B_{s}B_{\phi}/4\pi$ | $B_{s}^{m=0}B_{\phi}^{m=0}/4\pi$ | $v_{r}^{m\neq 0}v_{\phi}^{m\neq 0}/4\pi$ | $v_{s}^{m\neq 0}v_{\phi}^{m\neq 0}/4\pi$ | $l_{\mathrm{TI}}$ | $P_{\mathrm{cyc}}$ | $\alpha_{\mathrm{cyc}}$
---|---|---|---|---|---|---|---
| $[10^{-6}\times 4\pi\rho d^{2}\Omega_{o}^{2}]$ | $[10^{-6}\times 4\pi\rho d^{2}\Omega_{o}^{2}]$ | $[10^{-8}\times\rho d^{2}\Omega_{o}^{2}]$ | $[10^{-8}\times\rho d^{2}\Omega_{o}^{2}]$ | [km] | [s-1] |
Ro0.75s | $15$ | $2.9$ | $52$ | $76$ | $2.8$ | – | –
Pm1Pr0.1NO0.25 | $9.1$ | $2.4$ | $56$ | $89$ | $2.5$ | – | –
Pm1Pr0.1NO0.5 | $5.7$ | $1.8$ | $58$ | $110$ | $1.5$ | – | –
Pm1Pr0.1NO1 | $3.0$ | $0.84$ | $18$ | $35$ | $0.8$ | – | –
Pm2Pr0.1NO2 | $2.2$ | $1.7$ | $2.8$ | $10$ | $0.6$ | $27$ | $0.66$
Pm2Pr0.1NO4 | $1.5$ | $1.4$ | $7$ | $4.4$ | $0.4$ | $12$ | $0.59$
Pm2Pr0.1NO8 | $1.7$ | $1.0$ | $0.52$ | $5.3$ | $0.2$ | $13$ | $0.38$
Pm4Pr0.1NO4 | $0.45$ | $0.63$ | $0.16$ | $0.61$ | $0.4$ | $3$ | $0.62$
Pm4Pr0.1NO6 | $0.53$ | $0.67$ | $0.17$ | $0.94$ | $0.3$ | $16$ | $0.5$
Pm4Pr0.1NO8 | $0.65$ | $0.75$ | $0.12$ | $0.77$ | $0.2$ | $10$ | $0.47$
Pm4Pr0.1NO10 | $0.50$ | $0.52$ | $0.066$ | $0.46$ | $0.15$ | $6$ | $0.57$


IEEEexample:BSTcontrol

# Performance or Trust? Why Not Both. Deep AUC Maximization with Self-
Supervised Learning for COVID-19 Chest X-ray Classifications

Siyuan He National Research Council Canada Pengcheng Xi National Research
Council Canada Ashkan Ebadi National Research Council Canada Stéphane Tremblay
National Research Council Canada Alexander Wong University of Waterloo Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Effective representation learning is the key in improving model performance
for medical image analysis. In training deep learning models, a compromise
often must be made between performance and trust, both of which are essential
for medical applications. Moreover, models optimized with cross-entropy loss
tend to suffer from unwarranted overconfidence in the majority class and over-
cautiousness in the minority class. In this work, we integrate a new surrogate
loss with self-supervised learning for computer-aided screening of COVID-19
patients using radiography images. In addition, we adopt a new quantification
score to measure a model’s trustworthiness. Ablation study is conducted for
both the performance and the trust on feature learning methods and loss
functions. Comparisons show that leveraging the new surrogate loss on self-
supervised models can produce label-efficient networks that are both high-
performing and trustworthy.

##  1 Introduction

COVID-19 continues to affect our daily lives. In the fight against the
pandemic, computer-aided screening of patients using radiography images has
served as a complementary approach to standard polymerase chain reaction (PCR)
test. Despite our efforts in improving the performance of deep learning models
[1, 2], a compromise often must be made between the performance and model
trust [3].

Effective representation learning is the key in improving model performance.
The most common approach is through supervised learning on large-scale data
sets with labels. It can also be realized through un-supervised training, when
data labels are missing or expensive to collect. A series of self-supervised
models have achieved comparable performance to its supervised counterparts on
benchmark data sets [4, 5]. They learn image representations through
minimizing an embedding distance between image pairs derived from the same
image, while maximizing the distance between the pairs from different images.

Regarding model trust, a deep classification neural network optimized on
cross-entropy loss tends to be over-confident in its incorrect predictions of
the majority class, and overcautious in its correct predictions of the
minority class [6]. In this work, we investigate a new surrogate loss function
termed deep AUC Maximization [7]. We integrate it with self-supervised models
pre-trained with the MoCo framework [4] for improving model performance and
trust.

To our best knowledge, this is the first integration of the deep AUC
Maximization loss function with self-supervised learning. In addition, we
validate the models through quantitative comparisons to gain insight into the
models’ trustworthiness. Our assumption is that, by adopting the new surrogate
loss function to self-supervised models, we no longer need to sacrifice model
trust for performance but can achieve both. In summary, our contributions are
threefold:

  * •

We proposed the use of a new surrogate loss on self-supervised models to
improve representation learning and model trust pertinent to the task of
screening COVID-19 patients with deep learning;

  * •

We demonstrated with experimental results that the self-supervised models had
improved model performance over supervised ones on screening COVID-19
patients;

  * •

We showed that the use of a new surrogate loss can produce models that are
more trustworthy than those optimized with cross-entropy loss.

##  2 Literature Review

Self-supervised learning has gained momentum in learning visual
representations. It can be categorized into generative and discriminative
approaches. As a discriminative method, Momentum Contrast (MoCo) trains a
visual representation encoder by matching an encoded query to a dictionary of
encoded keys through a contrastive loss [4]. The query encoder is shared with
the key encoder, which receives slow updates in order to achieve consistency
in learning visual representations.

In medical AI, contrastive learning has led to improved representation
learning. In [8], a model named MoCo-CXR proved that linear models trained on
MoCo-CXR-pretrained representations outperform those without MoCo-CXR-
pretrained representations. Due to the scarcity of COVID-19 patient data, the
MoCo model has been applied to predicting patient deterioration based on chest
X-rays [9].

Area under the Receiver Operating Characteristic curve (AUC) is widely used in
medical image analysis for evaluating the performance of a neural network.
Recently, Yuan et al. [7] proposed a novel surrogate loss over the standard
cross-entropy loss to directly optimize for the AUC metric. AUC maximization,
as the authors claim, can lead to the largest increase in a network’s
performance. This new surrogate loss function was integrated with supervised
deep learning models and it achieved the first place in the Stanford CheXpert
competition [7].

##  3 Methodology

###  3.1 Model Architecture

Our model architecture is illustrated in Fig. 1. Our approach leverages deep
AUC maximization [7], a novel surrogate loss proposed for medical image
classification, with self-supervised pre-training to maximize label
efficiency, performance, and model trust. In our experiments, we compare the
loss function against traditional cross-entropy (CE) optimization on both
self-supervised and supervised models. The self-supervised model is built on
the MoCo framework [4] and it is pre-trained on MIMIC-CXR dataset [10]. All
models are then fine-tuned on the COVIDx8B dataset [11] for validation.
DenseNet-121 is chosen as the backbone architecture throughout our experiments
[12].

Figure 1: Model Architecture and Fine-tuning Pipeline.

###  3.2 Datasets

The MoCo model has been pre-trained on the MIMIC-CXR dataset for predicting
patient deterioration [9]. The dataset is composed of 377,110 chest
radiographs [10]. As the dataset was constructed before the COVID-19 pandemic,
it does not contain any positive chest X-ray samples of COVID-19.

Table 1: Data split for COVIDx8B Split | Negative | Positive | Total
---|---|---|---
Train | 13,794 | 2,158 | 15,952
Test | 200 | 200 | 400

We perform end-to-end fine tuning on the COVIDx8B dataset [11]. The latest
version of COVIDx8B consists of 15,952 chest radiographs for training and 400
for testing (Table 1). Each sample in the dataset is labelled as either
COVID-19 positive or negative. Stratified 5-fold cross validation is conducted
on the training split during the fine-tuning stage to evaluate model
performance.

###  3.3 Experiment Setups

The DenseNet-121 model pre-trained on MIMIC-CXR using the MoCo framework has a
projection dimension of 128, whereas the supervised model pre-trained on
ImageNet has a projection dimension of 1,000. For end-to-end fine tuning, the
parameters of the last fully connected layer of both pre-trained models are
replaced and randomly initialized with a single output neuron for binary
classification. We apply a sigmoid layer over the raw logits of the model to
obtain a probability distribution. All input images are resized to 224x224,
center cropped and normalized. Only random horizontal flipping is used for
data augmentation as further augmentations were noted to provide little
improvement for classification [9].

AUC Maximization. We adpot a novel surrogate loss function introduced by [7]
to maximize the area under the Receiver Operating Characteristic curve. For
end-to-end fine tuning, we use a learning rate of $0.1$ for all layers of the
DenseNet model. Then, we optimize the network with the new surrogate loss
function to maximize the AUC metric. Lastly, we train for 30 epochs while
decaying the learning rate by a factor of 10 at the 15th epoch.

CE Optimization. For standard end-to-end fine tuning, we set the learning rate
at $1e-3$ for all layers of the DenseNet model. Following similar procedures
to [9], we use cosine annealing learning rate decay to reduce the learning
rate. Finally, we use the SGD optimizer on cross-entropy loss with a momentum
of 0.9 and weight decay of $1e-4$ to fine tune the model for 30 epochs.

During each validation fold, we first compute an optimal threshold by
maximizing F1-score on the validation split. Then, we save the model
corresponding to the best validation accuracy. Finally, we evaluate the saved
models on the unseen test split.

##  4 Experimental Results

###  4.1 Supervised vs. Self-Supervised Pre-training

We first examine the performance difference between traditional supervised
pre-training on ImageNet and self-supervised contrastive pre-training on
MIMIC-CXR. Tables 2 and 3 show significant improvements in the precision
metric of the negative class and the sensitivity metric of the positive class
for both CE optimization and AUC maximization. In medical image analysis, this
improvement is key as maximizing the positive sensitivity score is necessary
to lower false-negatives.

However, this increase in performance comes at the cost of model trust. We
examine the trustworthiness of each model by calculating a trust score of the
positive class. As per the methods introduced in [3], we compute a score that
rewards well-placed confidence and penalizes undeserved overconfidence. In
Table 4, we notice that in the case of CE optimization, supervised models are
drastically more trustworthy than self-supervised models. Moreover, throughout
our CE optimization experiments, we observed that self-supervised models are
less confident in its correct predictions (overcautious) than its supervised
counterparts.

###  4.2 CE Optimization vs. AUC Maximization

Our comparisons of CE Optimization against AUC Maximization show improvements
across standard metrics as well as overall model trust-worthiness. Both Table
2 and Table 3 show improvements in the precision and sensitivity metrics
regarding the supervised models. Moreover, Fig. 2 demonstrates an increase in
the AUC scores of the supervised models. When examining self-supervised
models, AUC maximization still achieves performance comparable to CE
optimization.

More importantly, we observe significant gains in model trust scores,
especially in the context of self-supervised models. Table 4 shows a near 1%
increase in supervised pre-training and a near 6% increase in self-supervised
pre-training. Moreover, when using AUC maximization, we do not see the same
diminishment in model trust between supervised and self-supervised models.
Therefore, unlike CE optimization, we can freely leverage AUC maximization
with self-supervised pre-training to improve performance without sacrificing
model trust. As shown in Tables 2, 3 and 4, AUC maximization allows us to
achieve top metrics without trading off model trust for performance.

Table 2: Precision scores on the unseen COVIDx8B test split. The best metric out of each optimization method is bolded. The best metric across methods is denoted by *. Pre-trained Model | Negative | Positive
---|---|---
Supervised (CE Opt) | 0.8960 $\pm$ 1.6% | 0.9956 $\pm$ 0.4%
Self-Supervised (CE Opt) |  0.9295* $\pm$ 1.6% |  0.9978 $\pm$ 0.4%
Supervised (AUC Max) | 0.9134 $\pm$ 1.4% | 1.000
Self-Supervised (AUC Max) |  0.9251 $\pm$ 0.6% |  1.000*
Table 3: Sensitivity scores on the unseen COVIDx8B test split. The best metric out of each optimization method is bolded. The best metric across methods is denoted by *. Pre-trained Model | Negative | Positive
---|---|---
Supervised (CE Opt) | 0.9960 $\pm$ 0.3% | 0.8840 $\pm$ 2.1%
Self-Supervised (CE Opt) |  0.9980 $\pm$ 0.4% |  0.9240* $\pm$ 1.9%
Supervised (AUC Max) | 1.000 | 0.9050 $\pm$ 1.7%
Self-Supervised (AUC Max) |  1.000* |  0.9190 $\pm$ 0.7%
Table 4: Trust scores calculated from each experiment on the positive class. The best score overall is bolded. Cost Function | Supervised | Self-Supervised
---|---|---
CE Opt | 0.929 $\pm$ 0.9 % | 0.879 $\pm$ 1.4 %
AUC Max |  0.938 $\pm$ 1.0 % | 0.937 $\pm$ 1.0 %
Figure 2:  Area under the Receiver Operating Characteristic curve. Error bars
represent the standard deviation across cross-validation runs.

##  5 Conclusion

This work demonstrates that we no longer need to sacrifice model trust for
performance. Integrating AUC maximization can produce more trustworthy and
better performing models. By extending the AUC maximization paradigm [7] to
self-supervised pre-training, we showed that we can significantly improve key
metrics while also maintaining model trust.

We expect that our study on self-supervised learning with AUC maximization
will contribute to the classification of both COVID-19 and future illnesses.
More often than not, we cannot afford to collect large amount of labeled
samples at the onset of a pandemic. Therefore, it is important that we exploit
existing data, apply effective representation learning to maximizing model
performance, and gain optimal model confidence.

## Acknowledgments

We would like to thank Jianxing (Jason) Zhang for his preliminary work and
code preparation that led to the completion of this work. We also acknowledge
support from the Pandemic Response Challenge Program at the National Research
Council of Canada.

## References

  * [1] J. Zhang, P. Xi, A. Ebadi, H. Azimi, S. Tremblay, and A. Wong, “Covid-19 detection from chest x-ray images using imprinted weights approach,” _ICLR 2021 Workshop: Machine Learning for Preventing and Combating Pandemics_ , 2021.
  * [2] H. As’ ad, H. Azmi, P. Xi, A. Ebadi, S. Tremblay, and A. Wong, “Covid-19 detection from chest x-ray images using deep convolutional neural networks with weights imprinting approach,” _Journal of Computational Vision and Imaging Systems_ , vol. 6, no. 1, pp. 1–3, 2020.
  * [3] A. Wong, X. Y. Wang, and A. Hryniowski, “How much can we really trust you? towards simple, interpretable trust quantification metrics for deep neural networks,” 2021.
  * [4] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020\.
  * [5] T. Chen, S. Kornblith, M. Norouzi, and G. E. Hinton, “A simple framework for contrastive learning of visual representations,” 2020. [Online]. Available: https://arxiv.org/abs/2002.05709
  * [6] C. Esposito, G. A. Landrum, N. Schneider, N. Stiefl, and S. Riniker, “Ghost: Adjusting the decision threshold to handle imbalanced data in machine learning,” _Journal of Chemical Information and Modeling_ , vol. 61, no. 6, pp. 2623–2640, 2021, pMID: 34100609. [Online]. Available: https://doi.org/10.1021/acs.jcim.1c00160
  * [7] Z. Yuan, Y. Yan, M. Sonka, and T. Yang, “Large-scale robust deep auc maximization: A new surrogate loss and empirical studies on medical image classification,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021.
  * [8] H. Sowrirajan, J. Yang, A. Y. Ng, and P. Rajpurkar, “Moco pretraining improves representation and transferability of chest x-ray models,” in _Proceedings of the Fourth Conference on Medical Imaging with Deep Learning_ , ser. Proceedings of Machine Learning Research, M. Heinrich, Q. Dou, M. de Bruijne, J. Lellmann, A. Schläfer, and F. Ernst, Eds., vol. 143\. PMLR, 07–09 Jul 2021, pp. 728–744. [Online]. Available: https://proceedings.mlr.press/v143/sowrirajan21a.html
  * [9] A. Sriram, M. Muckley, K. Sinha, F. Shamout, J. Pineau, K. Geras, L. Azour, Y. Aphinyanaphongs, N. Yakubova, and W. Moore, “Covid-19 prognosis via self-supervised representation learning and multi-image prediction,” 2021.
  * [10] A. E. W. Johnson, T. J. Pollard, S. J. Berkowitz, N. R. Greenbaum, M. P. Lungren, C.-y. Deng, R. G. Mark, and S. Horng, “Mimic-cxr, a de-identified publicly available database of chest radiographs with free-text reports,” _Scientific Data_ , vol. 6, no. 1, p. 317, Dec 2019. [Online]. Available: https://doi.org/10.1038/s41597-019-0322-0
  * [11] L. Wang, Z. Q. Lin, and A. Wong, “Covid-net: a tailored deep convolutional neural network design for detection of covid-19 cases from chest x-ray images,” _Scientific Reports_ , vol. 10, no. 1, p. 19549, Nov 2020. [Online]. Available: https://doi.org/10.1038/s41598-020-76550-z
  * [12] G. Huang, Z. Liu, L. V. D. Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. Los Alamitos, CA, USA: IEEE Computer Society, jul 2017, pp. 2261–2269. [Online]. Available: https://doi.ieeecomputersociety.org/10.1109/CVPR.2017.243