© The Authors 2026

1. Introduction

The stability of magnetic fields in radiative stellar interiors is a long-standing problem in astrophysics. Observations show that a fraction of early-type stars, white dwarfs, and neutron stars harbor large-scale surface fields that can have a fossil origin (cf. Borra et al. 1982) or be the result of an unstable phase (Arlt & Rüdiger 2011). Asteroseismology has recently uncovered evidence for strong internal fields buried deep within the radiative cores of low-mass evolved stars. For instance, asteroseismic studies based on Kepler observations of red giants reveal suppressed dipole oscillation modes that imply magnetic fields of order 10⁵ G in the core (Fuller et al. 2015; Stello et al. 2016). Additionally, asteroseismology has revealed asymmetric splittings in the mixed modes frequency spectrum of low mass evolved stars that can be attributed to the presence of strong radial magnetic fields located deep inside the radiative region of the core (Li et al. 2023). Magnetohydrodynamic (MHD) instabilities can also induce angular momentum transport and mixing, potentially contributing to explain the slow rotation of red giant cores (see Aerts et al. 2019 for a review). Thus, understanding which field configurations can remain stable in a radiative zone, and the dynamics of unstable configurations, can have far-reaching implications for stellar evolution.

The stability of configurations containing a purely toroidal field has been extensively studied in the literature. A wide range of different instabilities have been identified, such as the shear-driven magnetorotational instability (Rüdiger et al. 2007, for toroidal fields) and the magnetic buoyancy instability (Parker 1966). Although these instabilities can have an impact on many astrophysical phenomena, in the context of the magnetism of radiative stellar cores, where stable stratification is dominant, it is generally accepted that the most prominent ones are those that can develop through almost horizontal displacements (Spruit 1999).

Tayler (1973) demonstrated that an axisymmetric toroidal field B_ϕ can develop a nonaxisymmetric kink instability, which grows on an Alfvén timescale and can occur for arbitrarily small values of the field strength in the absence of diffusive effects. This instability has later been identified in liquid metal experiments (Rüdiger et al. 2012; Seilmayer et al. 2012). Studies following up on Tayler’s work showed that rotation and stable stratification lower the growth rate but cannot suppress the instability entirely (Acheson & Gibbons 1978; Pitts & Tayler 1985; Spruit 1999; Kitchatinov & Rüdiger 2008; Bonanno & Urpin 2012; Skoutnev & Beloborodov 2024; Meduri et al. 2025). In the limit of strong stratification where the Brunt-Väisälä frequency dominates the toroidal Alfvén frequency, N ≫ ω_A, both global linear stability analyses and direct numerical simulations show that the growth rate scales as σ ∼ ω_A³/N² for perturbations with large radial wavelengths, which can be small but still nonzero, while small scale modes can continue to grow at the adiabatic growth rate ω_A = B_ϕ/(4πρ)^1/2R (Bonanno & Urpin 2012; Meduri et al. 2025).

Purely poloidal fields with field lines that close inside or outside the star can also be unstable (Markey & Tayler 1973; Flowers & Ruderman 1977). Purely toroidal and purely poloidal fields are therefore inherently unstable in stellar interiors. Mixed-field configurations are more likely to achieve stability, although analytic theory offers no simple general proof of global stability. Prendergast (1956) found a stable analytical mixed configuration in ideal MHD, while Duez & Mathis (2010) generalized this solution. However, recent work shows that the Prendergast solution can be destabilized by resistive effects, even under strongly stable stratification (Kaufman et al. 2022).

Direct numerical simulations show that a random initial field in a stably stratified star relaxes into a seemingly stable configuration with comparable toroidal and poloidal components (Braithwaite & Nordlund 2006). Such calculations provide valuable insight into the nonlinear evolution of candidate MHD equilibria. However, the absence of evident instabilities in the simulations does not necessarily imply stability, as it may result either from the damping of small-scale unstable modes by by explicit diffusivities much higher than those of real objects, or from insufficient numerical resolution due to numerical diffusivities. Our linear analysis aims at clarifying the stability conditions of a class of mixed poloidal-toroidal configurations and may guide future simulations.

Bonanno & Urpin (2011) showed that mixed configurations can support modes with high azimuthal order m, which remain unstable for any ratio of poloidal to toroidal field strength in the absence of viscous and magnetic diffusion. These modes satisfy local selection conditions and, due to their short azimuthal wavelengths, could be missed by numerical simulations. This type of instability was already identified in the context of plasma physics for a mixed axisymmetric magnetic configuration containing an azimuthal component B_ϕ and an axial one B_z. It occurs when the condition known as the Suydam criterion, which is necessary for stability,

$$\begin{array}{c}\hfill \frac{\mathit{dP}}{\mathit{ds}}+\frac{1}{8}s{B}_z^2[\frac{s{B}_z}{B_{\varphi }}\frac{d}{\mathit{ds}}(\frac{B_{\varphi }}{s{B}_z}){]}^2\ge 0\end{array}$$(1)

is violated (Suydam 1958). Here P is the local gas pressure and s represents the cylindrical radial coordinate. As we shall see in the following, the background field configuration chosen in the present study violates this criterion everywhere in the integration domain. The Suydam criterion neglects both stable stratification and diffusive effects.

In this work, we investigate this instability in the presence of stable stratification and all relevant diffusivities, both key ingredients of radiative stellar interiors that have not been considered before. Using a radially global linear stability analysis in spherical geometry, supported by local arguments and a semi-analytical approach in cylindrical coordinates, we show that simple mixed poloidal-toroidal equilibria can be subject to a type of selective instability even under strong stable stratification. Unlike the classical Tayler modes, which are dominated by unstable perturbations with azimuthal wavenumber m = 1, this instability can develop even for fields with a substantial poloidal component and occur for any azimuthal wavenumber, provided that the perturbation wavevector k satisfies k ⋅ B ≈ 0, where B is the magnetic field.

The rest of this paper is organized as follows. In Sect. 2 we describe the basic state adopted for our global linear analysis and the governing equations. Section 3 presents the results of our linear stability analysis in a spherical geometry. First, we discuss the diffusionless case with no gravity, then we introduce stable stratification and thermal diffusion, and finally we incorporate viscosity and magnetic diffusivity. Sect. 4 focuses on the region close to the axis, analyzing a complementary problem formulated in cylindrical coordinates and we compare the results obtained with the spherical case. Section 5 closes the paper with a summary of the results and discussion.

2. Equilibrium model and governing equations

We study an MHD basic state describing a fluid of constant density with a the background pressure P chosen to satisfy the magneto-hydrostatic equilibrium

$$\begin{array}{c}\hfill -\frac{\mathbf{\nabla }P}{\rho }+\mathit{g}+\frac{1}{4\pi \rho }(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}=\mathbf{0}.\end{array}$$(2)

Here, P = p(r)+h(r, θ) is the sum of a spherically symmetric component p(r), which balances gravity and depends only on the spherical radius r, and a nonspherically symmetric contribution needed to balance the Lorentz force h(r, θ), which also depends on the colatitude θ. ρ is the constant density and g is gravity which, for a gas sphere of uniform density, increases linearly with radius.

We study the stability of mixed magnetic field configurations where both poloidal and toroidal components are present. We linearize the MHD equations around the background state above by separating the unknowns into a mean and a fluctuating component, the latter denoted by a prime symbol.

We employ the Boussinesq approximation, neglecting all density perturbations except the ones related to the stabilizing buoyancy. The equation of state is dρ/ρ = −α dT/T, where α is the thermal expansion coefficient. By keeping only the linear terms in the perturbed quantities, we obtain

$$\begin{array}{cc}& \frac{\partial \mathit{u}\mathbf{\prime }}{\partial t}=-\frac{1}{\rho }\mathbf{\nabla }P\prime +\frac{1}{4\pi \rho }(\mathbf{\nabla }\times \mathit{B}\mathbf{\prime })\times \mathit{B}\hfill \\ \hfill & +\frac{1}{4\pi \rho }(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}\mathbf{\prime }+\frac{\rho \prime }{\rho }\mathit{g}+\nu {\mathrm{\nabla }}^2\mathit{u}\mathbf{\prime }\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{B}\mathbf{\prime }}{\partial t}=\mathbf{\nabla }\times (\mathit{u}\mathbf{\prime }\times \mathit{B})+\eta {\mathrm{\nabla }}^2\mathit{B}\mathbf{\prime }\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{\partial T\prime }{\partial t}+\mathit{u}\mathbf{\prime }·(\mathbf{\nabla }-{\mathbf{\nabla }}_{\mathrm{ad}})T=\kappa {\mathrm{\nabla }}^{2}T\prime \hfill \end{array}$$(5)

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{u}\mathbf{\prime }=0,\mathbf{\nabla }·\mathit{B}\mathbf{\prime }=0,\hfill \end{array}$$(6)

where the kinematic viscosity ν, magnetic diffusivity η, and thermal conductivity κ are constant. ∇_ad is the adiabatic temperature gradient.

We solve these equations following two different approaches. In the first one, we carry out a stability analysis in spherical coordinates (r, θ, ϕ), global in the radial direction and local in the latitudinal and azimuthal ones, including stable stratification and all of the diffusive effects. The perturbed quantities are assumed to be proportional to exp(σt − ilθ − imϕ). In this approximation, the perturbations vary on small scales in the latitudinal direction compared to the background field gradients, i.e. k_θ ≫ H_θ⁻¹, where k_θ = l/r and H_θ = r∂lnB/∂θ.

