Numerical calculations
We solve for the eigenstates of the linearized anelastic MHD equations30,31 in spherical-polar coordinates (r, θ, ϕ) = (radius, colatitude, longitude). Using R⊙ = 6.96 × 1010 cm for the solar radius, we simulate radii between r0 ≤ r ≤ r1 where r0/R⊙ = 0.89 and r1/R⊙ = 0.99. We place the top of the domain at 99% because several complicated processes quickly increase in importance between this region and the photosphere (for example, partial ionisation, radiative transport and much stronger convection effects). We use the anelastic MHD equations in an adiabatic background to capture the effects of density stratification on the background Alfvén velocities (density varies by roughly a factor of 100 across the NSSL, causing about a factor of 10 change in Alfvén speed) and asymmetries in velocity structures introduced by the density stratification by ∇ ⋅ (ρu). A key aspect of the anelastic approximation is that all entropy perturbations must be small, which is reasonable in the NSSL below 0.99R⊙. We do not use the fully compressible equations, as these linear instability modes do not have acoustic components. Future MRI studies incorporating buoyancy effects (for example, the deep MRI branches at high latitudes) should use a fully compressible (but low Mach number) model32.
Input background parameters
We include density stratification using a low-order polynomial approximation to the Model-S profile33. In units of g cm−3, with h = (r − r0)/(r1 − r0),
$${\rho }_{0}={\alpha }_{0}-{\alpha }_{1}\,h+{\alpha }_{2}\,{h}^{2}-{\alpha }_{3}\,{h}^{3}+{\alpha }_{4}\,{h}^{4},$$
(4)
$${\alpha }_{0}=0.031256,$$
(5)
$${\alpha }_{1}=0.053193,$$
(6)
$${\alpha }_{2}=0.033703,$$
(7)
$${\alpha }_{3}=0.023766,$$
(8)
$${\alpha }_{4}=0.012326,$$
(9)
which fits the Model-S data to better than 1% within the computational domain. The density at h = 1 is ρ0 = 0.000326 compared with 0.031256 at h = 0.
The density profile is close to an adiabatic polytrope with r−2 gravity and 5/3 adiabatic index. An adiabatic background implies that buoyancy perturbations diffuse independently of the MHD and decouple from the system.
We use a low-degree polynomial fit to the observed NSSL differential rotation profile. For μ = cos(θ),
$${{\bf{u}}}_{0}=\varOmega (r,\theta )\,r\sin (\theta )\,{{\bf{e}}}_{\phi },$$
(10)
$$\varOmega (r,\theta )={\varOmega }_{0}\,R(h)\,\Theta (\mu ),$$
(11)
where Ω0 = 466 nHz ≈ 2.92 × 10−6 s−1 and
$$R(h)=1+0.02\,h-0.01\,{h}^{2}-0.03\,{h}^{3},$$
(12)
$$\Theta (\mu )=1-0.145\,{\mu }^{2}-0.148\,{\mu }^{4}.$$
(13)
We use the angular fit from ref. 34. The radial approximation results from fitting the equatorial profile from ref. 29 shown in Fig. 1a. Below 60° latitude, the low-degree approximation agrees with the full empirical profile to within 1.25%. The high-latitude differential rotation profile is less constrained because of observational uncertainties.
We define the background magnetic field in terms of a vector potential,
$${{\bf{B}}}_{0}={\boldsymbol{\nabla }}\times {{\bf{A}}}_{0},$$
(14)
$${{\bf{A}}}_{0}=\frac{{\mathcal{B}}(r)}{2}\,r\sin (\theta )\,{{\bf{e}}}_{\phi },$$
(15)
where
$${\mathcal{B}}(r)={B}_{0}\,\left({(r/{r}_{1})}^{-3}-{(r/{r}_{1})}^{2}\right),$$
(16)
and B0 = 90 G. The r−3 term represents a global dipole. The r2 term represents a field with a similar structure but containing electric current,
$${{\bf{J}}}_{0}=\frac{{\boldsymbol{\nabla }}\times {{\bf{B}}}_{0}}{4{\rm{\pi }}}=\frac{5{B}_{0}}{4{\rm{\pi }}\,{r}_{1}^{2}}\,r\sin (\theta )\,{{\bf{e}}}_{\phi }.$$
(17)
The background field is in MHD force balance,
$${{\bf{J}}}_{0}\,\times \,{{\bf{B}}}_{0}={\boldsymbol{\nabla }}({{\bf{A}}}_{0}\cdot {{\bf{J}}}_{0}\,).$$
(18)
The MHD force balance generates magnetic pressure, which inevitably produces entropy, s′, and enthalpy, h′, perturbations using
$$\frac{{\boldsymbol{\nabla }}({{\bf{A}}}_{0}\cdot {{\bf{J}}}_{0})}{{\rho }_{0}}+{T}_{0}{\boldsymbol{\nabla }}{s}^{{\prime} }={\boldsymbol{\nabla }}{h}^{{\prime} },$$
(19)
where
$${s}^{{\prime} }=\frac{1}{{\varGamma }_{3}-1}\frac{{{\bf{A}}}_{0}\cdot {{\bf{J}}}_{0}}{{T}_{0}\,{\rho }_{0}},\quad {h}^{{\prime} }=\frac{{\varGamma }_{3}}{{\varGamma }_{3}-1}\frac{{{\bf{A}}}_{0}\cdot {{\bf{J}}}_{0}}{{\rho }_{0}},$$
(20)
and Γ3 is the third adiabatic index. However, the MRI is a weak-field instability, implying magnetic buoyancy and baroclinicity effects are subdominant. For the work presented here, we neglect the contributions of magnetism to entropy (magnetic buoyancy) and consider adiabatic motions. We expect this to be valid for MRI in the NSSL, but studies of MRI in the deep convection zone at high latitudes would need to incorporate these neglected effects.
We choose our particular magnetic field configuration rather than a pure dipole because the radial component \({{\bf{e}}}_{r}\cdot {{\bf{B}}}_{0}={\mathcal{B}}(r)\cos (\theta )\) vanishes at r = r1. The poloidal field strength in the photosphere is about 1 G, but measurements suggest sub-surface field strengths of about 50–150 G (ref. 9). The near-surface field should exhibit a strong horizontal (as opposed to radial) character. Magnetic pumping35 by surface granulation within the outer 1% of the solar envelope could account for filtering the outward radial field, with sunspot cores being prominent exceptions.
We also test pure dipoles and fields with an approximately 5% dipole contribution, yielding similar results. Furthermore, we test that the poloidal field is stable to current-driven instabilities. Our chosen confined field also has the advantage that eθ ⋅ B0 is constant to within 8% over r0 < r < r1. However, a pure dipole varies by about 37% across the domain. The RMS field amplitude is ∣B∣RMS ≈ 2B0 = 180 G, about 25% larger than the maximum-reported inferred dipole equivalent9. However, projecting our field onto a dipole template gives an approximately 70 G equivalent at the r = r1 equator. Overall, the sub-surface field is the least constrained input to our calculations, the details of which change over several cycles.
Model equations
Respectively, the linearized anelastic momentum, mass-continuity and magnetic induction equations are
$${\rho }_{0}({\partial }_{t}{\bf{u}}+{{\boldsymbol{\omega }}}_{0}\times {\bf{u}}+{\boldsymbol{\omega }}\times {{\bf{u}}}_{0}+{\boldsymbol{\nabla }}\varpi )=\nu {\boldsymbol{\nabla }}\cdot ({\rho }_{0}{\boldsymbol{\sigma }})+{\bf{j}}\times {{\bf{B}}}_{0}+{{\bf{J}}}_{0}\times {\bf{b}},$$
(21)
$${\boldsymbol{\nabla }}\cdot \left({\rho }_{0}{\bf{u}}\right)=0,$$
(22)
$${\partial }_{t}{\bf{b}}-\eta {\nabla }^{2}{\bf{b}}={\boldsymbol{\nabla }}\times \left({{\bf{u}}}_{0}\times {\bf{b}}+{\bf{u}}\times {{\bf{B}}}_{0}\right),$$
(23)
where the traceless strain rate
$${\boldsymbol{\sigma }}\,=\,{\boldsymbol{\nabla }}{\bf{u}}+{({\boldsymbol{\nabla }}{\bf{u}})}^{{\rm{\top }}}-\frac{2}{3}{\boldsymbol{\nabla }}\cdot {\bf{u}}\,{\bf{I}}.$$
(24)
To find eigenstates, ∂t → γ + iω, where γ is the real-valued growth rate, and ω is a real-valued oscillation frequency. The induction equation (23) automatically produces MRI solutions satisfying ∇ ⋅ b = 0.
Given the velocity perturbation, u, the vorticity ω = ∇ × u. Given the magnetic field (Gauss in cgs units), the current density perturbations j = ∇ × b/4π. At linear order, the Bernoulli function \(\varpi ={{\bf{u}}}_{0}\cdot {\bf{u}}+{h}^{{\prime} }\), where h′ represents enthalpy perturbations26.
The velocity perturbations are impenetrable (ur = 0) and stress-free (σrθ = σrϕ = 0) at both boundaries. For the magnetic field, we enforce perfect conducting conditions at the inner boundary (br = ∂rbθ = ∂rbϕ = 0). At the outer boundary, we test three different choices in common usage, as different magnetic boundary conditions have different implications for magnetic helicity fluxes through the domain, and these can affect global dynamo outcomes36. Two choices with zero helicity flux are perfectly conducting and vacuum conditions, and we find only modest differences in the results. We also test a vertical field or open boundary (that is, ∂rbr = bθ = bϕ = 0), which, although non-physical, explicitly allows a helicity flux. These open systems also had essentially the same results as the other two for growth rates and properties of eigenfunctions. We conduct most of our experiments using perfectly conducting boundary conditions, which we prefer on the same physical grounds as the background field.
We set constant and kinematic viscous and magnetic diffusivity parameters ν = η = 10−6 in units where Ω0 = R⊙ = 1. The magnetic Prandtl number ν/η = Pm = 1 assumes equal transport of vectors by the turbulent diffusivities. A more detailed analysis of the shear Reynolds numbers yields \({\rm{Re}}={\rm{Rm}}={U}_{0}\,{L}_{0}/\nu \approx \mathrm{1,500}\), where U0 ≈ 5,000 cm s−1 is the maximum shear velocity jump across the NSSL and L0 ≈ 0.06R⊙ is the distance between minimum and maximum shear velocity (see section ‘NSSL energetics and turbulence parameterization’ below).
We compute the following scalar-potential decompositions a posteriori,
$${\bf{u}}={u}_{\phi }\,{{\bf{e}}}_{\phi }+\frac{1}{{\rho }_{0}}{\boldsymbol{\nabla }}\times ({\rho }_{0}\,\psi \,{{\bf{e}}}_{\phi }),$$
(25)
$${\bf{b}}={b}_{\phi }\,{{\bf{e}}}_{\phi }+{\boldsymbol{\nabla }}\times ({a}_{\phi }\,{{\bf{e}}}_{\phi }),$$
(26)
where both the magnetic scalar potential, aϕ, and the streamfunction, ψ, vanish at both boundaries.
We, furthermore, compute the current helicity correlation relative to global RMS values,
$${\mathcal{H}}=\frac{{\bf{b}}\cdot {\bf{j}}}{| {\bf{b}}{| }_{{\rm{RMS}}}\,| \,{\bf{j}}{| }_{{\rm{RMS}}}}.$$
(27)
There is no initial helicity in the background poloidal magnetic field,
$${{\bf{B}}}_{0}={\boldsymbol{\nabla }}\times ({A}_{0}(r,\theta ){{\bf{e}}}_{\phi })\Rightarrow {{\bf{B}}}_{0}\cdot ({\boldsymbol{\nabla }}\times {{\bf{B}}}_{0})=0.$$
Linear dynamical perturbations, b(r, θ), will locally align with the background field and current. However, because the eigenmodes are wave-like, these contributions vanish exactly when averaged over hemispheres.
$$\langle {\bf{b}}\cdot ({\boldsymbol{\nabla }}\times {{\bf{B}}}_{0})\rangle =\langle {{\bf{B}}}_{0}\cdot ({\boldsymbol{\nabla }}\times {\bf{b}})\rangle =0.$$
The only possible hemispheric contributions arise when considering quadratic mode interactions,
$$\langle {\bf{b}}\cdot ({\boldsymbol{\nabla }}\times {\bf{b}})\rangle \ne 0.$$
This order is the first for which we could expect a non-trivial signal.
Finally, we also solve the system using several different mathematically equivalent equation formulations (for example, using a magnetic vector potential b = ∇ × a, or dividing the momentum equations by ρ0). In all cases, we find excellent agreement in the converged solutions. We prefer this formulation because of satisfactory numerical conditioning as parameters become more extreme.
Computational considerations
The Dedalus code25 uses general tensor calculus in spherical-polar coordinates using spin-weighted spherical harmonics in (θ, ϕ) (refs. 37,38). For the finite radial shell, the code uses a weighted generalized Chebyshev series with sparse representations for differentiation, radial geometric factors and non-constant coefficients (for example, ρ0(r)). As the background magnetic field and the differential rotation are axisymmetric and they contain only a few low-order separable terms in latitude and radius, these two-dimensional non-constant coefficients have a low-order representation in a joint expansion of spin-vector harmonics and Chebyshev polynomials. The result is a two-dimensional generalized non-Hermitian eigenvalue problem Ax = λBx, where x represents the full system spectral-space state vector. The matrices, A and B, are spectral-coefficient representations of the relevant linear differential and multiplication operators. Cases 1 and 2 use 384 latitudinal and 64 radial modes (equivalently spatial points). The matrices A and B remain sparse, with respective fill factors of about 8 × 10−4 and 4 × 10−5.
The eigenvalues and eigenmodes presented here are converged to better than 1% relative absolute error (comparing 256 and 384 latitudinal modes). We also use two simple heuristics for rejecting poorly converged solutions. First, because λ0 is complex valued, the resulting iterated solutions do not automatically respect Hermitian-conjugate symmetry, which we often find violated for spurious solutions. Second, the overall physical system is reflection symmetric about the equator, implying the solutions fall into symmetric and anti-symmetric classes. Preserving the desired parity is a useful diagnostic tool for rejecting solutions with mixtures of the two parities, which we check individually for each field quantity. The precise parameters and detailed implementation scripts are available at GitHub (https://github.com/geoffvasil/nssl_mri).
Analytic and semi-analytic estimates
Local equatorial calculation
Our preliminary estimates of the maximum poloidal field strength involve solving a simplified equatorial model of the full perturbation equations, setting the diffusion parameters ν, η → 0. Using a Lagrangian displacement vector, ξ, in Eulerian coordinates
$${\bf{u}}={{\rm{\partial }}}_{t}{\boldsymbol{\xi }}+{{\bf{u}}}_{0}\cdot {\boldsymbol{\nabla }}{\boldsymbol{\xi }}-{\boldsymbol{\xi }}\cdot {\boldsymbol{\nabla }}{{\bf{u}}}_{0},$$
(28)
$${\bf{b}}={\boldsymbol{\nabla }}\times ({\boldsymbol{\xi }}\times {{\bf{B}}}_{0}).$$
(29)
In local cylindrical coordinates near the equator (r, ϕ, z), we assume all perturbations are axis-symmetric and depend harmonically \(\propto {{\rm{e}}}^{{\rm{i}}({k}_{z}z-\omega t)}\). The cylindrical assumption simplifies the analytical calculations while allowing a transference of relevant quantities from the more comprehensive spherical model. That is, we assume a purely poloidal background field with the same radial form as the full spherical computations, B0 = Bz(r)ez. We use the same radial density and angular rotation profiles, ignoring latitudinal dependence. The radial displacement, ξr, determines all other dynamical quantities,
$${\xi }_{\phi }=-\frac{2{\rm{i}}\omega \,\varOmega }{{\omega }^{2}-{k}_{z}^{2}{v}_{{\rm{A}}}^{2}}{\xi }_{r},$$
(30)
$${\xi }_{z}=\frac{{\rm{i}}}{{k}_{z}\,r\,{\rho }_{0}}\frac{{\rm{d}}(r{\rho }_{0}{\xi }_{r})}{{\rm{d}}r}$$
(31)
$$\varpi ={v}_{{\rm{A}}}^{2}\frac{{B}_{z}^{{\prime} }}{{B}_{z}}{\xi }_{r}+\frac{{\omega }^{2}}{{k}_{z}^{2}\,r\,{\rho }_{0}}\frac{{\rm{d}}(r{\rho }_{0}{\xi }_{r})}{{\rm{d}}r},$$
(32)
where \({v}_{{\rm{A}}}(r)={B}_{z}(r)/\sqrt{4{\rm{\pi }}{\rho }_{0}(r)}\). The radial momentum equation gives a second-order two-point boundary-value problem for ξr(r). The resulting real-valued differential equation depends on ω2; the instability transitions directly from oscillations to exponential growth when ω = 0. We eliminate terms containing \({\xi }_{r}^{{\prime} }(r)\) with the Liouville transformation \(\varPsi (r)=\sqrt{r}{B}_{z}(r){\xi }_{r}(r)\). The system for the critical magnetic field reduces to a Schrödinger-type equation,
$$-{\varPsi }^{{\prime\prime} }(r)+{k}_{z}^{2}\,\varPsi (r)+V(r)\,\varPsi (r)\,=\,0,$$
(33)
with boundary conditions
$$\varPsi (r={r}_{0})\,=\,\varPsi (r={r}_{1})\,=\,0$$
(34)
and potential,
$$V=\frac{r}{{v}_{{\rm{A}}}^{2}}\frac{{\rm{d}}{\varOmega }^{2}}{{\rm{d}}r}+\frac{r{\rho }_{0}}{{B}_{z}}\frac{{\rm{d}}}{{\rm{d}}r}\,\,\left(\,\frac{1}{r{\rho }_{0}}\frac{{\rm{d}}{B}_{z}}{{\rm{d}}r}\,\right)+\frac{3}{4{r}^{2}}.$$
(35)
Upper bound
The maximum background field strength occurs in the limit kz → 0. With fixed functional forms for Ω(r), ρ0(r), we suppose
$$\begin{array}{c}{B}_{z}(r)\,=\,{B}_{1}\frac{1+4{(r/{r}_{1})}^{5}}{5{(r/{r}_{1})}^{3}},\end{array}$$
(36)
with B1 = Bz(r1) setting and overall amplitude and \(1/{B}_{1}^{2}\) serving as a generalized eigenvalue parameter. We solve the resulting system with Dedalus using both Chebyshev and Legendre series for 64, 128 and 256 spectral modes, all yielding the same result, B1 = 1,070 G. The results are also insensitive to detailed changes in the functional form of the background profile.
Growth rate
We use a simplified formula for the MRI exponential growth, proportional to eγt, in a regime not extremely far above onset22. That is,
$${\gamma }^{2}\,\approx \,\frac{{\alpha }^{2}{\omega }_{{\rm{A}}}^{2}\,(2\varOmega S-{\omega }_{{\rm{A}}}^{2}\,(1+{\alpha }^{2}))}{{\omega }_{{\rm{A}}}^{2}+4{\varOmega }^{2}},$$
(37)
where α = 2H/L ≈ 0.2–0.3 is the mode aspect ratio with latitudinal wavelength, L ≈ 20°–30°R⊙, and NSSL depth H ≈ 0.05R⊙. The main text defines all other parameters. In the NSSL, S ≈ Ω. Therefore, γ/Ω ≈ 0.1, when α ≈ 0.3 and ωA/Ω ≈ 1; and γ/Ω ≈ 0.01, when α ≈ 0.2 and ωA/Ω ≈ 0.1.
Saturation amplitude
We use non-dissipative quasi-linear theory22 to estimate the amplitude of the overall saturation. In a finite-thickness domain, the MRI saturates by transporting mean magnetic flux and angular momentum radially. Both quantities are (approximately) globally conserved; however, the instability shifts the magnetic flux inward and angular momentum outward, so the potential from equation (35) is positive everywhere in the domain.
Given the cylindrical radius, r, the local angular momentum and magnetic flux density
$$L={\rho }_{0}r{u}_{\phi },\,M={\rho }_{0}r{a}_{\phi }.$$
(38)
The respective local flux transport
$${{\rm{\partial }}}_{t}L+{\boldsymbol{\nabla }}\cdot (L{\bf{u}})={\boldsymbol{\nabla }}\cdot (r\,{b}_{\phi }{\bf{b}}),$$
(39)
$${{\rm{\partial }}}_{t}M+{\boldsymbol{\nabla }}\cdot (M{\bf{u}})=0.$$
(40)
For quadratic-order feedback,
$${\partial }_{t}({\rho }_{0}{r}^{2}\delta {u}_{\phi })={\partial }_{r}({r}^{2}({b}_{\phi }{b}_{r}-{\rho }_{0}{u}_{\phi }{u}_{r}))+{\partial }_{z}({r}^{2}({b}_{\phi }{b}_{z}-{\rho }_{0}{u}_{\phi }{u}_{z})),$$
(41)
$${{\rm{\partial }}}_{t}({\rho }_{0}{r}^{2}\delta {a}_{\phi })=-{{\rm{\partial }}}_{r}({r}^{2}{\rho }_{0}{a}_{\phi }{u}_{r})-{{\rm{\partial }}}_{z}({r}^{2}{\rho }_{0}{a}_{\phi }{u}_{z}).$$
(42)
For linear meridional perturbations,
$${u}_{r}=-{{\rm{\partial }}}_{z}\psi ,\,{u}_{z}=\frac{{{\rm{\partial }}}_{r}(r{\rho }_{0}\psi )}{r{\rho }_{0}},$$
(43)
$${b}_{r}=-{{\rm{\partial }}}_{z}{a}_{\phi },\,{b}_{z}=\frac{{{\rm{\partial }}}_{r}(r{a}_{\phi })}{r}.$$
(44)
For the angular components,
$${\partial }_{t}{u}_{\phi }={\partial }_{z}\,\left(\,(2\varOmega -S)\,\psi +\frac{{B}_{z}}{4{\rm{\pi }}{\rho }_{0}}{b}_{\phi }\right),$$
(45)
$${\partial }_{t}{a}_{\phi }={\partial }_{z}({B}_{z}\psi ),$$
(46)
$${{\rm{\partial }}}_{t}{b}_{\phi }={{\rm{\partial }}}_{z}({B}_{z}{u}_{\phi }+S\,{a}_{\phi }).$$
(47)
Using the linear balances, we time integrate to obtain the latitudinal-mean rotational and magnetic feedback,
$$\delta \varOmega =\frac{1}{{r}^{3}{\rho }_{0}}{\partial }_{r}\left({r}^{2}{\rho }_{0}\,{\mathcal{L}}\right),$$
(48)
$$\delta A=\frac{1}{{r}^{2}{\rho }_{0}}{\partial }_{r}\left({r}^{2}{\rho }_{0}\,\Phi \right).$$
(49)
where angle brackets represent z averages and
$${\mathcal{L}}=\frac{2{B}_{z}\langle {a}_{\phi }{u}_{\phi }\rangle -(2\varOmega -S)\,\langle {a}_{\phi }^{2}\rangle }{2{B}_{z}^{2}},$$
(50)
$$\Phi =\frac{\langle {a}_{\phi }^{2}\rangle }{2{B}_{z}}.$$
(51)
The dynamic shear and magnetic corrections,
$$\delta S=-r{{\rm{\partial }}}_{r}\delta \varOmega ,\,\delta {B}_{z}=\frac{1}{r}{{\rm{\partial }}}_{r}(r\delta A).$$
(52)
We derive an overall amplitude estimate by considering the functional
$${\mathcal{F}}=\int (V| \varPsi {| }^{2}+| {\boldsymbol{\nabla }}\varPsi {| }^{2}){\rm{d}}r,$$
(53)
which results from integrating equation (33) with respect to Ψ*(r). The saturation condition is
$$\delta {\mathcal{F}}=-{\mathcal{F}}.$$
(54)
The left-hand side includes all linear-order perturbations in the potential, δV, and wavefunction, δΨ, where
$$\begin{array}{l}\delta V=\frac{2r}{{v}_{{\rm{A}}}^{2}}\frac{{\rm{d}}(\varOmega \delta \varOmega )}{{\rm{d}}r}-2\frac{\delta {B}_{z}}{{B}_{z}}\frac{r}{{v}_{{\rm{A}}}^{2}}\frac{{\rm{d}}{\varOmega }^{2}}{{\rm{d}}r}\\ \,+\,\frac{r{\rho }_{0}}{{B}_{z}}\frac{{\rm{d}}}{{\rm{d}}r}\,\,\left(\,\frac{1}{r{\rho }_{0}}\frac{{\rm{d}}\delta {B}_{z}}{{\rm{d}}r}\,\right)-\frac{\delta {B}_{z}}{{B}_{z}}\frac{r{\rho }_{0}}{{B}_{z}}\frac{{\rm{d}}}{{\rm{d}}r}\,\,\left(\,\frac{1}{r{\rho }_{0}}\frac{{\rm{d}}{B}_{z}}{{\rm{d}}r}\,\right),\end{array}$$
(55)
$$\delta \varPsi =\frac{\delta {B}_{z}}{{B}_{z}}\,\varPsi .$$
(56)
All reference and perturbation quantities derive from the full sphere numerical eigenmode calculation. We translate to cylindrical coordinates by approximating z averages with latitudinal θ averages. The spherical eigenmodes localize near the equator, and the NSSL thickness is only about 5% of the solar radius, justifying the cylindrical approximation in the amplitude estimate.
Empirically, the first δV term dominates the overall feedback calculation, owing to the shear corrections \(\propto \,{\rm{d}}\delta \varOmega /{\rm{d}}r \sim 1/{H}_{r}^{2}\). Isolating the shear effect produces the simple phenomenological formula in equation (3).
NSSL energetics and turbulence parameterization
We estimate that the order-of-magnitude energetics of the NSSL are consistent with the amplitudes of torsional oscillations. The torsional oscillations comprise |Ω′| ≈ 1 nHz rotational perturbation, relative to the Ω0 ≈ 466 nHz equatorial frame rotation rate. However, the NSSL contains ΔΩ ≈ 11 nHz mean rotational shear estimated from the functional form in equations (10)–(13). In terms of velocities, the shear in the NSSL has a peak contrast of roughly U0 ≈ 5 × 103 cm s−1 across a length scale L0 ≈ 0.06R⊙. The relative amplitudes of the torsional oscillations to the NSSL background, |Ω′|/ΔΩ, are thus about 10%. Meanwhile, the radial and latitudinal global differential rotations have amplitudes of the order of about 100 nHz. The relative energies are approximately the squares of these, implying that the ΔKE of the torsional oscillations is about 0.01% to the differential rotation and about 1% to the NSSL. These rough estimates show that the NSSL and the differential rotation can provide ample energy reservoirs for driving an MRI dynamo, and the amplitude of the torsional oscillations is consistent with nonlinear responses seen in classical convection-zone dynamos17.
Vigorous hydrodynamic convective turbulence probably establishes the differential rotation of the NSSL. The large reservoir of shear in the solar interior plays the analogue part of gravity and Keplerian shear in accretion disks. The details of solar convection are neither well understood nor well constrained by observations. There are indications, however, that the maintenance of the NSSL is separate from the solar cycle because neither the global differential rotation nor the NSSL shows substantial changes during the solar cycle other than in the torsional oscillations.
Strong dynamical turbulence in the outer layers of the Sun is an uncertainty of our MRI dynamo framework, but scale separation gives hope for progress. From our linear instability calculations, the solar MRI operates relatively close to the onset and happens predominantly on large scales. If the fast turbulence of the outer layers of the Sun acts mainly as an enhanced dissipation, then the solar MRI should survive relatively unaffected. Treating scale-separated dynamics in this fashion has good precedent: large-scale baroclinic instability in the atmosphere of Earth gracefully ploughs through the vigorous moist tropospheric convection (thunderstorms). Scale-separated dynamics are particularly relevant because the MRI represents a type of essentially nonlinear dynamo, which cyclically reconfigures an existing magnetic field using kinetic energy as a catalyst. From previous work, it is clear that the deep solar convection zone can produce global-scale fields, but these fields generally have properties very different from the observed fields17. Essentially nonlinear MHD dynamos have analogues in pipe turbulence, and, similar to those systems, the self-sustaining process leads to an attractor in which the dynamo settles into a cyclic state independent of its beginnings.
A full nonlinear treatment of turbulence in the NSSL-MRI setting awaits future work. Here we adopt a simplified turbulence model using enhanced dissipation. To model the effects of turbulence, we assume that the viscous and magnetic diffusivities are enhanced such that the turbulent magnetic Prandtl number Pm = 1 (with no principle of turbulence suggesting otherwise). The momentum and magnetic Reynolds numbers are Re = Rm ≈ 1.5 × 103. These values are vastly more dissipative than the microphysical properties of solar plasma (that is, Re ~ 1012), and the microphysical Pm ≪ 1, implying that Rm ≪ Re. The studies conducted here find relative independence in the MRI on the choices of Re within a modest range. By contrast, other instabilities (for example, convection) depend strongly on Re. We compute sample simulations down to Re ≈ 50 with qualitatively similar results, although they match the observed patterns less well and require somewhat stronger background fields. Our adopted value of Re ≈ 1,500 strikes a good balance for an extremely under-constrained process. Our turbulent parameterizations also produce falsifiable predictions: our proposed MRI dynamo mechanism would face severe challenges if future helioseismic studies of the Sun suggest that the turbulent dissipation is much larger than expected (for example, if the effective Re ≪ 1). However, it is difficult to imagine how any nonlinear dynamics would happen in this scenario.
