The solar dynamo begins near the surface

The magnetic dynamo cycle of the Sun features a distinct pattern: a propagating region of sunspot emergence appears around 30° latitude and vanishes near the equator every 11 years (ref. 1). Moreover, longitudinal flows called torsional oscillations closely shadow sunspot migration, undoubtedly sharing a common cause2. Contrary to theories suggesting deep origins of these phenomena, helioseismology pinpoints low-latitude torsional oscillations to the outer 5–10% of the Sun, the near-surface shear layer3,4. Within this zone, inwardly increasing differential rotation coupled with a poloidal magnetic field strongly implicates the magneto-rotational instability5,6, prominent in accretion-disk theory and observed in laboratory experiments7. Together, these two facts prompt the general question: whether the solar dynamo is possibly a near-surface instability. Here we report strong affirmative evidence in stark contrast to traditional models8 focusing on the deeper tachocline. Simple analytic estimates show that the near-surface magneto-rotational instability better explains the spatiotemporal scales of the torsional oscillations and inferred subsurface magnetic field amplitudes9. State-of-the-art numerical simulations corroborate these estimates and reproduce hemispherical magnetic current helicity laws10. The dynamo resulting from a well-understood near-surface phenomenon improves prospects for accurate predictions of full magnetic cycles and space weather, affecting the electromagnetic infrastructure of Earth.