Our second approach adopts cylindrical coordinates (s, ϕ, z) and expands on Bonanno & Urpin (2011) by including stable stratification along the vertical direction z and thermal diffusion. The analysis is global in the cylindrical radius s and local in the vertical and azimuthal directions, and the perturbations are assumed to be proportional to exp(σt − ik_zz − imϕ). Since the background state has cylindrical symmetry, there is no condition between the background and perturbed quantities to be satisfied as in the spherical case.

3. Spherical coordinates

We first consider spherical coordinates (r, θ, ϕ) and a mixed poloidal-toroidal background magnetic field configuration of the form

$$\begin{array}{c}\hfill \mathit{B}=B_{\varphi }(r,\theta ){\widehat{\mathit{e}}}_{\varphi }+B_{z}cos\theta {\widehat{\mathit{e}}}_r-B_{z}sin\theta {\widehat{\mathit{e}}}_{\theta },\end{array}$$(7)

which satisfies the magnetohydrostatic equilibrium in Eq. (2). The toroidal field B_ϕ is axisymmetric,

$$B_{\varphi }=B_{\varphi 0}f(r,\theta ),$$

where f is a function of order unity which we define below, while B_z is the (constant) axial field strength. We express the perturbed quantities u′ and B′ in terms of toroidal and poloidal potentials as

$$\begin{array}{cc}\hfill {\mathit{u}}^{\prime }& =\mathbf{\nabla }\times (\mathbf{\nabla }\times W{\widehat{\mathit{e}}}_r)+\mathbf{\nabla }\times Z{\widehat{\mathit{e}}}_r\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {\mathit{B}}^{\prime }& =\mathbf{\nabla }\times (\mathbf{\nabla }\times \mathrm{\Phi }{\widehat{\mathit{e}}}_r)+\mathbf{\nabla }\times \mathrm{\Psi }{\widehat{\mathit{e}}}_r\hfill \end{array}$$(9)

in order to reduce the number of unknowns and have solutions that are divergence free by construction.

The equation of evolution of the poloidal field potential Φ is obtained from the radial component of Eq. (4), and the equations for the toroidal flow potential Z and the toroidal field potential Ψ are obtained from the radial components of the curled Eqs. (3) and (4), respectively. Lastly, the radial component of Eq. (3) curled twice gives the equation for the poloidal flow W.

The system comprises four coupled equations, whose coefficients involve terms of various orders in the product l m. For each equation, and for each derivative order, we retain only the highest-order terms in l m appearing in the corresponding coefficient. Therefore, we do not keep unnecessary terms under the l ≫ 1 approximation, but we also do not filter out possible solutions with m ≫ 1. Under these assumptions, we obtain an equation of evolution for the poloidal flow potential

$$\begin{array}{cc}& \sigma k_{\perp }^{4}W-\sigma k_{\perp }^2\frac{d^{2}W}{d{r}^2}=-k_{\perp }^2\nu \frac{d^{4}W}{d{r}^4}\hfill \\ \hfill & +2k_{\perp }^2\nu \frac{d^{2}W}{d{r}^2}-\frac{4k_{\perp }^4\nu }{r}\frac{\mathit{dW}}{\mathit{dr}}-k_{\perp }^6\nu W\hfill \\ \hfill & -\frac{1}{4\pi \rho }\{k_{\perp }^2{B}_{z}cos\theta \frac{d^3\mathrm{\Phi }}{d{r}^3}+k_{\perp }^4\frac{m{B}_{\varphi }csc\theta -l{B}_{z}sin\theta }{r}\frac{d^2\mathrm{\Phi }}{d{r}^2}\hfill \\ \hfill & +[k_{\perp }^4{B}_{z}cos\theta +\frac{2lmcsc\theta }{r^4}(B_{\varphi }cot\theta -\frac{d{B}_{\varphi }}{d\theta }+rcot\theta \frac{d{B}_{\varphi }}{\mathit{dr}}\hfill \\ \hfill & -r\frac{d^2{B}_{\varphi }}{drd\theta })]\frac{d\mathrm{\Phi }}{\mathit{dr}}-k_{\perp }^4\frac{2B_{z}cos\theta +i(m{B}_{\varphi }csc\theta -l{B}_{z}sin\theta )}{r}\mathrm{\Phi }\hfill \\ \hfill & +\frac{im{B}_z}{r^2}\frac{d^2\mathrm{\Psi }}{d{r}^2}-\frac{csc\theta }{r^3}(2l^2{B}_{\varphi }cos\theta +2m^{2}csc\theta \frac{d{B}_{\varphi }}{d\theta })\frac{d\mathrm{\Psi }}{\mathit{dr}}\hfill \\ \hfill & -i{k}_{\perp }^2\frac{2l{B}_{\varphi }+m{B}_{z}cos\theta }{r^2}\mathrm{\Psi }\}+g\alpha k_{\perp }^2\frac{T\prime }{T},\hfill \end{array}$$(10)

the toroidal flow potential

$$\begin{array}{cc}& \sigma k_{\perp }^{2}Z=k_{\perp }^2\nu \frac{d^{2}Z}{d{r}^2}-k_{\perp }^4\nu Z\hfill \\ \hfill & +\frac{1}{4\pi \rho }\{\frac{\mathit{im}}{r^2}{B}_{z}sin\theta \frac{d^2\mathrm{\Phi }}{d{r}^2}-\frac{k_{\perp }^2}{r}(B_{\varphi }cot\theta -\frac{d{B}_{\varphi }}{d\theta })\frac{d\mathrm{\Phi }}{\mathit{dr}}\hfill \\ \hfill & -\frac{i{k}_{\perp }^2}{r^2}[l(B_{\varphi }+r\frac{d{B}_{\varphi }}{\mathit{dr}})+m{B}_z]\mathrm{\Phi }+k_{\perp }^2{B}_{z}cos\theta \frac{d\mathrm{\Psi }}{\mathit{dr}}\hfill \\ \hfill & +\frac{i{k}_{\perp }^2}{r}(l{B}_{z}sin\theta -m{B}_{\varphi }csc\theta )\mathrm{\Psi }\},\hfill \end{array}$$(11)

the poloidal field potential

$$\begin{array}{cc}\hfill \sigma k_{\perp }^2\mathrm{\Phi }& =k_{\perp }^2{B}_{z}cos\theta \frac{\mathit{dW}}{\mathit{dr}}-\frac{i{k}_{\perp }^2}{r}(m{B}_{\varphi }csc\theta -l{B}_{z}sin\theta )W\hfill \\ \hfill & -\frac{im{B}_z}{r^2}Z+k_{\perp }^2\eta \frac{d^2\mathrm{\Phi }}{d{r}^2}-k_{\perp }^4\eta \mathrm{\Phi },\hfill \end{array}$$(12)

the toroidal field potential

$$\begin{array}{cc}\hfill \sigma k_{\perp }^2\mathrm{\Psi }& =-\frac{im{B}_z}{r^2}\frac{d^{2}W}{d{r}^2}\hfill \\ \hfill & -\frac{l^2-m^{2}csc\theta }{r^3}(B_{\varphi }cot\theta -\frac{d{B}_{\varphi }}{d\theta })\frac{\mathit{dW}}{\mathit{dr}}\hfill \\ \hfill & -\frac{i{k}_{\perp }^2}{r^2}[l(B_{\varphi }-r\frac{d{B}_{\varphi }}{\mathit{dr}})+m{B}_z]W+k_{\perp }^2{B}_{z}cos\theta \frac{\mathit{dZ}}{\mathit{dr}}\hfill \\ \hfill & -\frac{i{k}_{\perp }^2}{r}(m{B}_{\varphi }csc\theta -l{B}_{z}sin\theta )Z+k_{\perp }^2\eta \frac{d^2\mathrm{\Psi }}{d{r}^2}-k_{\perp }^4\eta \mathrm{\Psi },\hfill \end{array}$$(13)

and the temperature perturbations

$$\begin{array}{c}\hfill \sigma T\prime =-\frac{\kappa }{r^2}\frac{d}{\mathit{dr}}(r^2\frac{dT\prime }{\mathit{dr}})-\kappa k_{\perp }^{2}T\prime -\frac{N^{2}T}{\alpha g_{\mathrm{in}}}{k}_{\perp }^{2}W,\end{array}$$(14)

where g_in is the acceleration of gravity calculated at the inner boundary r = r_in. The perpendicular wavenumber k_⊥ is k_⊥² = l²/r² + m²/r²sin²θ and the Brunt-Väiälä frequency is defined by

$$N^2=g_{\mathrm{in}}\frac{\alpha }{T}(\mathbf{\nabla }-{\mathbf{\nabla }}_{\text{ad}})T·{\widehat{\mathit{e}}}_r.$$

At both boundaries, the perturbations are imposed to be zero, which yields

$$\begin{array}{c}\hfill W=\frac{\mathit{dW}}{\mathit{dr}}=0,Z=0,\mathrm{\Phi }=\frac{d\mathrm{\Phi }}{\mathit{dr}}=0,\mathrm{\Psi }=0,T=0\end{array}$$(15)

at r = r_in and r = r_out.