surface-shear layer (NSSL) within the Sun's outer 5-10% containing strong inwardly increasing angular velocity fostering the Magneto-Rotational Instability (MRI).Despite progress, prevailing theories have distinct limitations.Interface dynamos (proposed within the tachocline 8 ) preferentially generate high-latitude fields 12 , and would produce severe shear disruptions 13 which are not observed 14 .Mean-field dynamos offer qualitative insights but suffer from the absence of first principles 15 , and are contradicted by observed meridional circulations 16 .
Global convection-zone models often misalign with critical solar observations, require conditions diverging from solar reality [17][18][19] , and fail to provide a theoretical dynamical understanding.
Borrowing from well-established ideas in accretion-disk physics 5,6 , we propose an alternative hypothesis that produces clear predictions and quantitatively matches key observations.For electrically conducting plasma like the Sun, the local axisymmetric linear instability criterion for the MRI is 5,6 : where the local Alfvén frequency and shear are The system control parameters are the background poloidal magnetic field strength (B 0 in cgs units), the atmospheric density (ρ 0 ), the smallest non-trivial radial wavenumber that will fit in the domain (k r ≈ π/H r , where H r is a relevant layer depth or density-scale height), bulk rotation rate (Ω), and the differential rotation, or shear (S > 0 in the NSSL).An adiabatic density stratification holds to a good approximation within the solar convection zone, eliminating buoyancy modifications to the stability condition (1).
The MRI is essential for generating turbulent angular momentum transport in magnetized astrophysical disks 6 .Previous work 21 postulated the NSSL as the possible seat of the global dynamo without invoking the MRI.Another kinematic-dynamo study 22 dismissed the relevance of NSSL for kinematic dynamos but did not allow for full magnetohydrodynamic (MHD) instabilities (such as the MRI).Modern breakthroughs in our understanding of large-scale MRI physics 23,24 have not been applied in a solar context, and local MRI studies of the Sun 25 have only considered small scales.To our knowledge, no work has yet considered the large-scale MRI dynamics relevant to the observed features of the global dynamo.We, therefore, describe here a potential MRI-driven solar dynamo cycle.
The start of the solar cycle is the period surrounding the sunspot minimum when there is no significant toroidal field above the equator and a maximal poloidal field below the photosphere.
This configuration is unstable to the axisymmetric MRI, which generates a dynamically active toroidal field in the outer convection zone.The observed torsional oscillations are the longitudinal flow perturbations arising from the MRI.The relative energetics are consistent with nonlinear dynamo estimates (see Methods).As the cycle progresses, the toroidal field can undergo several possible MHD instabilities contributing to poloidal-field regeneration, e.g., the helical MRI, nonaxis-symmetric MRI, the "clamshell" instability, and several more, including a surface Babcock-Leighton process.We hypothesise that the axisymmetric subsurface field and torsional oscillations constitute a nonlinear MRI travelling wave.The instability saturates via radial transport of (globally conserved) mean magnetic flux (B 0 ) and angular momentum (Ω, S), which neutralise the instability criterion eq. ( 1) (see Methods).
Empirical timescales of the torsional oscillations imply an approximate growth rate, γ, for the MRI and, hence, a relevant poloidal field strength.Background shear modification dominates the MRI saturation mechanism (see Methods), roughly where Ω ′ represents the dynamic changes in differential rotation.For S ≈ Ω ≈ ω A , |Ω ′ | ≈ 7 nHz, roughly consistent with the observed torsional oscillation amplitude; see fig.1(d).
We compute a suite of growing global perturbations using Dedalus 26 to model the initial phase of the solar cycle with quasi-realistic solar input parameters (see Methods).Fig. 2 shows representative solutions.
We find two distinct cases: (i) a "fast branch" with direct growth rates, γ, comparable to a priori estimates; and (ii) a "slow branch" with longer but relevant growth times and oscillation periods.
The eigenmodes are confined to the NSSL, reaching from the surface to r/R ⊙ ≈ 0.90-0.95,where the background shear becomes MRI stable.
For case (i), γ/Ω 0 ≈ 6•10 −2 (given Ω 0 = 466 nHz) with corresponding e-folding time, t e ≈ 60 days and no discernible oscillation frequency.The pattern comprises roughly one wave period between the equator and ≈ 20 • latitude, similar to the rotation perturbations seen in the torsional oscillations.
In addition to cases (i) and (ii), we find 34 additional purely growing "fast-branch" modes, two additional oscillatory modes and one intermediate exceptional mode (see Extended Data figs.1-3).
Using the full numerical MHD eigenstates, we compute a systematic estimate for the saturation amplitude using quasi-linear theory (see Methods): |Ω ′ | ≈ 6 nHz for case (i) and |Ω ′ | ≈ 3 nHz for the slow-branch case (ii); both comparable to the observed torsional oscillation amplitude and the simple analytical estimates from eq. ( 3).The true saturated state would comprise an interacting superposition of the full spectrum of modes.
Notably, the slow-branch current helicity, H ∝ b • ∇×b, follows the hemispherical sign rule 10 , with H < 0 in the north and H > 0 in the south.The slow-branch modes appear rotationally constrained, consistent with their low Rossby number 27 , perhaps providing a pathway for understanding the helicity sign rule.
New helioseismic data analyses could test our predictions.The MRI would not operate if the poloidal field is too strong, nor would it explain the torsional oscillations if it is too weak.We predict correlations between the flow perturbations and magnetic fields, which time-resolved measurements could test, constraining joint helioseismic inversions of flows and magnetic fields.
An MRI-driven dynamo may also explain the formation and cessation of occasional grand minima 28 (e.g., Maunder).An "essentially nonlinear" MRI dynamo does not start from an infinitesimal seed field upon each new cycle (see Methods).Rather, a moderate poloidal field exists at the solar minimum, and the MRI processes it into a toroidal configuration.If the self-sustaining poloidal-totoroidal regeneration sometimes happens imperfectly, then subsequent solar cycles could partially fizzle, leading to weak subsurface fields and few sunspots.Eventually, noise could push the system back onto its normal cyclic behaviour, as in the El Nino Southern Oscillation 29 .
Finally, our simulations intentionally contain reduced physics to isolate the MRI as a critical agent in the dynamo process, filtering out large-scale baroclinic effects, small-scale convection and nonlinear dynamo feedback.Modelling strong turbulent processes is arduous: turbulence can simultaneously act as dissipation, drive large-scale flows like the NSSL, produce mean electromotive forces, and excite collective instabilities.While sufficiently strong turbulent dissipation could eventually erase all large-scale dynamics, the mere presence of the solar torsional oscillations implies much can persist within the roiling background.ii)   H(r, θ).The timescales t e , P represent the instability e-folding time and oscillation period, respectively.Top row: Case (i) shows a typical directly growing "fast-branch" mode with no oscillation and growth rates γ ≈ 0.06 Ω 0 .The bottom row: Case (ii) shows a typical large-scale "slow-branch" mode with a roughly five-year period.In each case, we fix the overall amplitude to 1 nHz for the rotational perturbations, with all other quantities taking their corresponding relative values.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 (e.g., the deep MRI branches at high latitudes) should employ a fully compressible (but low Mach number) model 32 .
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(θ), where Ω 0 = 466 nHz ≈ 2.92•10 −6 s −1 and Z(z) = 1 + 0.02 z − 0.01 z 2 − 0.03 z 3 (12) We use the angular fit from Howe 34 .The radial approximation results from fitting the equatorial profile from Larson and Schou 20 shown in fig.1(a).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, where and B 0 = 90 G.The r −3 term represents a global dipole.The r 2 term represents a field with a similar structure but containing electric current, The background field is in MHD force balance, The MHD force balance generates magnetic pressure, which inevitably produces entropy, s ′ , and enthalpy, h ′ , perturbations via, where and Γ 3 is the 3rd adiabatic index.However, the MRI is a weak-field instability, implying magnetic buoyancy/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 e r • B 0 = B(r) cos(θ) vanishes at r = r 1 .The poloidal field strength in the photosphere is ≈ 1 G but measurements suggest sub-surface field strengths ≈ 50-150 G 9 .The near-surface field should exhibit a strong horizontal (as opposed to radial) character.Magnetic pumping 36 via 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 ≈ 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 θ • B 0 is constant to within 8% over r 0 < r < r 1 .
However, a pure dipole varies by ≈ 37% across the domain.The RMS field amplitude is |B| RMS ≈ 2 B 0 = 180 G, about 25% larger than the maximum reported inferred dipole equivalent 9 .However, projecting our field onto a dipole template gives an ≈ 70 G equivalent at the r = r 1 equator.
Overall, the subsurface field is the least constrained input to our calculations, the details of which surely change over multiple cycles.
Model equations: Respectively, the linearised anelastic momentum, mass-continuity, and magnetic induction equations are: where the traceless strain rate To find eigenstates, ∂ t → γ + i ω, where γ is the real-valued grown rate, and ω is a real-valued oscillation frequency.The induction equation ( 23) automatically produces MRI solutions satisfying 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 , where h ′ represent enthaply perturbations 27 .
The velocity perturbations are impenetrable (u r = 0) and stress-free (σ rθ = σ rϕ = 0) at both boundaries.For the magnetic field, we enforce perfect conducting conditions at the inner boundary . 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 outcomes 35 .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 (i.e., ∂ r b r = b θ = b ϕ = 0) which, though 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 Re = Rm = U 0 L 0 /ν ≈ 1500 where U 0 ≈ 5000 cm/s is the maximum shear velocity jump across the NSSL and L 0 ≈ 0.06R ⊙ is the distance between min/max shear velocity (see "NSSL energetics and turbulence parameterisation" below).
We a posteriori compute the following scalar-potential decompositions, 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, There is no initial helicity in the background poloidal magnetic field, 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.
The only possible hemispheric contributions arise when considering quadratic mode interactions, This order is the first where we could expect a non-trivial signal.
Finally, we also solve the system using multiple different mathematically equivalent equation formulations (e.g., 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 code 26 employs general tensor calculus in sphericalpolar coordinates using spin-weighted spherical harmonics in (θ, ϕ) 37,38 .For the finite radial shell, the code uses a weighted generalised Chebyshev series with sparse representations for differentiation, radial geometric factors and non-constant coefficients (e.g., ρ 0 (r)).Since the background magnetic field and the differential rotation are axisymmetric and they only contain a few low-order separable terms in latitude and radius, these two-dimensional non-constant coefficients have a loworder representation in a joint expansion of spin-vector harmonics and Chebyshev polynomials.
The result is a two-dimensional generalised non-Hermetian eigenvalue problem A x = λ B x, where x represents the full system spectral-space state vector.The eigenvalues and eigenmodes presented here are converged to better than 1% relative absolute error (comparing 256 versus 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 https://github.com/geoffvasil/nssl_mri.
the potential from eq. ( 35) is positive everywhere in the domain.
Given the cylindrical radius, r, the local angular momentum and magnetic flux density The respective local flux transport For quadratic-order feedback, For linear meridional perturbations, For the angular components, The modes along the vertical axis appear to lie on a continuum, accumulating at a lower value.The isolate modes appear to be point spectra.The red circle represents the case (i) "fast branch" from the main text.The purple circle (with its complex conjugate) represent the case (ii) presented "slow branch".

1 -Figure 1 :
Figure 1: Measured internal solar rotation profiles.a, Heliosesimic differential rotation profile, Ω(θ, r), using publicly available data from [T.P. Larson and J. Schou.Global-mode analysis of full-disk data from the Michelson Doppler imager and the helioseismic and magnetic imager.Solar Physics, 293(2), 2018] 20 .b, c The respective latitudinal and radial shear gradients r sin(θ)∇Ω(θ, r); computed via a non-uniform 4th-order centered finite-difference scheme.The latitudinal mean of tachocline shear is ≈ 200 nHz and peak amplitudes are below ≈ 350 nHz.Conversely, the near-surface shear averages ≈ 400-600 nHz (with rapid variation in depth) and peak values over 1200 nHz.Bottom row: d, Helioseismic measurements of solar torsional oscillations.The red shows positive residual rotation rates and blue shows negative residual rotation rates after removing the 1996 annual mean of Ω(r, θ).Each slice shows the rotational perturbations 1, 2, 3 and 4 years after the approximate solar minimum.The notation "min+1Yr -min" means taking the profile at 1 year past solar minimum and subtracting the profile at solar minimum.The colour table saturates at 1 nHz, corresponding to about ≈ 400 cm/s surface flow amplitude.Further contour lines show 1 nHz increments within the saturated regions.Diagram in d reprinted from figure 2 of [S.V. Vorontsov, J. Christensen-Dalsgaard, J. Schou, V. N. Strakhov, and M. J. Thompson.Helioseismic measurement of solar torsional oscillations.Science, 296(5565):101-103, April 2002.] 3 with permission from AAAS.
The matrices, A and B, are spectralcoefficient representations of the relevant linear differential and multiplication operators.Cases (i) & (ii) 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 .

∂Figure 3 :
Figure 3: Full diagram of complex-valued eigen-spectrum.The time dependence sends ∂ t → γ +iω, with each real/imaginary component measured in terms of Ω 0 = 466 nHz.The modes along the vertical axis appear to lie on a continuum, accumulating at a lower value.The isolate modes appear to be point spectra.The red circle represents the case (i) "fast branch" from the main text.The purple circle (with its complex conjugate) represent the case (ii) presented "slow branch".

Figure 4 :
Figure 4: The complete collection of "fast branch" modes.The growth rates correspond to the vertical axis of Extended Data fig. 1.Each case contains no discernible oscillations.For completeness, we show (boxed in grey) the t e = 60 day fast-branch case (i) presented in the main text fig.2(a)