To nondimensionalize the equations above, we scale time in units of the inverse of the toroidal Alfvén frequency ω_A0 = B_ϕ0/(4πρ)^1/2r_in, lengths in units of the inner radius r_in, and magnetic field intensities in units of the background toroidal field strength B_ϕ0. The background temperature gradient is $(\mathbf{\nabla }-{\mathbf{\nabla }}_{\text{ad}})T\equiv A(r){\widehat{\mathit{e}}}_r$ and the temperature perturbations are nondimensionalized with A(r_in)r_in. In this scaling scheme, there are five key dimensionless parameters. The first is the ratio of the Brunt-Väisälä frequency at the inner boundary, N₀ = N(r_in), to the Alfvén frequency ω_A0,

$$\delta =\frac{N_0}{{\omega }_{\mathrm{A}0}}$$

and quantifies the strength of stable stratification. The second is the ratio of the axial to the toroidal field

$$\mu =\frac{B_z}{B_{\varphi }}·$$

The remaining parameters, measuring the relative importance of diffusive effects, are: the ratio of the thermal diffusion rate κ/r_in² to the Alfvén frequency,

$$\begin{array}{c}\hfill ϵ=\frac{\kappa }{r_{\mathrm{in}}^2{\omega }_{A0}},\end{array}$$

the Lundquist number,

$$\begin{array}{c}\hfill Lu=\frac{{\omega }_{A0}{r}_{\mathrm{in}}^2}{\eta },\end{array}$$

and the magnetic Prandtl number

$$Pm=\frac{\nu }{\eta }·$$

Nondimensional variables are denoted with a tilde hereafter and the system of dimensionless equations is

$$\begin{array}{cc}& \mathrm{\Gamma }{\stackrel{\sim }{k}}_{\perp }^4\stackrel{\sim }{W}-\mathrm{\Gamma }{\stackrel{\sim }{k}}_{\perp }^2\frac{d^2\stackrel{\sim }{W}}{d{\stackrel{\sim }{r}}^2}=-{\stackrel{\sim }{k}}_{\perp }^2\frac{\mathit{Pm}}{\mathit{Lu}}\frac{d^4\stackrel{\sim }{W}}{d{\stackrel{\sim }{r}}^4}\hfill \\ \hfill & +2{\stackrel{\sim }{k}}_{\perp }^2\frac{\mathit{Pm}}{\mathit{Lu}}\frac{d^2\stackrel{\sim }{W}}{d{\stackrel{\sim }{r}}^2}-\frac{4{\stackrel{\sim }{k}}_{\perp }^4}{\stackrel{\sim }{r}}\frac{\mathit{Pm}}{\mathit{Lu}}\frac{d\stackrel{\sim }{W}}{d\stackrel{\sim }{r}}-{\stackrel{\sim }{k}}_{\perp }^6\frac{\mathit{Pm}}{\mathit{Lu}}\stackrel{\sim }{W}\hfill \\ \hfill & -{\stackrel{\sim }{k}}_{\perp }^2{\mu }_{0}cos\theta \frac{d^3\stackrel{\sim }{\mathrm{\Phi }}}{d{\stackrel{\sim }{r}}^3}+{\stackrel{\sim }{k}}_{\perp }^4\frac{mfcsc\theta -l{\mu }_{0}sin\theta }{\stackrel{\sim }{r}}\frac{d^2\stackrel{\sim }{\mathrm{\Phi }}}{d{\stackrel{\sim }{r}}^2}\hfill \\ \hfill & +[{\stackrel{\sim }{k}}_{\perp }^4{\mu }_{0}cos\theta +\frac{2lmcsc\theta }{{\stackrel{\sim }{r}}^4}(fcot\theta -\frac{\mathit{df}}{d\theta }+\stackrel{\sim }{r}cot\theta \frac{\mathit{df}}{d\stackrel{\sim }{r}}\hfill \\ \hfill & -\stackrel{\sim }{r}\frac{d^{2}f}{d\stackrel{\sim }{r}d\theta })]\frac{d\stackrel{\sim }{\mathrm{\Phi }}}{d\stackrel{\sim }{r}}-{\stackrel{\sim }{k}}_{\perp }^4\frac{2{\mu }_{0}cos\theta +i(mfcsc\theta -l{\mu }_{0}sin\theta )}{\stackrel{\sim }{r}}\stackrel{\sim }{\mathrm{\Phi }}\hfill \\ \hfill & +\frac{im{\mu }_0}{{\stackrel{\sim }{r}}^2}\frac{d^2\stackrel{\sim }{\mathrm{\Psi }}}{d{\stackrel{\sim }{r}}^2}-\frac{csc\theta }{{\stackrel{\sim }{r}}^3}(2l^{2}fcos\theta +2m^{2}csc\theta \frac{\mathit{df}}{d\theta })\frac{d\stackrel{\sim }{\mathrm{\Psi }}}{d\stackrel{\sim }{r}}\hfill \\ \hfill & -i{\stackrel{\sim }{k}}_{\perp }^2\frac{2lf+m{\mu }_{0}cos\theta }{{\stackrel{\sim }{r}}^2}\stackrel{\sim }{\mathrm{\Psi }}+{\delta }^2{\stackrel{\sim }{k}}_{\perp }^2\stackrel{\sim }{g}\stackrel{\sim }{T}\prime ,\hfill \\ \hfill & \mathrm{\Gamma }{\stackrel{\sim }{k}}_{\perp }^2\stackrel{\sim }{Z}={\stackrel{\sim }{k}}_{\perp }^2\frac{\mathit{Pm}}{\mathit{Lu}}\frac{d^2\stackrel{\sim }{Z}}{d{\stackrel{\sim }{r}}^2}-{\stackrel{\sim }{k}}_{\perp }^4\frac{\mathit{Pm}}{\mathit{Lu}}\stackrel{\sim }{Z}\hfill \\ \hfill & +\frac{\mathit{im}}{{\stackrel{\sim }{r}}^2}{\mu }_{0}sin\theta \frac{d^2\stackrel{\sim }{\mathrm{\Phi }}}{d{\stackrel{\sim }{r}}^2}-\frac{{\stackrel{\sim }{k}}_{\perp }^2}{\stackrel{\sim }{r}}(fcot\theta -\frac{\mathit{df}}{d\theta })\frac{d\stackrel{\sim }{\mathrm{\Phi }}}{d\stackrel{\sim }{r}}\hfill \\ \hfill & -\frac{i{\stackrel{\sim }{k}}_{\perp }^2}{{\stackrel{\sim }{r}}^2}[l(f+\stackrel{\sim }{r}\frac{\mathit{df}}{d\stackrel{\sim }{r}})+m{\mu }_0]\stackrel{\sim }{\mathrm{\Phi }}+{\stackrel{\sim }{k}}_{\perp }^2{\mu }_{0}cos\theta \frac{d\stackrel{\sim }{\mathrm{\Psi }}}{d\stackrel{\sim }{r}}\hfill \end{array}$$(16)

$$\begin{array}{cc}& +\frac{i{\stackrel{\sim }{k}}_{\perp }^2}{\stackrel{\sim }{r}}(l{\mu }_{0}sin\theta -mfcsc\theta )\stackrel{\sim }{\mathrm{\Psi }},\hfill \\ \hfill & \mathrm{\Gamma }{\stackrel{\sim }{k}}_{\perp }^2\stackrel{\sim }{\mathrm{\Phi }}={\stackrel{\sim }{k}}_{\perp }^2{\mu }_{0}cos\theta \frac{d\stackrel{\sim }{W}}{d\stackrel{\sim }{r}}-\frac{i{\stackrel{\sim }{k}}_{\perp }^2}{\stackrel{\sim }{r}}(mfcsc\theta -l{\mu }_{0}sin\theta )\stackrel{\sim }{W}\hfill \\ \hfill & -\frac{im{\mu }_0}{{\stackrel{\sim }{r}}^2}\stackrel{\sim }{Z}+{\stackrel{\sim }{k}}_{\perp }^2\frac{1}{\mathit{Lu}}\frac{d^2\stackrel{\sim }{\mathrm{\Phi }}}{d{\stackrel{\sim }{r}}^2}-{\stackrel{\sim }{k}}_{\perp }^4\frac{1}{\mathit{Lu}}\stackrel{\sim }{\mathrm{\Phi }},\hfill \end{array}$$(17)

$$\begin{array}{cc}& \mathrm{\Gamma }{\stackrel{\sim }{k}}_{\perp }^2\stackrel{\sim }{\mathrm{\Psi }}=-\frac{im{\mu }_0}{{\stackrel{\sim }{r}}^2}\frac{d^2\stackrel{\sim }{W}}{d{\stackrel{\sim }{r}}^2}\hfill \\ \hfill & -\frac{l^2-m^{2}csc\theta }{{\stackrel{\sim }{r}}^3}(fcot\theta -\frac{\mathit{df}}{d\theta })\frac{d\stackrel{\sim }{W}}{d\stackrel{\sim }{r}}\hfill \\ \hfill & -\frac{i{\stackrel{\sim }{k}}_{\perp }^2}{{\stackrel{\sim }{r}}^2}[l(f-\stackrel{\sim }{r}\frac{\mathit{df}}{d\stackrel{\sim }{r}})+m{\mu }_0]\stackrel{\sim }{W}+{\stackrel{\sim }{k}}_{\perp }^2{\mu }_{0}cos\theta \frac{d\stackrel{\sim }{Z}}{d\stackrel{\sim }{r}}\hfill \end{array}$$(18)

$$\begin{array}{cc}& -\frac{i{\stackrel{\sim }{k}}_{\perp }^2}{\stackrel{\sim }{r}}(mfcsc\theta -l{\mu }_{0}sin\theta )\stackrel{\sim }{Z}+{\stackrel{\sim }{k}}_{\perp }^2\frac{1}{\mathit{Lu}}\frac{d^2\stackrel{\sim }{\mathrm{\Psi }}}{d{\stackrel{\sim }{r}}^2}-{\stackrel{\sim }{k}}_{\perp }^4\frac{1}{\mathit{Lu}}\stackrel{\sim }{\mathrm{\Psi }},\hfill \end{array}$$(19)

$$\begin{array}{cc}& \mathrm{\Gamma }\stackrel{\sim }{T}\prime =-\frac{ϵ}{{\stackrel{\sim }{r}}^2}\frac{d}{d\stackrel{\sim }{r}}({\stackrel{\sim }{r}}^2\frac{d\stackrel{\sim }{T}\prime }{d\stackrel{\sim }{r}})-ϵ{\stackrel{\sim }{k}}_{\perp }^2\stackrel{\sim }{T}\prime -{\stackrel{\sim }{N}}^2{\stackrel{\sim }{k}}_{\perp }^2\stackrel{\sim }{W},\hfill \end{array}$$(20)

where

$$\stackrel{\sim }{g}=g/g_{\mathrm{in}},{\stackrel{\sim }{N}}^2=N^{2}T/\alpha g_{\mathrm{in}}A,{\mu }_0=B_z/B_{\varphi 0}.$$

The above equations, together with their boundary conditions, define an eigenvalue problem that we solved numerically using an eighth-order finite difference scheme for the radial derivatives. The radial domain is $1\le \stackrel{\sim }{r}=r/r_{\text{in}}\le 2$. The eigenvalues of the resulting finite matrix are evaluated numerically through the software MATLAB, employing a Krylov-Schur algorithm (Stewart 2001) to solve the eigensystem. Gravity varies linearly with radius, g = g_inr/r_in, and we employ a temperature gradient dT/dr ∝ 1/r², which satisfies a purely conductive heat equation.

The background toroidal field varies spatially following

$$f(r,\theta )=s/r_{\text{in}}=(rsin\theta )/r_{\text{in}}.$$

The mixed background field adopted here is not meant to describe a realistic global fossil-field equilibrium across the star, but rather to provide a minimal local representation of a regular mixed-field geometry near the symmetry axis. Regularity near the axis requires the toroidal component of an axisymmetric field to vanish at most linearly with cylindrical radius, and a smooth poloidal component is approximately uniform to leading order over a sufficiently small region. Our approach is complementary to direct numerical simulations that explore stable relaxed states. Such simulations cannot exclude the presence of locally unstable configurations at high azimuthal wavenumbers. A further motivation for this choice of background field is that the instability is found to grow preferentially near the axis.

For our background field configuration, B_ϕ ∝ s and B_z = const, the second term in Eq. (1) vanishes, so the Suydam stability condition reduces to dP/ds ≥ 0. Since the inward toroidal hoop stress implies dP/ds < 0 everywhere for magnetostatic balance, the present equilibrium configuration lies on the unstable side of the Suydam criterion.

To enable comparison with previous work, we first considered a purely toroidal field (μ₀ = 0), for which we recover the linear results reported by Meduri et al. (2025). In the limit of very short radial wavelengths, k_r ≫ k_θ, we also recover the fully local dispersion relation of Skoutnev & Beloborodov (2024) as expected. Introducing even a weak axial field changes the stability properties substantially. The configuration becomes unstable for all values of the azimuthal wavenumber m, but with vertical wavenumbers satisfying

$$\begin{array}{c}\hfill \mathit{k}·\mathit{B}\approx 0,\end{array}$$(21)

as we shall see in the following.

In the remainder of this section we investigate the properties of the instability first in the diffusionless limit and with no gravity, demonstrating the spectrally selective nature of the instability for μ₀ = 0.1, 0.2, and 1. Then, we consider stable stratification and vary δ from 0 to 1000, again for selected values of μ₀. For each of these configurations, we focus first on the polar regions and then on lower latitudes. We subsequently examine the effects of viscosity and magnetic diffusivity, for a fixed Lundquist number of Lu = 10⁴ and for no viscosity (Pm = 0) and a finite viscosity of Pm = 10⁻³. Finally, we incorporate all physical effects and derive the instability boundaries exploring the ranges 10² ≤ Lu ≤ 10⁶ and 0 ≤ δ ≤ 1000.

3.1. Diffusionless case with no gravity

We first consider the diffusionless case with no gravity (δ = 0) and discuss the properties of the instability for three values of the poloidal-to-toroidal field ratio, μ₀ = 0.1, 0.2, and 1. We calculate the unstable eigenmodes and the corresponding growth rates for different azimuthal wavenumbers m, first at high latitudes and then at the equator.

Close to the axis, where k_r ≈ k_z, the selection condition (21) reads

$$\begin{array}{c}\hfill k_r{B}_z+\frac{m}{rsin\theta }{B}_{\varphi }\approx 0,\end{array}$$(22)

which, for our choice of the toroidal field profile, yields

$$\begin{array}{c}\hfill {\stackrel{\sim }{k}}_r=k_r{r}_{\text{in}}\approx \frac{m}{{\mu }_0}·\end{array}$$(23)

The radial scale of the instability is then expected to decrease with m and to increase with the magnetic field ratio. Figure 1 shows the most unstable radial velocity eigenmodes at a colatitude of θ = 10° for μ₀ = 0.1 (left panels) and μ₀ = 0.2 (right panels) for different azimuthal wavenumbers m as indicated in the legend insets. The eigenmodes were computed for a latitudinal wavenumber l = 10. Choosing a different value of l has little impact on the results, as expected from the selection condition above, which is independent of this parameter. The growth rates of all these eigenmodes are of the order of the toroidal Alfvén frequency of the background field. For both values of the field ratio μ₀, the number of eigenmode’s nodes increases linearly with increasing m, as expected from Eq. (23). For a fixed m, the eigenmode radial lengthscale increases when increasing μ₀ from 0.1 to 0.2 as expected. A Fourier analysis of the eigenmodes shows that the power spectrum is sharply peaked at the radial wavenumber expected from the condition (23).

Fig. 1. Most unstable eigenmodes for the radial velocity perturbations ${\stackrel{\sim }{u}}_r\prime$ in the unstratified case (δ = 0) at a colatitude of θ = 10° for l = 10 and various azimuthal wavenumbers m. The poloidal-to-toroidal field ratio is μ0 = 0.1 and μ0 = 0.2 in the left and right panels, respectively.

At lower latitudes, the properties of the instability also comply well with predictions obtained from the selection condition (21). At the equator (θ = 90°), this condition reads

$$\begin{array}{c}\hfill -\frac{l}{r}{B}_z+\frac{m}{r}{B}_{\varphi }\approx 0,\end{array}$$(24)

which yields

$$\begin{array}{c}\hfill l\approx \frac{m{B}_{\varphi }}{B_z}=\frac{m\stackrel{\sim }{r}}{{\mu }_0}·\end{array}$$(25)

Figure 2 presents the nondimensional growth rates Γ = σ/ω_A0 at the equator as a function of the latitudinal wavenumber l in the case μ₀ = 0.1 and for different azimuthal wavenumbers (m = 1, 2, 4, and 10). The spectrally selective nature of the instability is evident from the localization of the unstable latitudinal wavenumbers l in a range determined by Eq. (25), where the eigenmodes grow at a rate of the order of the toroidal Alfvén frequency ω_A0. For m = 10 (dotted line), for example, we find that the spectrum of unstable modes roughly ranges between 100 and 200, with the selection condition above predicting 100 ≲ l ≲ 200, when using the domain range $1\le \stackrel{\sim }{r}\le 2$.

Fig. 2. Nondimensional growth rate Γ = σ/ωA0 as a function of the latitudinal wavenumber l at the equator for different azimuthal wavenumbers m. The other parameters are δ = 0 and μ0 = 0.1.

The structure of the spectrum for the azimuthal modes m = 1 and m = 2 shows two peaks (solid and dashed lines in Fig. 2), similarly to what reported by Bonanno & Urpin (2011). As we increase the azimuthal wavenumber, the two peaks tend to blend into each other, leading to a flat spectrum. This is an effect related to our choice of the basis functions. Increasing m leads to unstable modes that correspond to higher l. Since k_θ = l/r, both peaks of the instability can eventually fit in the integration domain for high enough values of l, and we have $d{k}_{\theta }\sim \frac{l}{r^2}dr$, with dk_θ being the physical width of the instability. The corresponding region of the integration domain that can be unstable, dr, or, equivalently, the separation in space between the two peaks, decreases as we increase l.

Figure 3 presents the most unstable eigenmodes ${\stackrel{\sim }{u}}_r\prime$ for different azimuthal wavenumbers at the equator and for μ₀ = 1. Each eigenmode has a latitudinal wavenumber l satisfying the selection condition (25) for $\stackrel{\sim }{r}=1.5$, that is for m = 10, 20, 40, 100, and 200 we fixed l to 15, 30, 60, 150, and 300 respectively. For this choice of l and μ₀, the two-peaks structure is absent in the spectrum and is instead incorporated in the eigenmodes, which are localized around $\stackrel{\sim }{r}=1.5$ as expected. Additionally, as we increase m, the eigenmodes become progressively more peaked around this radius for the reasons mentioned above.

Fig. 3. Eigenmodes at the equator for δ = 0 and μ0 = 1. The azimuthal wavenumbers m = 10, 20, 40, 100, and 200 are shown, with corresponding latitudinal wavenumber l = 15, 30, 60, 150, and 300 derived from the selection condition (25) at $\stackrel{\sim }{r}=1.5$.

3.2. Effect of stable stratification

To investigate how the instability is affected by the stabilizing buoyancy, we fixed ϵ = 0.01 and varied the stratification parameter δ from 10 to 1000. Viscosity and magnetic diffusion were neglected. In the following, we discuss the unstable eigenmodes and the corresponding growth rates only for μ₀ = 0.1 and 0.2, first in the polar region and subsequently at the equator, following the approach of the previous section. For these values of μ₀, we find that gravity tames the unstable modes in a way which is similar to the one of the classical Tayler instability of a purely toroidal field, as we shall see in the following.

Figure 4 shows the growth rate Γ as a function of δ for a latitudinal wavenumber l = 10 and various azimuthal modes m, evaluated at a colatitude of θ = 10°. The cases μ₀ = 0.1 and μ₀ = 0.2 are presented in the top and bottom panels, respectively. The unstable modes with azimuthal wavenumbers m ≥ 100 have small radial scales as a consequence of Eq. (23). These modes therefore develop through almost horizontal displacements that resist the stabilizing buoyancy and can grow at the adiabatic rate for almost all values of δ explored (green and cyan lines).

Fig. 4. Nondimensional growth rate Γ = σ/ωA0 as a function of the stratification parameter δ for l = 10 at colatitude θ = 10° for two values of μ0 = 0.1 (top panel) and 0.2 (bottom panel), and for various azimuthal wavenumbers m. The dashed black line indicates the scaling Γ ∝ δ−2.

Lower azimuthal wavenumbers are associated with unstable modes with larger radial length scales that can be tamed by gravity. An order of magnitude estimate for the minimum unstable radial wavenumber presenting adiabatic growth, k_c, can be obtained by equating the timescale for the buoyancy response t_N = κk_r⁴/k_⊥²N² to the timescale of adiabatic growth 1/σ_ad, where σ_ad is the growth rate in the absence of gravity, and of viscous and resistive effects (Skoutnev & Beloborodov 2024). In dimensionless form, this critical wavenumber is ${\stackrel{\sim }{k}}_c=k_c{r}_{\text{in}}={({\stackrel{\sim }{k}}_{\perp }\delta )}^{1/2}/{ϵ}^{1/4}{\mathrm{\Gamma }}_{\text{ad}}^{1/4}$, where Γ_ad = σ_ad/ω_A0. For the moderate values of the poloidal-to-toroidal field ratio μ₀ considered here, Γ_ad ≈ 1. However, as we shall see in Sect. 4, the adiabatic growth rate can be small for a field configuration dominated by the poloidal component.

Using the selection condition (23) with ${\stackrel{\sim }{k}}_r={\stackrel{\sim }{k}}_c$ provides a critical value of the stratification parameter,

$$\begin{array}{c}\hfill {\delta }_c={\mathrm{\Gamma }}_{\text{ad}}^{1/2}{ϵ}^{1/2}{m}^2/{\stackrel{\sim }{k}}_{\perp }{\mu }_0^2,\end{array}$$(26)

beyond which the growth rate starts to become smaller than its adiabatic value. The predicted values of δ_c are generally in good agreement with those inferred from Fig. 4. For μ₀ = 0.1 and m = 50 (purple line in the top panel), for example, we obtain δ_c ≈ 90.

For δ ≫ δ_c, the radial wavenumber of the unstable modes becomes significantly smaller than ${\stackrel{\sim }{k}}_c$, and buoyancy can tame the instability effectively. In this regime, we observed that the growth rate scales as Γ ∝ δ⁻² (dashed lines in Fig. 4). This scaling is also reported by Meduri et al. (2025) for the global, radially large-scaled unstable modes of a purely toroidal field. We remark here that, in contrast to the classical m = 1 Tayler modes of a purely toroidal field, since δ_c decreases with decreasing m, the lowest azimuthal wavenumbers of the spectrally selective instability explored here are already influenced by buoyancy even at weak stratification. As a result, the m = 1 mode is the least dominant (blue lines in Fig. 4).

Finally, when increasing the toroidal-to-poloidal field ratio μ₀ from 0.1 to 0.2, δ_c decreases. Gravity therefore starts to tame the instability at smaller values of δ for a fixed value of m, in agreement with the behavior observed in Fig. 4.

At the equator (θ = 90°), the stabilizing buoyancy affects radially large-scale eigenmodes in a manner similar to the polar region. Figure 5 presents the growth rate of the most unstable eigenmodes with azimuthal wavenumbers m = 1 and m = 5 at the equator (dashed lines) and at θ = 10° (solid lines), illustrating that the scaling Γ ∝ δ⁻² for strong stratification is the same in both regions. The increase of the equatorial growth rates with increasing m is also less pronounced than at higher latitudes. In the limit of strong stratification, which is generally relevant to stellar interiors, the instability then grows first close to the axis. Therefore, we focus only on this latter region when discussing the effect of viscosity and resistivity in the next two sections.

Fig. 5. Nondimensional growth rate Γ = σ/ωA0 as a function of the stratification parameter δ for μ0 = 0.1. Dashed lines refer to the equatorial growth rates, while solid lines to the ones at a colatitude θ = 10°.

3.3. Effect of resistivity and viscosity

In this section, we explore the effect of resistivity and viscosity on the instability, focusing on the region near the axis and neglecting stable stratification (δ = 0). We fixed the poloidal-to-toroidal field ratio μ₀ to 0.1, the latitudinal wavenumber l to 10 and the Lundquist number Lu to 10⁴. We first show the results when viscosity is neglected (Pm = 0) and then consider a finite viscosity with Pm = 10⁻³.

For Pm = 0, we find that magnetic diffusion alone cannot fully quench the instability. Instead, it only reduces the growth rate once the characteristic wavenumber of the eigenmodes approaches the resistive wavenumber ${\stackrel{\sim }{k}}_{\eta }=k_{\eta }{r}_{\mathrm{in}}={({\mathrm{\Gamma }}_{\text{ad}}Lu)}^{1/2}$, which we obtained by equating the resistive timescale t_η = (ηk²)⁻¹ with the adiabatic growth timescale 1/σ_ad. Figure 6 (top panel) presents the growth rate Γ as a function of the azimuthal wavenumber m at a colatitude θ = 10°. Modes with m ≲ 10 grow essentially at the adiabatic rate since their characteristic total wavenumbers, which we observed to be compatible with those given by Eq. (23), are significantly smaller than the resistive wavenumber, which is ${\stackrel{\sim }{k}}_{\eta }=100$ for Γ_ad ≈ 1.

Fig. 6. Nondimensional growth rate Γ = σ/ωA0 as a function of the azimuthal wavenumber m close to the axis (θ = 10°) for l = 10, μ0 = 0.1, δ = 0, and Lu = 104. The top and bottom panels present Pm = 0 and Pm = 10−3, respectively. The orange line in the top panel shows the power law Γ ∝ m−2 of the resistive regime. The vertical dashed line in the bottom panel presents the diffusive azimuthal wavenumber (see the main text for details).

Approximating the total non-dimensional wavenumber as ${\stackrel{\sim }{k}}^2\approx {\stackrel{\sim }{k}}_r^2+{\stackrel{\sim }{k}}_{\varphi }^2$ and using the selection condition (23), we estimated the critical azimuthal wavenumber beyond which resistive effects become important by considering

$$\begin{array}{c}\hfill {\stackrel{\sim }{k}}_r^2+{\stackrel{\sim }{k}}_{\varphi }^2={\stackrel{\sim }{k}}_{\eta }^2.\end{array}$$

This equation leads to $m\approx {(1/{\mu }_0^2+1/{\stackrel{\sim }{r}}^2{sin}^2\theta )}^{-1/2}{\stackrel{\sim }{k}}_{\eta }\simeq 9$ for $\stackrel{\sim }{r}\approx 1$ and θ = 10°. Consistently, the top panel of Fig. 6 shows that the growth rate starts to decrease significantly beyond this value, as expected when resistive effects dominate. In this resistive regime, we found that the growth rate scales as Γ ∝ m⁻² (orange line in Fig. 6).

We now include viscosity and consider the case Pm = 10⁻³. Including viscosity introduces an additional timescale into the system, the viscous time t_ν = 1/(νk²). In the low-Pm limit, unstable modes are fully suppressed beyond the diffusive wavenumber ${\stackrel{\sim }{k}}_{\nu \eta }^2\approx {\mathrm{\Gamma }}_{\text{ad}}Lu/\sqrt{\mathit{Pm}}$, obtained by setting the reduced growth rate of the instability Γ ∼ Γ_ad Lu / k² equal to the inverse of the viscous timescale t_ν.

The bottom panel of Fig. 6 shows that modes with m ≲ 10 grow approximately adiabatically, since, for these modes, the wavenumber satisfies $\stackrel{\sim }{k}<{\stackrel{\sim }{k}}_{\eta }$, as discussed in the purely resistive case above. For Γ_ad ≈ 1, the diffusive wavenumber is then ${\stackrel{\sim }{k}}_{\nu \eta }\approx 560$.

The azimuthal wavenumber associated with ${\stackrel{\sim }{k}}_{\nu \eta }$ can be estimated from

$$\begin{array}{c}\hfill {\stackrel{\sim }{k}}^2\approx {\stackrel{\sim }{k}}_r^2+{\stackrel{\sim }{k}}_{\varphi }^2\sim {\stackrel{\sim }{k}}_{\nu \eta }^2.\end{array}$$(27)

Estimating ${\stackrel{\sim }{k}}_r$ from the selection condition (23) and solving for m, we obtained m ≈ 50, shown as a vertical dashed line in the figure, in very good agreement with the numerical results. In the intermediate range of azimuthal wavenumbers corresponding to the modes with ${\stackrel{\sim }{k}}_{\eta }<{\stackrel{\sim }{k}}_r<{\stackrel{\sim }{k}}_{\nu \eta }$, we recovered the scaling Γ ∝ m⁻², as in the Pm = 0 case. The response of this spectrally selective instability to resistive and viscous effects described above closely parallels the behavior reported for the classical Tayler instability of purely toroidal fields (Skoutnev & Beloborodov 2024).

3.4. Combined effects of viscosity, resistivity, and stratification

Under the combined influence of stable stratification, magnetic diffusivity, and viscosity, the behavior of the instability is determined by the relative strength of these effects. In this section, we focus on the case μ₀ = 0.1, Pm = 10⁻³, and ϵ = 0.01. We discuss the instability properties at a colatitude θ = 10° and for a fixed latitudinal wavenumber l = 10.

In order to grow at their adiabatic rate, unstable modes have to satisfy ${\stackrel{\sim }{k}}_r\gg {\stackrel{\sim }{k}}_c$ and $\stackrel{\sim }{k}\ll {\stackrel{\sim }{k}}_{\nu \eta }$, where – as discussed above – ${\stackrel{\sim }{k}}_c$ is the minimum unstable radial wavenumber set by stratification and ${\stackrel{\sim }{k}}_{\nu \eta }$ is the cutoff wavenumber set by viscous and resistive dissipation. At high latitudes, ${\stackrel{\sim }{k}}^2\approx {\stackrel{\sim }{k}}_r^2+{\stackrel{\sim }{k}}_{\varphi }^2$ and we can express the second condition in terms of the radial wavenumber as

$$\begin{array}{c}\hfill {\stackrel{\sim }{k}}_r^2\ll \frac{{\mathrm{\Gamma }}_{\text{ad}}Lu}{\sqrt{\mathit{Pm}}(1+{\mu }^2)}·\end{array}$$

When k_c is larger or equal than the limit radial wavenumber in the condition above, no unstable modes can grow. Stability then occurs when these two limit wavenumbers are equal, that is

$$\begin{array}{c}\hfill \frac{\frac{m}{\stackrel{\sim }{r}sin\theta }\delta }{{\mathrm{\Gamma }}_{\text{ad}}^{1/2}{ϵ}^{1/2}}=\frac{{\mathrm{\Gamma }}_{\text{ad}}Lu}{\sqrt{\mathit{Pm}}(1+{\mu }^2)}·\end{array}$$(28)

We can remove the m dependence from Eq. (28) using the selection condition (23) with ${\stackrel{\sim }{k}}_r={\stackrel{\sim }{k}}_c$ to obtain

$$\begin{array}{c}\hfill \delta =\frac{{\mathrm{\Gamma }}_{\text{ad}}}{\mu }{\left[\frac{ϵLu}{P{m}^{1/2}(1+{\mu }^2)}\right]}^{1/2}.\end{array}$$(29)

Figure 7 shows the instability diagram obtained exploring numerically the parameter space (δ, Lu). Crosses denote stable eigenmodes and circles unstable ones, with colors coding the eigenmode growth rate Γ and sizes representing its azimuthal wavenumber m. The solid blue line corresponds to the stability condition given by Eq. (29) for Γ_ad = 0.8. The numerical results agree well with our theoretical prediction. We indeed found no unstable modes below the predicted stability line (gray shaded region). Weakly unstable eigenmodes close to the stability line show low m_max, while moving deeper into the unstable region selects progressively higher m disturbances. This behavior is consistent with the azimuthal modes’ selection arguments presented above. We remind here that, in contrast to the Tayler instability of purely toroidal fields, the m = 1 mode is the dominant mode only in a few weakly unstable cases close to the stability line.

Fig. 7. Instability diagram showing the Lundquist number Lu as a function of the stratification parameter δ. Numerical results are represented by symbols, where crosses indicate stability while circles mark unstable modes color coded with their instability growth rates Γ. The size of the markers is proportional to the azimuthal wavenumbers m of the most unstable eigenmode, which range between m = 1 and m ∼ 25. The solid blue line presents the stability condition Eq. (29) for Γad = 0.8.

4. Cylindrical coordinates

The linear analysis and results discussed in Sect. 3 were obtained using a local approximation in the latitudinal direction. In this section, we relax this assumption – which may lose accuracy at high latitudes where the instability is most prominent – by performing an analysis in cylindrical coordinates (s, ϕ, z). We compare the results with those obtained for moderate poloidal-to-toroidal field ratios in the previous section and extend the analysis to higher ratios, μ₀ > 1. The linear analysis is global in the cylindrical radial direction s, and neglects viscosity and resistivity. As in Sect. 3, the background magnetic field configuration has a toroidal component increasing linearly with the cylindrical radius,

$$\begin{array}{c}\hfill B_{\varphi }=B_{\varphi 0}s/s_{\text{in}},\end{array}$$

where s_in is the inner cylindrical boundary radius, and a poloidal component characterized by a constant axial field, B_z. The gravitational acceleration is $\mathit{g}=-g_{\text{in}}{\widehat{\mathit{e}}}_z$, where ${\widehat{\mathit{e}}}_z$ is the unit vector in the vertical direction.

Following the approach of Bonanno & Urpin (2011), the system of linearized MHD Equations (3)–(6) with ν = η = 0 can be reduced to two coupled ordinary differential equations for the radial velocity perturbations, u_s′, and the temperature perturbations, T′, which read

$$\begin{array}{cc}& \frac{d}{\mathit{ds}}[\frac{1}{\lambda }({\sigma }^2+{\omega }_A^2)(\frac{d{u}_s\prime }{\mathit{ds}}+\frac{u_s\prime }{s})]-k_z^2({\sigma }^2+{\omega }_A^2)u_s\prime \hfill \\ \hfill & -\frac{2{\omega }_B}{\lambda }\{[(k_z^2-\frac{m^2}{s^2}){\omega }_B-\frac{m{\omega }_{\mathit{Az}}}{s^2}](1-\frac{s}{B_{\varphi }}\frac{d{B}_{\varphi }}{\mathit{ds}})-\frac{2m}{\lambda s^2}{\omega }_A\}{u}_s\prime \hfill \\ \hfill & +\frac{4k_z^2{\omega }_A^2{\omega }_B^2}{\lambda ({\sigma }^2+{\omega }_A^2)}{u}_s\prime =i{k}_z\sigma (1-\frac{m^2}{k_z^2{s}^2\lambda })\frac{\alpha g}{T}\frac{dT\prime }{\mathit{ds}}\hfill \\ \hfill & +\frac{2im\sigma }{\lambda s}(\frac{m}{k_z{s}^2\lambda }+\frac{k_z{\omega }_A{\omega }_B}{{\sigma }^2+{\omega }_A^2})\alpha g\frac{T\prime }{T},\hfill \\ \hfill & \sigma T\prime +\frac{1}{i{k}_z}[\frac{1}{s}\frac{d}{\mathit{ds}}(s{u}_s\prime )-\frac{\mathit{im}}{s}(\frac{m}{i{k}_z^2\lambda s^2}\frac{d}{\mathit{ds}}(s{u}_s\prime )\hfill \\ \hfill & -\frac{2i{\omega }_A{\omega }_B}{\lambda ({\sigma }^2+{\omega }_A^2)}{u}_s\prime -\frac{m\sigma \alpha g}{k_{z}s({\sigma }^2+{\omega }_A^2)}T\prime )]\frac{N^{2}T\prime }{\alpha g_0}\hfill \end{array}$$(30)

$$\begin{array}{cc}& =\kappa [\frac{1}{s}\frac{d}{\mathit{ds}}(s\frac{dT\prime }{\mathit{ds}})-\frac{m^2}{s^2}T\prime -k_z^{2}T\prime ],\hfill \end{array}$$(31)

where

$$\begin{array}{cc}\hfill {\omega }_{\mathit{Az}}& =k_z{B}_z/\sqrt{4\pi \rho },{\omega }_B=B_{\varphi }/s\sqrt{4\pi \rho },\hfill \\ \hfill {\omega }_A& ={\omega }_{\mathit{Az}}+m{\omega }_B,\lambda =1+\frac{m^2}{s^2{k}_z^2}·\hfill \end{array}$$

Here k_z is the vertical wavenumber of the perturbations, which are proportional to exp(σt − ik_zz − imϕ). The Brunt-Väisälä frequency N is defined by

$$\begin{array}{c}\hfill N^2=g_{\mathrm{in}}\frac{\alpha }{T}{\widehat{\mathit{e}}}_z·(\mathbf{\nabla }-{\mathbf{\nabla }}_{\text{ad}})T,\end{array}$$

and the constant background temperature gradient is $(\mathbf{\nabla }-{\mathbf{\nabla }}_{\text{ad}})T\equiv A{\widehat{\mathit{e}}}_z$. We imposed u_s′ = 0 and T′ = 0 at both the cylindrical boundaries s_in and s_out.

To nondimensionalize the equations above, we scaled time in units of the inverse of the toroidal Alfvén frequency ω_A0 = B_ϕ0/(4πρ)^1/2s_in, lengths in units of the inner cylindrical radius s_in, and the magnetic field in units of the background toroidal field amplitude B_ϕ0. The temperature perturbations were nondimensionalized with A s_in. In this scaling scheme, the nondimensional equations read

$$\begin{array}{cc}& \frac{d}{d\stackrel{\sim }{s}}[\frac{1}{\lambda }({\mathrm{\Gamma }}^2+{\stackrel{\sim }{\omega }}_{\mathrm{A}}^2)(\frac{d{\stackrel{\sim }{u}}_s\prime }{d\stackrel{\sim }{s}}+\frac{{\stackrel{\sim }{u}}_s\prime }{\stackrel{\sim }{s}})]-{\stackrel{\sim }{k}}_z^2({\mathrm{\Gamma }}^2+{\stackrel{\sim }{\omega }}_A^2){\stackrel{\sim }{u}}_s\prime \hfill \\ \hfill & -\frac{2f}{\stackrel{\sim }{s}\lambda }\{[({\stackrel{\sim }{k}}_z^2-\frac{m^2}{{\stackrel{\sim }{s}}^2})\frac{f}{\stackrel{\sim }{s}}-\frac{m{\stackrel{\sim }{k}}_z{\mu }_0}{{\stackrel{\sim }{s}}^2}](1-\frac{\stackrel{\sim }{s}}{f}\frac{\mathit{df}}{d\stackrel{\sim }{s}})-\frac{2m}{\lambda {\stackrel{\sim }{s}}^2}{\stackrel{\sim }{\omega }}_A\}{\stackrel{\sim }{u}}_s\prime \hfill \\ \hfill & +\frac{4{\stackrel{\sim }{k}}_z^2{\stackrel{\sim }{\omega }}_A^2{f}^2}{{\stackrel{\sim }{s}}^2\lambda ({\mathrm{\Gamma }}^2+{\stackrel{\sim }{\omega }}_A^2)}{\stackrel{\sim }{u}}_s\prime =i{\stackrel{\sim }{k}}_z\mathrm{\Gamma }(1-\frac{m^2}{{\stackrel{\sim }{k}}_z^2{\stackrel{\sim }{s}}^2\lambda }){\delta }^2\stackrel{\sim }{g}\frac{dT\prime }{d\stackrel{\sim }{s}}\hfill \\ \hfill & +\frac{2im\mathrm{\Gamma }}{\lambda \stackrel{\sim }{s}}(\frac{m}{{\stackrel{\sim }{k}}_z{\stackrel{\sim }{s}}^2\lambda }+\frac{{\stackrel{\sim }{k}}_z{\stackrel{\sim }{\omega }}_{A}f}{{\mathrm{\Gamma }}^2+{\stackrel{\sim }{\omega }}_A^2}){\delta }^2\stackrel{\sim }{g}\stackrel{\sim }{T}\prime ,\hfill \end{array}$$(32)

$$\begin{array}{cc}& \mathrm{\Gamma }\stackrel{\sim }{T}\prime +\frac{1}{i{\stackrel{\sim }{k}}_z}[\frac{1}{\stackrel{\sim }{s}}\frac{d}{d\stackrel{\sim }{s}}(\stackrel{\sim }{s}{\stackrel{\sim }{u}}_s\prime )-\frac{\mathit{im}}{\stackrel{\sim }{s}}(\frac{m}{i{\stackrel{\sim }{k}}_z^2\lambda {\stackrel{\sim }{s}}^2}\frac{d}{d\stackrel{\sim }{s}}(\stackrel{\sim }{s}{\stackrel{\sim }{u}}_s\prime )\hfill \\ \hfill & -\frac{2i{\stackrel{\sim }{\omega }}_{A}f}{\stackrel{\sim }{s}\lambda ({\mathrm{\Gamma }}^2+{\stackrel{\sim }{\omega }}_A^2)}{\stackrel{\sim }{u}}_s\prime -\frac{m\mathrm{\Gamma }\stackrel{\sim }{g}}{{\stackrel{\sim }{k}}_z\stackrel{\sim }{s}({\mathrm{\Gamma }}^2+{\stackrel{\sim }{\omega }}_A^2)}{\delta }^2\stackrel{\sim }{T\prime })]{\stackrel{\sim }{N}}^2\hfill \\ \hfill & =ϵ[\frac{1}{\stackrel{\sim }{s}}\frac{d}{d\stackrel{\sim }{s}}(\stackrel{\sim }{s}\frac{d\stackrel{\sim }{T\prime }}{d\stackrel{\sim }{s}})-\frac{m^2}{{\stackrel{\sim }{s}}^2}\stackrel{\sim }{T\prime }-{\stackrel{\sim }{k}}_z^2\stackrel{\sim }{T\prime }].\hfill \end{array}$$(33)

There are three dimensionless parameters that control the system: the stratification parameter

$$\begin{array}{c}\hfill \delta =\frac{N}{{\omega }_{\mathrm{A}0}},\end{array}$$

the ratio of the thermal diffusion rate to the Alfvén frequency

$$\begin{array}{c}\hfill ϵ=\frac{\kappa }{{\omega }_{A0}{s}_{\text{in}}^2},\end{array}$$

and the vertical to azimuthal field ratio

$$\begin{array}{c}\hfill {\mu }_0=\frac{B_z}{B_{\varphi 0}}·\end{array}$$

The cylindrical radial domain is $1\le \stackrel{\sim }{s}\le 2$, where $\stackrel{\sim }{s}=s/s_{\text{in}}$. To solve the eigenvalue problem above, we discretized the radial derivatives in the above equations using an eight-order finite difference scheme and used the same eigenvalue solver as the one of Sect. 3.

In this section, we cover a broader range of values for the poloidal-to-toroidal field ratio μ₀ than the one explored in the spherical case, spanning values from 0.1 to 100. The following sections discuss first the case with no gravity and then the effect of stable stratification.

4.1. Unstratified case

We first consider the unstratified case (δ = 0) and explore the instability for some fixed values of the azimuthal wavenumber m, varying the field ratio μ₀ in the range 0.1 − 100.

Figure 8 shows the growth rate Γ = σ/ω_A0 as a function of the vertical wavenumber ${\stackrel{\sim }{k}}_z=k_z{r}_{\mathrm{in}}$ for three values of μ₀ = 0.1, 1, and 10 (top to bottom panels) and for the azimuthal wavenumbers m = 1, 10, and 100 (left to right panels). In all cases, the spectrum shows a structure with two clear peaks around the value of the vertical wavenumber obtained from the selection condition (23) (with k_r = k_z) and shown by the vertical dashed lines in the figure. Table 1 reports the vertical wavenumbers associated with the maximum growth rates for the case μ₀ = 0.1. They are very similar to those that were found close to the axis for the analysis in spherical coordinates of the previous section and that are also reported in the table. The spectral width of the unstable band depends on the field ratio μ₀. Larger values of the field ratio yield a narrower range of unstable vertical wavenumbers, whereas smaller values broaden the band and allow a wider range of modes.

Fig. 8. Nondimensional growth rate Γ as a function of the vertical wavenumber ${\stackrel{\sim }{k}}_z$ for μ0 = 0.1 (top panels), μ0 = 1 (middle panels), μ0 = 10 (bottom panels) and different values of the azimuthal wavenumber m. The vertical dashed line in each panel marks the wavenumber ${\stackrel{\sim }{k}}_z$ obtained from Eq. (23) (see the main text for details).

For μ₀ ≤ 1, the maximum growth rates are of the order of the adiabatic rate ω_A0 and tend to increase with the azimuthal wavenumber. For μ₀ > 1, the growth rates become significantly reduced for all values of m explored (see the bottom panels in Fig. 8 for the case at μ₀ = 10). All of the instability properties described above are consistent with the results of the linear analysis presented in Bonanno & Urpin (2011) for our choice of the toroidal field profile.

Figure 9 presents the growth rates of the eigenmodes associated to the most unstable vertical wavenumber as a function of μ₀ for different values of m. The growth rates systematically decrease with μ₀, and rapidly approach the scaling Γ ∝ μ₀⁻¹ (dashed line) when the poloidal component dominates the background field configuration, μ₀ > 1. Eigenmodes with higher azimuthal order m systematically exhibit larger growth rates across the entire range of μ₀ explored. For m ≳ 100, however, the growth rates approach saturation. That is, further increasing the azimuthal wavenumber does not lead to appreciable changes in the growth rate (cf. m = 100 and m = 1000 in Fig. 9). The vanishing growth rate at large μ₀ should not be interpreted as evidence that purely poloidal fields are stable in general, but simply reflects that the selective instability considered here ceases to operate in the limit of a vanishing toroidal component. This does not contradict the classical results of Markey & Tayler (1973) and Flowers & Ruderman (1977) on the instability of purely poloidal fields. Those studies consider different field geometries, in which magnetic field lines possess neutral points either inside or outside the radiative region, and therefore do not correspond to the simple axial configuration explored here. In contrast, our poloidal field is current-free everywhere.

Fig. 9. Growth rate of the most unstable vertical wavenumbers as a function of on the poloidal-to-toroidal field ratio μ0 for different azimuthal wavenumbers m. The dashed line shows the scaling Γ ∝ μ0−1.

4.2. Effect of stable stratification

We now discuss the effect of stable stratification, varying its strength in the range 10 ≤ δ ≤ 10³, similarly to Sect. 3.2, and for μ₀ = 0.1, 1, and 10. Figure 10 presents the results for the case at μ₀ = 1.

Fig. 10. Growth rate of the most unstable vertical wavenumbers as a function of the stratification parameter δ for μ0 = 1 and different values of m as indicated in the legend inset.

Lower m modes, due to the condition (23), are characterized by larger vertical scales and can therefore be significantly affected by stable stratification. Their growth rates follow the approximate scaling Γ ∝ δ⁻². However, for any given stratification strength, there are modes with sufficiently short radial wavelengths or, due to the spectrally selective nature of the instability, sufficiently large m, growing with nearly Alfvénic rates. These results closely parallel the calculations performed in spherical coordinates at high latitudes in Sect. 3.2.

Figure 11 presents the growth rate Γ, for the azimuthal wavenumber m = 100, as a function of δ and μ₀ = 0.1, 1, and 10. For the weak field ratio of μ₀ = 0.1, the growth rate remains essentially Alfvénic over the full range of δ explored. For μ₀ = 1, it decreases significantly, except at low δ, while for μ₀ = 10 the reduction is substantial for all stratification values. As μ₀ increases, we expect, from Eq. (23), the radial wavenumber of the unstable modes to decrease. Modes with small radial wavenumbers are strongly dampened by the stabilizing buoyancy and the critical δ_c below which adiabatic growth is possible becomes smaller (cf. Eq. (26)). For δ ≫ δ_c, the growth rate scales as Γ ∝ δ⁻² (dashed line in Fig. 11).

Fig. 11. Nondimensional growth rate Γ of the most unstable vertical wavenumbers as a function of the stratification parameter δ for m = 100 and different values of μ0. The dashed black line shows the scaling law Γ ∝ δ−2.

5. Summary and conclusions

In this work, we performed a radially global linear analysis of a spectrally selective magnetic instability arising for mixed poloidal–toroidal axisymmetric configurations. The background magnetic field satisfies the magnetohydrostatic equilibrium expected in radiative stellar interiors and consists of an azimuthal component B_ϕ increasing linearly with the cylindrical radius s and a homogeneous axial component B_z.

Our analysis relies on two complementary approaches. First, we adopted spherical coordinates and a poloidal-toroidal decomposition for the perturbations, including stable stratification and all diffusivities. We explored moderate values of the poloidal-to-toroidal ratio μ₀ = B_z/B_ϕ0 < 1, where B_ϕ0 is the azimuthal field strength at the inner boundary on the equator. As the instability is found to be prominent at high latitudes, we then performed a complementary analysis in cylindrical coordinates, valid near the stellar axis and neglecting viscous and resistive effects. This approach, similar to that of Bonanno & Urpin (2011) but including gravity, allowed us to explore the regime μ₀ > 1 where the poloidal field dominates.

Despite sharing several qualitative properties with the classical m = 1 Tayler instability of purely toroidal fields – namely an Alfvénic growth rate for small μ₀, partial resilience to stable stratification, and sensitivity to magnetic diffusion – the instability explored here is fundamentally distinct. It is spectrally selective, that is its unstable wavenumbers k satisfy the selection condition k ⋅ B ≈ 0. For μ₀ ≤ 1 and strongly stable stratification, growth rates of the order of the Alfvén frequency of the background azimuthal field, ω_A0, are reached for azimuthal modes m > 1, which are associated with high radial wavenumbers by the selection condition. When the poloidal field dominates, μ₀ > 1, the instability growth rate is reduced and scales as σ ∝ ω_A/μ₀. We also find that, in contract to the classical Tayler instability of purely toroidal fields, the growth rate σ increases with the azimuthal wavenumber m, but eventually saturates for m ≥ 100. Unstable modes with low radial wavenumbers (or low m) are tamed by the stabilizing buoyancy, following the scaling σ ∝ ω_A³/N², where N is the Brunt-Väisälä frequency. Under the combined effects of stable stratification, thermal diffusion, resistivity, and viscosity, we estimated, using standard WKB assumptions, the marginal stability line below which the instability is suppressed, δ² = Γ_ad² ϵ Lu / Pm^1/2 μ² (1 + μ²). We found excellent agreement between this estimate and our numerical results obtained for μ₀ = 0.1, Pm = 10⁻³, and ϵ = 0.01.

We acknowledge that rotation is dynamically important in radiative stellar interiors and is known to affect the classical Tayler instability of toroidal fields (Pitts & Tayler 1985; Spruit 1999; Skoutnev & Beloborodov 2024). However, its impact on the spectrally selective instability is beyond the scope of this work and will be addressed in a future study including the Coriolis force in the linearized equations.

In conclusion, the global linear analysis presented here provides a comprehensive characterization of a spectrally selective instability of mixed poloidal-toroidal magnetic field configurations, exploring the strongly stably stratified regime relevant to radiative stellar interiors and including all relevant diffusive effects. Our results can provide guidance for future direct numerical simulations and for stellar evolution models aiming to parametrize the effect of magnetic fields on long timescales, as this instability may contribute to the angular momentum transport in radiative interiors.

Acknowledgments

RA acknowledges the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through grant number AR 355/13-1 “Stability and generation of magnetic fields in red giants: towards a unified picture”. DGM acknowledges support from the European Union – Next Generation EU within the framework of the Italian National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.2, through the “Bando Young Researcher 2024” project number SOE2024_0000097 “Elucidating the Transport of Angular momentum by MAgnetic fields in red GIAnts (ETAMAGIA)”, CUP C63C24001200006. AB acknowledges support from the European Union – Next Generation EU RRF M4C2 1.1 n: 2022HY2NSX “CHRONOS: adjusting the clock(s) to unveil the CHRONO-chemo-dynamical Structure of the Galaxy” and from the research grant “Unveiling the magnetic side of the Stars”, funded under the INAF national call for Fundamental Research 2023.

References

All Tables

All Figures

	Fig. 1. Most unstable eigenmodes for the radial velocity perturbations ${\stackrel{\sim }{u}}_r\prime$ in the unstratified case (δ = 0) at a colatitude of θ = 10° for l = 10 and various azimuthal wavenumbers m. The poloidal-to-toroidal field ratio is μ0 = 0.1 and μ0 = 0.2 in the left and right panels, respectively.
In the text

	Fig. 2. Nondimensional growth rate Γ = σ/ωA0 as a function of the latitudinal wavenumber l at the equator for different azimuthal wavenumbers m. The other parameters are δ = 0 and μ0 = 0.1.
In the text

	Fig. 3. Eigenmodes at the equator for δ = 0 and μ0 = 1. The azimuthal wavenumbers m = 10, 20, 40, 100, and 200 are shown, with corresponding latitudinal wavenumber l = 15, 30, 60, 150, and 300 derived from the selection condition (25) at $\stackrel{\sim }{r}=1.5$.
In the text

	Fig. 4. Nondimensional growth rate Γ = σ/ωA0 as a function of the stratification parameter δ for l = 10 at colatitude θ = 10° for two values of μ0 = 0.1 (top panel) and 0.2 (bottom panel), and for various azimuthal wavenumbers m. The dashed black line indicates the scaling Γ ∝ δ−2.
In the text

	Fig. 5. Nondimensional growth rate Γ = σ/ωA0 as a function of the stratification parameter δ for μ0 = 0.1. Dashed lines refer to the equatorial growth rates, while solid lines to the ones at a colatitude θ = 10°.
In the text

Fig. 6. Nondimensional growth rate Γ = σ/ωA0 as a function of the azimuthal wavenumber m close to the axis (θ = 10°) for l = 10, μ0 = 0.1, δ = 0, and Lu = 104. The top and bottom panels present Pm = 0 and Pm = 10−3, respectively. The orange line in the top panel shows the power law Γ ∝ m−2 of the resistive regime. The vertical dashed line in the bottom panel presents the diffusive azimuthal wavenumber (see the main text for details).

In the text

Fig. 7. Instability diagram showing the Lundquist number Lu as a function of the stratification parameter δ. Numerical results are represented by symbols, where crosses indicate stability while circles mark unstable modes color coded with their instability growth rates Γ. The size of the markers is proportional to the azimuthal wavenumbers m of the most unstable eigenmode, which range between m = 1 and m ∼ 25. The solid blue line presents the stability condition Eq. (29) for Γad = 0.8.

In the text

	Fig. 8. Nondimensional growth rate Γ as a function of the vertical wavenumber ${\stackrel{\sim }{k}}_z$ for μ0 = 0.1 (top panels), μ0 = 1 (middle panels), μ0 = 10 (bottom panels) and different values of the azimuthal wavenumber m. The vertical dashed line in each panel marks the wavenumber ${\stackrel{\sim }{k}}_z$ obtained from Eq. (23) (see the main text for details).
In the text

	Fig. 9. Growth rate of the most unstable vertical wavenumbers as a function of on the poloidal-to-toroidal field ratio μ0 for different azimuthal wavenumbers m. The dashed line shows the scaling Γ ∝ μ0−1.
In the text

	Fig. 10. Growth rate of the most unstable vertical wavenumbers as a function of the stratification parameter δ for μ0 = 1 and different values of m as indicated in the legend inset.
In the text

	Fig. 11. Nondimensional growth rate Γ of the most unstable vertical wavenumbers as a function of the stratification parameter δ for m = 100 and different values of μ0. The dashed black line shows the scaling law Γ ∝ δ−2.
In the text