¹

1 The Solar Tachocline↩︎

The solar tachocline is a region of primarily radial shear at the base of the solar convection zone (CZ), where strong latitudinal differential rotation transitions to nearly solid-body rotation in the underlying radiative zone (RZ). The tachocline is observed helioseismically to be centered at \(r_{t,\odot}\approx0.69R_\odot\) (which roughly coincides with the base of the CZ) and to have a thickness of \(\Gamma_\odot\lesssim0.05R_\odot\) (\(\Gamma_\odot\) is too small to be helioseismically resolved, implying that it has an upper bound roughly equal to the helioseismic inversion kernel width; e.g., [1]). Some measurements estimate a wider tachocline (\(\Gamma_\odot\lesssim0.10R_\odot\); e.g., [2], [3]) or a narrower tachocline (\(\Gamma_\odot\lesssim0.02R_\odot\); e.g., [4], [5]).

Regardless of the true tachocline thickness, even the most liberal estimates for \(\Gamma_\odot\) pose a major dynamical problem for solar physics. It is hypothesized [6] that the CZ’s differential rotation should spread into the RZ by a process similar to circulation “burrowing" in rotating stably stratified shear flows (e.g., [7], [8]), thus widening the tachocline. A shear flow in a rotating system (i.e., differential rotation) is usually accompanied by a horizontal temperature gradient due to thermal wind balance (e.g., [9], [10]). This gradient tends to spread (burrow) further into the stable layer via thermal conduction, carrying with it the circulation and differential rotation associated with the thermal wind. In the Sun, the dominant thermal diffusion is radiative, and [6] showed that burrowing (now referred to as”radiative spread") should have increased \(\Gamma_\odot\) to \(\sim\) \(0.4R_\odot\) by the current age of the Sun.

[6]’s original argument that the solar tachocline should radiatively spread assumes axisymmetry and linearized fluid equations. Under those conditions, radiative spread occurs “hyperdiffusively" (governed by \(\nabla^4\) instead of \(\nabla^2\)) on the solar Eddington-Sweet time \({P_{ {\rm ES}, \odot}}\). In the hyperdiffusive case, \(\Gamma(t)/\Gamma_\odot\sim(t/{P_{ {\rm ES}, \odot}})^{1/4}\), where \(\Gamma(t)\) is the time-dependent tachocline thickness and \(t\) is the time since initial confinement [i.e., \(\Gamma(0)=\Gamma_\odot\)]. Since the Eddington-Sweet time is so long for the Sun (\({P_{ {\rm ES}, \odot}}\approx600\) Gyr; see Table 6), this hyperdiffusive property is essential for the burrowing to be significant on time-scales as small as the solar age (\(\sim\)​5 Gyr). Recent 3D fully nonlinear simulations have shown that circulation burrowing does indeed occur in more realistic settings (as long as the time-scales are properly ordered; see [11], [12]). But whether realistic solar burrowing would be hyperdiffusive is still an open question and requires further investigation.

If circulation burrowing is indeed significant for the Sun, it is obvious that there must be a confining (or “rigidifying") torque in the RZ to keep \(\Gamma_\odot\) under the helioseismically constrained upper bound. There are currently two dominant tachocline confinement scenarios that postulate the origin of this torque. The first, proposed by [6], is essentially hydrodynamic. It is supposed that hydrodynamic shear instabilities associated with the differential rotation create turbulence with predominantly horizontal motion, owing to the strong convectively stable stratification of the RZ. The Reynolds stresses from this horizontal turbulence then act like an enhanced horizontal viscosity, causing preferentially horizontal angular momentum transport, thereby eliminating any burrowing shear on the relatively fast time-scale of months to years. Hence, this scenario is often also called the”fast confinement scenario" (e.g., [13], [14]).

However, similar horizontal turbulence in the Earth’s stratosphere is theorized to be “anti-diffusive," that is, transporting angular velocity up the rotation gradient instead of down it and driving the system away from solid-body rotation (e.g., [15], [16]). In any event, angular momentum transport by stratified turbulence in a solar-like system is likely more complicated than simply”diffusive or anti-diffusive." For example, [17] argue that horizontal turbulence in the presence of a weak toroidal magnetic field creates Maxwell stresses that nearly exactly cancel the Reynolds stresses, yielding zero net momentum transport. Finally, it remains unclear exactly how anisotropic stratified turbulent transport really is. For example, recent 3D direct numerical simulations [18], [19] show that meanders of the streamwise flow (in a sufficiently turbulent regime) can vary on small vertical length-scales until secondary vertical shear instabilities (and associated vertical momentum transport) develop.

[20] proposed an alternative, magnetic confinement scenario. They argued that a weak (minimum \(\sim\)​1 G) poloidal magnetic field in the RZ could resist the shearing motion of any imposed differential rotation via magnetic tension. This magnetic torque would be generated on the time-scale of radiative spread, namely some fraction of \({P_{ {\rm ES}, \odot}}\). Hence, [20]’s scenario is sometimes called the “slow confinement scenario." Note that the fast confinement scenario is mostly hydrodynamic (with magnetism possibly playing a secondary role in modifying the primary baroclinic and shear instabilities), while the slow confinement scenario is fundamentally magnetic.

Finally, a “fast magnetic confinement scenario" has been proposed and modeled in 1D (e.g., [21], [22]). Here, source of the RZ’s confining poloidal field is the cycling solar dynamo (with the cycle period of \(\sim\)​11 yr, i.e., fast compared to \({P_{ {\rm ES}, \odot}}\), but slow compared to time-scales associated with most hydrodynamic instabilities) diffusing downward to a skin depth.

In prior global simulations of solar-like CZ–RZ systems, the chosen parameters have made radiative spread insignificant on the time-scales the simulations can be run. Nevertheless, significant viscous spread occurs (see Section 4.2) and simulated tachoclines have been confined against this viscous spread through a variety of mechanisms. [23] used combined mechanical and thermal forcing to explicitly impose a steady-state tachocline in a simulation using the ASH code. Further simulations using the ASH and Rayleigh codes—which are direct numerical simulation (DNS) codes—have implemented temporary, slowly spreading tachoclines through significantly lowered values of the viscosity in the RZ compared to the CZ (e.g., [24]–[26]). Finally, the implicit large-eddy simulation (ILES) code EULAG ensures very small effective numerical viscosity in stable regions owing to the nature of the ILES time-stepping algorithm MPDATA [27]. On the time-scales for which EULAG simulations are run, both viscous and radiative spread are thus negligible and tachoclines that are effectively steady can occur (in both magnetic and purely hydrodynamic cases) without an explicit confinement mechanism being necessary (e.g., [28], [29]).

We recently presented ([30]; hereafter ; see also [31]), the first 3D, spherical-shell simulation (in our case, a DNS) to achieve a steady-state tachocline that was self-consistently confined against explicit viscous spread. The source of the confinement was magnetic torque, which was in turn generated by a non-axisymmetric, quasi-periodic dynamo. In the CZ, the magnetism was topologically similar to the “partial wreaths" (longitudinally elongated bands of intense toroidal magnetism, with alternating polarity in longitude) identified in our prior CZ-only dynamos [32], [33]. We showed in [32] that the partial wreaths in the CZ-only case tended to form a long-lasting magnetic structure that more or less rotated rigidly in a preferred frame. In the combined CZ–RZ tachocline systems considered in the current work, the partial wreaths rotate with respect to the RZ below. As far as the rigidly-rotating RZ is concerned, the partial wreaths above resemble a periodically reversing poloidal field and therefore the field diffusively imprints from the overshoot layer to a depth in the RZ consistent with the electromagnetic skin effect.

The main conclusion of the present paper is that the confinement mechanism identified in can be regarded as a more general version of the fast magnetic confinement scenario that stays robust in a wider parameter space (containing multiple cycling frequency components of the dynamo) and in a 3D geometry with a fully coupled CZ and RZ. Furthermore, a rotating, large-scale non-axisymmetric poloidal field structure takes the place of the reversing axisymmetric magnetism (“full wreaths") typically invoked in connection with the observed solar cycle, or magnetic butterfly diagram. Our evidence consists of a family of solutions related to the one from , but with a range of magnetic Prandtl numbers \({\rm{Pr_m}}\). One key effect of varying \({\rm{Pr_m}}\) (while keeping the other control parameters fixed) is to achieve a range of magnetic field strengths in the saturated dynamo state, while keeping other key diagnostic parameters (like the Reynolds and Rossby numbers) relatively unchanged. The rest of this paper is structured as follows. In Section 2, we describe our equation set and control parameters. In Section 3, we describe the three solution regimes (weak-, medium-, and strong-field) that our dynamos achieve. In Section 4, we present the degree of tachocline confinement, as well as the associated torque balance, for our simulations. In Section 5, we describe the two distinct types of magnetic cycle exemplified by the weak- and strong-field regimes. In Section 6, we show that for all cases, the poloidal magnetic field strength in the RZ is consistent with diffusive imprinting of the CZ’s poloidal field according to the electromagnetic skin effect. In Section 7, we highlight the distinctions between axisymmetric and non-axisymmetric polarity reversals. Finally, in Section 8, we discuss our results in the context of the solar tachocline confinement problem.

2 Numerical Scheme & Simulation Parameters↩︎

We evolve the 3D magnetohydrodynamic (MHD) equations in spherical shells using the open-source Rayleigh code [34]–[36]. We make use of both spherical coordinates [\(r\) (radius), \(\theta\) (colatitude), and \(\phi\) (azimuth angle)] and cylindrical coordinates [\(\lambda=r\sin\theta\) (cylindrical radius), \(\phi\) (azimuth angle), and \(z=r\cos\theta\) (axial coordinate)]. The symbol \(\hat{\boldsymbol{e}}\) denotes a unit vector. The equations are solved in a frame rotating with the constant angular velocity \(\boldsymbol{\Omega}_0=\Omega_0\hat{\boldsymbol{e}}_z\). The Coriolis force is kept but the oblateness and centrifugal force are ignored. Each shell extends from an inner radius \(r_{\rm{in}}\) to an outer radius \(r_{\rm{out}}\). We divide the shell into two layers of equal depth, separated at \(r_0\equiv(r_{\rm{in}}+r_{\rm{out}})/2\). The top half (\(r_0\) to \(r_{\rm{out}}\); the CZ) is nominally convectively unstable and the bottom half (\(r_{\rm{in}}\) to \(r_0\); the RZ) convectively stable. Rayleigh solves the anelastic MHD equations, which allow significant density contrast across the shell, but disallow sound waves (e.g., [37]–[40]). The anelastic approximation consists of assuming a solenoidal mass flux [see Equation 1 ] and thermodynamic perturbations that are small relative to a well-chosen “background" or”reference" state. In Rayleigh, the background state is always spherically symmetric and time-independent (e.g., [34]). We choose a background entropy gradient \(d\overline{S}/dr\) that changes from stable to unstable near \(r=r_0\) over the transition width \(\delta\) and a gravitational acceleration \(\overline{g}=GM_\odot/r^2\) (where \(G=6.67\times 10^{-8}\;\rm{cm^3\;g^{-1}\;s^{-2}}\) is the universal gravitational constant and \(M_\odot=1.99\times 10^{33}\;\rm g\) the solar mass). If we further assume a hydrostatic, ideal gas [with constant specific heats \(c_{\rm{v}}\) (at constant volume) and \(c_{\rm{p}}\) (at constant pressure)], the choices for \(d\overline{S}/dr\) and \(\overline{g}\) determine the background density \(\overline{\rho}\), temperature \(\overline{T}\), and squared buoyancy frequency \(\overline{N^2}\equiv(\overline{g}/c_{\rm{p}})d\overline{S}/dr\) (we use \(\overline{N^2}\) in favor of \(d\overline{S}/dr\) in the equations).

We choose all diffusivities (kinematic viscosity \(\overline{\nu}\), thermal diffusivity \(\overline{\kappa}\), and magnetic diffusivity \(\overline{\eta}\)) to increase with height like \(1/\overline{\rho}^{1/2}\). We choose an internal heating function \(\overline{Q}\) (representing radiative heating from below) that deposits thermal energy preferentially in roughly the bottom third of the CZ and drives convection. In the RZ, we set \(\overline{Q}=0\), tapered from its profile in the CZ over a width \(\delta_{\rm{heat}}\). We fully describe our reference state in Appendix 9 and its analogy to the Sun in Appendix 10.

The dimensional equations of motion are \[\begin{align} \nabla\cdot(\overline{\rho}\boldsymbol{u}) &= 0, \label{eq:contdim} \end{align}\tag{1}\] \[\begin{align} \nabla\cdot\boldsymbol{B} &= 0, \label{eq:divb0dim} \end{align}\tag{2}\] \[\begin{align} \overline{\rho}\left(\frac{D\boldsymbol{u}}{Dt} \right) &= -2\overline{\rho}\boldsymbol{\Omega}_0\times\boldsymbol{u}-\overline{\rho}\nabla \left(\frac{P}{\overline{\rho}}\right) +\frac{\overline{\rho}\,\overline{g}S}{c_{\rm{p}}}\hat{\boldsymbol{e}}_r\nonumber\\ &\;\;\;+ \nabla\cdot \boldsymbol{D} + \frac{1}{\mu}(\nabla\times\boldsymbol{B})\times\boldsymbol{B}\nonumber\\\tag{3}\\ \text{where}\ \ \ \ \ D_{ij} &\equiv 2\overline{\rho}\,\overline{\nu}\left[e_{ij} - \frac{1}{3}(\nabla\cdot\boldsymbol{u}) \delta_{ij} \right]\tag{4}\\ \text{and}\ \ \ \ \ e_{ij} &\equiv\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)\tag{5}, \end{align}\] \[\begin{align} \overline{\rho}\overline{T}\left(\frac{DS}{Dt} \right) =\;&Q - \overline{\rho}\overline{T}\frac{d\overline{S}}{dr}u_r + \nabla\cdot\left(\overline{\kappa}\,\overline{\rho}\overline{T}\nabla S\right]\nonumber \\ &+ D_{ij}e_{ij} + \frac{\eta}{4\pi}|\nabla\times\boldsymbol{B}|^2, \label{eq:endim} \end{align}\tag{6}\] and \[\begin{align} \frac{\partial\boldsymbol{B}}{\partial t} =\;&\nabla\times(\boldsymbol{u}\times\boldsymbol{B} - \eta\nabla\times\boldsymbol{B}). \label{eq:inddim} \end{align}\tag{7}\] Here, \(D/Dt\equiv\partial/\partial t+\boldsymbol{u}\cdot\nabla\) is the material derivative and \(\mu\) the vacuum permeability (\(\mu=4\pi\) in Gaussian units).

Rayleigh was originally run by solving these dimensional equations. In this work, however, we discuss only the equivalent non-dimensional simulations. Length is scaled by the CZ (or RZ) thickness \(H\equiv(r_{\rm{out}}-r_{\rm{in}})/2\) and time by the rotational time-scale \(\Omega_0^{-1}\). The velocity \(\boldsymbol{u}\) is scaled by \([\boldsymbol{u}]\equiv\Omega_0H\) and the vorticity \(\boldsymbol{\omega}\equiv\nabla\times\boldsymbol{u}\) by \(\Omega_0\) (we use square brackets to denote the unit of each fluid variable). Each background-state profile is scaled by its volume-average over the CZ (denoted by a tilde, e.g., \(\tilde{\rho}\)), except for \(\overline{N^2}\), which is scaled by its volume-average over the RZ (denoted by \(\langle \overline{N^2}\rangle_{\rm{RZ}}\)), and \(\overline{Q}\), which is scaled as described below. The pressure perturbation \(P\) is scaled by \([P]\equiv\tilde{\rho}(\Omega_0H)^2\) and the magnetic field \(\boldsymbol{B}\) by \([\boldsymbol{B}]\equiv\sqrt{\mu\tilde{\rho}}(\Omega_0H)\).

As noted by [41], the chosen non-dimensionalization omits the diffusivities from the scales for time and the magnetic field. This is helpful in extending scaling relationships to stellar regimes, where diffusive effects are not believed to play a large role (although we note at the outset that such scaling relationships are likely not present in this work, where diffusive effects do play a large role). An added benefit of this non-dimensionalization is that \(\boldsymbol{u}\) and \(\boldsymbol{B}\) appear with order-unity coefficients in the momentum equation, so their relative importance (to both the force balance and the partition of kinetic and magnetic energy) can be inferred directly from their non-dimensional values.

The internal heating \(\overline{Q}\), coupled with the thermal boundary conditions described below, drives convection by establishing sharp entropy gradients in a thermal boundary layer near the top of the CZ (e.g., [34], [42]). This convection (and conduction, especially in the boundary layer), must carry a “non-radiative" energy flux \(\overline{F_{\rm nr}}\equiv(1/r^2)\int_{r_0}^r\overline{Q}(x)x^2dx\) in the statistically steady state. The entropy perturbation \(S\) is thus scaled by its estimated difference across the conductive boundary layer (\([S]=\Delta S\equiv\widetilde{F_{\rm nr}}H/\tilde{\rho}\tilde{T}\tilde{\kappa}\)) and \(\overline{Q}\) by \(\widetilde{F_{\rm nr}}/H\).

With these scaling choices, the non-dimensional equations of motion are \[\begin{align} \label{eq:mom} \overline{\rho}\left(\frac{D\boldsymbol{u}}{Dt}\right) &= -2\overline{\rho}\hat{\boldsymbol{e}}_z\times\boldsymbol{u}-\overline{\rho}\nabla\left(\frac{P}{\overline{\rho}} \right) +{\rm{Ra}}_{\rm{F}}^*\overline{\rho}\, \overline{g}S\hat{\boldsymbol{e}}_r\nonumber\\ &\;\;\;+{\rm{Ek}}\nabla\cdot\boldsymbol{D}+(\nabla\times\boldsymbol{B})\times\boldsymbol{B} \end{align}\tag{8}\] \[\begin{align} \label{eq:heat} \overline{\rho}\overline{T}\frac{DS}{Dt} = &\frac{{\rm{Ek}}}{{\rm{Pr}}} \overline{Q}- \frac{{\rm{Bu}}}{{\rm{Ra}}_{\rm{F}}^*} \overline{\rho}\overline{T}\frac{\overline{N^2}}{\overline{g}} u_r+ \frac{{\rm{Ek}}}{{\rm{Pr}}} \nabla\cdot(\overline{\rho}\overline{T}\overline{\kappa}\nabla S) \nonumber\\ &+ \frac{{\rm{Di}}{\rm{Ek}}}{{\rm{Ra}}_{\rm{F}}^*} D_{ij}e_{ij} + \frac{{\rm{Di}}{\rm{Ek}}}{{\rm{Pr_m}}{\rm{Ra}}_{\rm{F}}^*} \overline{\eta}|\nabla\times\boldsymbol{B}|^2, \end{align}\tag{9}\] \[\begin{align} \label{eq:ind} \text{and}\ \ \ \ \ \frac{\partial\boldsymbol{B}}{\partial t} = \nabla\times(\boldsymbol{u}\times\boldsymbol{B}) - \frac{{\rm{Ek}}}{{\rm{Pr_m}}} \nabla\times(\overline{\eta}\nabla\times\boldsymbol{B}). \end{align}\tag{10}\] Here, Equations 1 and 2 still apply and are unchanged, and \(D_{ij}\) and \(e_{ij}\) are defined exactly as in Equations 4 and 5 , respectively. All field variables (\(\boldsymbol{u}\), \(\boldsymbol{B}\), \(S\), and \(P\)), spatial quantities (\(r\), \(t\), \(\lambda\), \(z\), and \(\nabla\)), and background-state profiles now denote their non-dimensional values. The non-dimensional input numbers (definitions and values) are given in Table 1.

The reference-state control parameters are the ratio of specific heats \(\gamma\), the CZ-to-RZ aspect ratio \(\alpha\), the CZ aspect ratio \(\beta\), the number of scale heights across the CZ \(N_\rho\), and the transition widths \(\delta\) and \(\delta_{\rm{heat}}\). This reference state (except for the diffusivity profiles) is reasonably solar-like and describes the upper 2.1 density scale-heights of the solar RZ and the lower 3 density scale-heights of the solar CZ (see Appendix 9). In units of \(H\), the non-dimensional solar radius is \(R_\odot=4.39\) (see Table 6). We plot radial profiles as functions of \(r/R_\odot\), to more easily compare to prior work.

The fluid control parameters are the Prandtl number \({\rm{Pr}}\), the magnetic Prandtl number \({\rm{Pr_m}}\), the modified Rayleigh number \({\rm{Ra}}_{\rm{F}}^*\), the Ekman number \({\rm{Ek}}\), and the buoyancy number \(\rm{Bu}\). The dissipation number \({\rm{Di}}\equiv\tilde{g}H/(c_{\rm{p}}\tilde{T})=1.72\) for the cases here and is not a control parameter in our convention, being a function of \(\gamma\), \(\beta\), and \(N_\rho\), which we deem reference state control parameters (see Appendix 9 and [43]). Some additional parameters (that can be derived from the input parameters given in Table 1) are given in Table 5.

Equations 1 –7 are discretized in space. For all simulations, we use three sets of stacked Chebyshev collocation points in \(r\) (\(N_r/3=64\) points in each domain), \(N_\theta=384\) Legendre collocation points in \(\theta\), and \(N_\phi=2N_\theta=768\) uniformly spaced collocation points in \(\phi\). The Chebyshev points cluster near each domain’s boundaries. We require increased resolution in the overshoot layer (i.e., in the vicinity of \(r_0/R_\odot=0.719\)), and so we set the radial domain boundaries to lie at \(r/R_\odot=\{0.491, 0.669, 0.719, 0.947\}\) (or equivalently, \(r-r_0=\{-1.000, -0.219, 0.000, 1.000\}\)). Nonlinear terms and the Coriolis force are evaluated in physical space (i.e., on the discretized spatial grid), while the remaining linear terms are evaluated in spectral space, using Chebyshev polynomials in each \(r\) sub-domain and spherical harmonics in \(\theta\) and \(\phi\). The variables in physical space are de-aliased using the 2/3 rule: the maximum Chebyshev degree (in each \(r\) sub-domain) is \(n_{\rm{max}}=42\) and the maximum spherical harmonic degree is \({\ell_{\rm{max}}}=255\). For more details, see [44] and [40], who pioneered Rayleigh’s pseudo-spectral algorithm.

Each magnetic simulation differs only in the choice of \({\rm{Pr_m}}\), which ranges from 1 to 8. We also consider a purely hydrodynamic simulation (referred to as “Case H"), which has all the parameters listed in Table 1, but no magnetic field. We refer to each magnetic case by its value of \({\rm{Pr_m}}\) rounded to two decimal places: e.g.,”Case 1.08" means \({\rm{Pr_m}}=1.076\). Cases H and 4.00 were analyzed in . All chosen values of \({\rm{Pr_m}}\) are listed in Table 2.

At both boundaries, we use stress-free and impenetrable conditions on \(\boldsymbol{u}\), potential-field-matching conditions on \(\boldsymbol{B}\), and fixed-entropy-gradient conditions on \(S\). Specifically, we set \(\partial S/\partial r\) to zero at the bottom boundary (thus allowing no conductive flux in or out) and set it to a latitudinally independent negative value at the top boundary, such that the energy conducted out the top is equal to the energy injected by \(\overline{Q}\) (e.g., [42]). The convection is initialized by introducing weak noise in \(S\) (amplitude \(\sim\) \(10^{-3}\)), randomly distributed in space throughout the entire shell. For the magnetic cases, we further introduce weak noise in \(\boldsymbol{B}\) (amplitude \(\sim\) \(10^{-6}\)), randomly distributed in space throughout the CZ only. The other field variables (\(\boldsymbol{u}\) and \(P\)) are initialized to zero in all space.

We use several types of averages in this work. Let \(\psi=\psi(r,\theta,\phi,t)\) denote a scalar quantity (or a single component of a vector quantity) dependent on position and time. Then \(\left\langle\psi\right\rangle_{\phi}\), \(\left\langle\psi\right\rangle_{\rm{sph}}\), \(\left\langle\psi\right\rangle_{\rm{CZ}}\), \(\left\langle\psi\right\rangle_{\rm{RZ}}\), and \(\left\langle\psi\right\rangle_{\rm{full}}\) denote instantaneous averages of \(\psi\) over longitude, spherical surfaces, the CZ (volume-average from \(r_0\) to \(r_{\rm{out}}\)), the RZ (volume-average from \(r_{\rm{in}}\) to \(r_0\)), and the full shell (volume-average from \(r_{\rm{in}}\) to \(r_{\rm{out}}\)), respectively. An additional temporal average (over the “equilibrated state"; see the following section) is denoted by appending a”\(t\)" to the subscript in the average: e.g., \(\left\langle\psi\right\rangle_{\phi,t}\). Subtracting the instantaneous longitudinal average is denoted by a prime: \(\psi^\prime\equiv\psi-\left\langle\psi\right\rangle_{\phi}\). We also colloquially refer to \(\left\langle\psi\right\rangle_{\phi}\) and \(\psi^\prime\) as the “mean and fluctuating" components of \(\psi\), respectively.

3 Dynamo Regimes↩︎

All the magnetic cases presented here yield sustained large-scale dynamos. As convection and dynamo action become significant, the field variables grow from their initially small values to amplitudes of order unity. We quantify this growth in terms of the kinetic energy density of the differential rotation, \(\rm KE_{DR}\), and the magnetic energy density, \(\rm ME\): \[\begin{align} \label{eq:drke95and95me} {\rm KE_{DR}} \equiv\frac{1}{2}\overline{\rho}\left\langle u_\phi\right\rangle_{\phi}^2 \ \ \ \ \ \text{and}\ \ \ \ \ {\rm ME}\equiv\frac{1}{2}\left\langle\boldsymbol{B}^2\right\rangle_{\phi}. \end{align}\tag{11}\]

In Equation 11 , the energy densities are functions of \(r\), \(\theta\), and \(t\), which we average further in the subsequent analysis. In Figure 1, we show the growth and long-term behavior of the full-shell-averaged energy densities for some representative simulations. After a certain time (which we call \(t={t_{\rm{eq}}}\)), the system achieves a “statistically steady" or”equilibrated" state, in which the volume-averaged magnitude of each field variable fluctuates about a well-defined temporal mean. We choose \({t_{\rm{eq}}}\) (fairly roughly) by eye from plots like Figure 1 (a). For example, we choose \({t_{\rm{eq}}}=2000{P_{\rm{rot}}}\) for Case 1.06, \({t_{\rm{eq}}}=1000{P_{\rm{rot}}}\) for Case 2.00, and \({t_{\rm{eq}}}=600{P_{\rm{rot}}}\) for Case 4.00 (see Table 2 for all values of \({t_{\rm{eq}}}\)).

Figure 1 (a) suggests three basic dynamo regimes. The low-\({\rm{Pr_m}}\) solution (Case 1.06; blue curves) lies in a “weak-field regime", characterized by \(\left\langle\rm ME\right\rangle_{\rm{full}}\) always being orders of magnitude weaker than \(\left\langle\rm KE_{DR}\right\rangle_{\rm{full}}\). There is a regular magnetic energy cycle, with a period of roughly \(750\;{P_{\rm{rot}}}\). The high-\({\rm{Pr_m}}\) solution (Case 4.00; red curves) lies in a”strong-field regime", characterized by \(\left\langle\rm ME\right\rangle_{\rm{full}}\) about in equipartition with \(\left\langle\rm KE_{DR}\right\rangle_{\rm{full}}\), while \(\left\langle\rm KE_{DR}\right\rangle_{\rm{full}}\) itself is much weaker than in the weak-field case. Finally, the intermediate-\({\rm{Pr_m}}\) solution (Case 2.00; orange curves) lies in a “medium-field regime." Case 2.00 has properties similar to those of a strong-field dynamo some of the time, but occasionally \(\left\langle\rm ME\right\rangle_{\rm{full}}\) falls below its strong-field value and then \(\left\langle\rm KE_{DR}\right\rangle_{\rm{full}}\) steadily increases above its strong-field value (representing an increase in differential rotation or partial disappearance of the tachocline). After \(\left\langle\rm KE_{DR}\right\rangle_{\rm{full}}\) reaches a critical level, \(\left\langle\rm ME\right\rangle_{\rm{full}}\) grows rapidly, lowering \(\left\langle\rm KE_{DR}\right\rangle_{\rm{full}}\) back to its lower, strong-field value. There is no clear cycling behavior, obvious physical trigger, or general predictability for the medium-field cases’ temporary epochs of strong differential rotation. Figure 1 (b) shows the equilibrated levels of the kinetic energy in the differential rotation, \(\left\langle\rm KE_{DR}\right\rangle_{{\rm full},t}\), and the magnetic energy, \(\left\langle\rm ME\right\rangle_{{\rm full},t}\), for all simulations. The weak-field regime (for which the differential rotation has the same magnitude as in Case H) sits in the narrow range of roughly \(1.00\lesssim{\rm{Pr_m}}\lesssim1.06\). The medium-field regime (for which the differential rotation is substantially weakened compared to the weak-field cases but intermittently becomes stronger) occupies roughly \(1.08\lesssim{\rm{Pr_m}}\lesssim2.5\). The strong-field regime (lowest differential rotation and highest magnetic energy) occupies \({\rm{Pr_m}}\gtrsim2.5\). Note that these identified regimes and their boundaries are only suggestive, given our limited resolution in \({\rm{Pr_m}}\)-space.

The weak-field dynamos tend to be more axisymmetric (magnetism dominated by azimuthal wavenumber \(m=0\)) than the strong-field dynamos. Figure 2 shows the toroidal magnetic field projected on spherical surfaces for four solutions at different \({\rm{Pr_m}}\) and therefore in the different regimes. For Case 1.00 (the weak-field regime), there is a strong \(m=0\) component, both in the CZ and even more so in the RZ. For higher values of \({\rm{Pr_m}}\) (the medium- and strong-field regimes), the field in the CZ becomes increasingly dominated by small scales, but retains a large-scale (\(m=0,1,2\)) envelope. For all cases, the RZ appears to act as a low-pass filter for the spatial scales of the field, letting only the low \(m\)’s survive. This is especially apparent for Case 8.00 [Figures 2 (g,h)].

As a whole, the dominant \(m=0\) field structures at low \({\rm{Pr_m}}\) resemble what has previously been called magnetic “wreaths"—toroidal bands of strong magnetism looping the full sphere in a given hemisphere (e.g., [45]–[47]). Such wreaths are often invoked in connection to the magnetic butterfly diagram, as interior reservoirs of toroidal field from which smaller loops can potentially break off and buoyantly rise to form sunspot pairs (e.g., [48]–[52]).

For higher \({\rm{Pr_m}}\), the strong \(m=1,2\) components of \(\boldsymbol{B}\) resemble the “partial wreaths" discussed at some length by [32], [33]. The cases from that work contained a dominant \(m=1\) field structure that appeared to be two opposite-polarity full wreaths tilted into each other, or possibly linked. On a spherical slice, the tilted full wreaths showed up as two opposite-polarity”partial wreaths", extending in longitude by about \(180^\circ\), with central longitudes on opposite sides of the sphere. We should note that “partial wreath" is really a placeholder for lack of a better term. The 3D field-line tracings of these non-axisymmetric structures tend to be quite difficult to interpret and it remains unclear exactly what topology (linked wreaths, tilted wreaths, or even open field-lines) is leading to the two-dimensional projections shown in Figure 2.

To quantify the non-axisymmetry in our dynamos more precisely, we partition the magnetic energy according to \(m\)-value. We define the \(m\)-component of \(\boldsymbol{B}\) through \[\tag{12} \begin{align} \boldsymbol{B}_m &\equiv \left\langle\boldsymbol{B}e^{-im\phi}\right\rangle_{\phi}\tag{13}\\ \text{or}\ \ \ \ \ \boldsymbol{B}&\equiv \sum_m \boldsymbol{B}_m e^{im\phi}.\tag{14} \end{align}\] Note that since \(\left\langle\cdots\right\rangle_{\phi}\) is computed by averaging over the uniformly spaced \(\phi\)-grid, Equation 12 represents the forward and inverse discrete Fourier transforms (DFTs) in \(\phi\). For the chosen normalization, Parseval’s theorem takes the form \[\begin{align} \label{eq:parceval} \left\langle\boldsymbol{B}^2\right\rangle_{\phi} = \sum_m |\boldsymbol{B}_m|^2. \end{align}\tag{15}\] Note that by definition, \(\boldsymbol{B}_0=\left\langle\boldsymbol{B}\right\rangle_{\phi}\). Different components of the magnetic energy can thus be attributed to different \(m\)-components of \(\boldsymbol{B}\).

Table 2 shows some basic simulation properties (\({\rm{Pr_m}}\), regime, total run time \({t_{\rm{max}}}\), and equilibration time \({t_{\rm{eq}}}\)), as well as the partitions of kinetic and magnetic energy for each simulation.² We define the magnetic energy in the “small-scale fields" by \({\rm ME_{\geq3}}\equiv(1/2)\sum_{|m|\geq3}|\boldsymbol{B}_m|^2\). With increasing \({\rm{Pr_m}}\), the fraction of magnetic energy in the small-scale fields (\(\rm ME_{\geq3}/ME\)) increases, as might be expected from the more prominent small-scale structures seen in Figure 2 at higher \({\rm{Pr_m}}\). The deficit, i.e., the fraction of energy in the large-scale fields (\(m=0,1,2\)) decreases, but the partition between each \(m\)-component is complex.

For the weak-field cases, in both the CZ and RZ, the power in the axisymmetric field dominates over the power in the \(m=1,2\) components. But for the medium- and strong-field cases, there seems to be no general rule for whether \(m=0,1,\) or \(2\) dominates. One robust feature is that in all regimes, the magnetic energy in the RZ is stored primarily in the large-scale fields. Even for the strongest-field Case 8.00, the small-scale (\(|m|\geq3\)) field components account for only \(\sim\)​30% of the magnetic energy in the RZ. These results strengthen the earlier idea that the RZ acts as a low-pass filter, letting in only the lowest-\(m\) components of the field. Such behavior is expected if the field evolution in the RZ is primarily governed by diffusion, an idea we return to in Section 6, where we discuss the skin-effect behavior of the poloidal field.

4 Tachocline Confinement↩︎

4.1 Tachocline Appearance↩︎

All dynamos in the medium- and strong-field regimes sustain steady-state tachoclines. To describe them, we define the rotation rate \(\Omega\) (as measured in the rotating frame): \[\begin{align} \label{eq:rotrate} \Omega(r,\theta)\equiv \frac{\left\langle u_\phi\right\rangle_{\phi,t}}{\lambda}. \end{align}\tag{16}\] Figure 3 shows \(\Omega\) for a weak-, medium-, and strong-field case, as well as the Sun. The weak-field case does not have a tachocline and has rotation rate nearly identical to that of Case H, that is, there is strong latitudinal rotation contrast in the CZ (\(\Delta\Omega_{\rm{CZ}}\sim0.2\), similar to the solar value), which imprints throughout the entire RZ (for Case H’s rotation profile, see ).

By contrast, the medium- and strong-field cases all have tachoclines, that is, there is (weak) differential rotation in the CZ, but nearly solid-body rotation in the RZ. The radial transition from differential to solid-body rotation appears to be quite abrupt at most latitudes from the color plots in Figures 3 (b,c), suggesting thin simulated tachoclines. However, this visual abruptness is at least partly due to our choice of bi-linear colormap for the asymmetric values of \(\Omega\) about 0. This colormap deepens the blue tones in the CZ at high latitudes. As we now demonstrate, the simulated tachoclines (after they are fit systematically) are in fact about twice as thick as the solar one.

We define the radially varying latitudinal differential rotation contrast \(\Delta\Omega(r)\): \[\begin{align} \label{eq:drotrate} \Delta\Omega(r)\equiv \Omega(r,\pi/2) - \frac{1}{2}[\Omega(r,\pi/6) + \Omega(r,5\pi/6)], \end{align}\tag{17}\] i.e., the difference in rotation rate between the equator and the average rate of 60\(^\circ\) latitude north and south. We define the rotation contrasts in the CZ and RZ, and their ratio: \[\label{eq:dom} \begin{align} \Delta\Omega_{\rm{CZ}}&\equiv\left\langle\Delta\Omega\right\rangle_{\rm{CZ}},\\ \Delta\Omega_{\rm{RZ}}&\equiv\left\langle\Delta\Omega\right\rangle_{\rm{RZ}},\\ \text{and}\ \ \ \ \ f& \equiv\frac{ \Delta\Omega_{\rm{RZ}}}{ \Delta\Omega_{\rm{CZ}}}, \end{align}\tag{18}\] respectively. Because the medium- and strong-field RZs often rotate like solid bodies, we define the (volume-averaged) constant rotation rate of the RZ: \[\begin{align} \label{eq:omrz} \Omega_{\rm{RZ}}\equiv\left\langle\Omega\right\rangle_{\rm{RZ}} \end{align}\tag{19}\]

Figure 4 shows line plots of the rotation rate along radial lines for the Sun and Case 4.00. Clearly in the solar case, the tachocline is confined to a relatively narrow radial layer, with strong differential rotation in most of the CZ and very little in the RZ, indicating a large \(\Delta\Omega_{\rm{CZ}}\), a small \(\Delta\Omega_{\rm{RZ}}\), and therefore small \(f\) for the true solar case. By contrast, in Case 4.00, while \(\Delta\Omega_{\rm{RZ}}\) is severely diminished, thus indicating that the RZ rotates nearly as a solid body, \(\Delta\Omega_{\rm{CZ}}\) is also significantly diminished.

A further deviation from the solar case is that our simulated rotation profiles have most of the differential rotation contrast confined to a low-latitude band between about \(\pm30^\circ\). The result is relatively strong radial shear distributed far more evenly throughout the CZ in the simulations than in the Sun. Equivalently, each simulated tachocline is not thin, but basically occupies the whole convective layer and is centered near mid-CZ, well above \(r_0\).

In order to define the location and width of the simulated tachoclines, we define: \[\begin{align} \label{eq:psi} \psi(r) &\equiv\frac{\Delta\Omega(r) - {\rm min}(\Delta\Omega)}{{\rm max}(\Delta\Omega) - {\rm min}(\Delta\Omega)} - \frac{1}{2}. \end{align}\tag{20}\] The shape function \(\psi(r)\) is normalized to vary between \(-1/2\) where \(\Delta\Omega\) obtains its minimum value (always in the RZ) and \(+1/2\) where \(\Delta\Omega\) obtains its maximum value (always in the CZ). We define a given tachocline’s centroid \(r_t\) and thickness \(\Gamma\) as the parameters in the function \((1/2)\tanh[2(r-r_t)/\Gamma]\) which is the best fit to \(\psi(r)\).³

Figure 5 shows the \(\psi\) profiles for the Sun and some medium- and strong-field cases, along with the corresponding best-fit tanh functions. As we would expect, the solar \(\psi\) (or best-fit tanh) profile has a centroid near \(r=r_0\) and a relatively narrow width. For each simulated tachocline, the distributed radial shear in the CZ both widens the \(\psi\) (or tanh) profile and pushes its centroid close the middle of the CZ.

Table 3 shows the tachocline parameters, as well as \(\Delta\Omega_{\rm{CZ}}\), \(\Delta\Omega_{\rm{RZ}}\), \(f\), and \(\Omega_{\rm{RZ}}\) for our simulations and the Sun. Clearly stronger fields reduce rotation contrast everywhere, but significantly more so in the RZ. Interestingly, the “tachocline contrast ratio" \(f\) is a non-monotonic function of regime and appears to be minimized (to a value of roughly 0.11) near \({\rm{Pr_m}}=4\). The nominal solar value is quite a bit higher: \(f_\odot\sim0.3\) from Table 3. However, this is mostly due to the uncertain tachocline width. The solar”\(\Delta\Omega_{\rm{RZ}}\)" thus contains substantial contrast from the lower half of the tachocline. Deeper in the RZ \((r/R_\odot\lesssim0.6)\) the helioseismic inversion gives \(\Delta\Omega\sim0.014\) or \(f_\odot\sim0.07\), which is probably closer to the value simulations should tend towards to be considered sufficiently solar-like.

In contrast to the Sun, all the simulated tachoclines have centroids well within the CZ and are roughly twice as thick. There is no clear scaling of the tachocline centroid with regime. However, the tachoclines in the strong-field regime are very slightly thinner than the tachoclines in the medium-field regime. Overall, it seems unlikely that further increasing the magnetic field strength will push the tachoclines to be more solar-like. Furthermore, there does not appear to be a solar-like “sweet spot" (say, in between the medium- and weak-field regimes), wherein the RZ rotates like a solid body, but strong rotation contrast is sustained in the CZ.

4.2 Torque Balance↩︎

In , we explicitly showed that the magnetic torque from the cycling, non-axisymmetric dynamo field was responsible for confining the tachocline in Case 4.00. This remains true for all the medium- and strong-field cases here. In the equilibrated state, the zonally and temporally averaged \(\phi\)-component of the momentum Equation 8 yields \[\label{eq:torque} \begin{align} &\underbrace {-\nabla\cdot[\overline{\rho}r \sin\theta \left\langle u^\prime_\phi \boldsymbol{u}_{\rm{pol}}^\prime\right\rangle_{\phi,t}]}_{\tau_{\rm{rs}}\text{ (Reynolds stress)}} -\underbrace{\overline{\rho}\left\langle \left\langle\boldsymbol{u}_{\rm{pol}}\right\rangle_{\phi}\cdot\nabla\mathcal{L} \right\rangle_{t} }_{\tau_{\rm{mc}}\text{ (meridional circulation)}}+\nonumber\\ &+\underbrace{{\rm{Ek}}\nabla\cdot\left[\overline{\rho}\overline{\nu}r^2\sin^2\theta\nabla\Omega\right]}_{\tau_{\rm{v}}\text{ (viscous)}} + \underbrace{\nabla\cdot \left[r\sin\theta \left\langle B_\phi^\prime\boldsymbol{B}_{\rm{pol}}^\prime\right\rangle_{\phi,t}\right]}_{\tau_{\rm{ms}}\text{ (Maxwell stress)}}\nonumber\\ &+ \underbrace{\nabla\cdot \left[r\sin\theta \left\langle \left\langle B_\phi\right\rangle_{\phi}\left\langle\boldsymbol{B}_{\rm{pol}}\right\rangle_{\phi} \right\rangle_{t}\right]}_{\tau_{\rm{mm}}\text{ (mean magnetic)}} = 0,\\ &\text{where}\ \ \ \ \ \underbrace{\mathcal{L}\equiv r\sin\theta(r\sin\theta + \left\langle u_\phi\right\rangle_{\phi}) }_{\text{angular momentum density}}, \end{align}\tag{21}\] where for a vector field \(\boldsymbol{A}\) we define its poloidal component \(\boldsymbol{A}_{\rm pol}\equiv A_r\hat{\boldsymbol{e}}_r+ A_\theta\hat{\boldsymbol{e}}_\theta\). The torques can be attributed to the physical processes labeled underneath each term (see [56], [57]).

Figure 6 shows the full steady-state torque balance (in the CZ and RZ separately) for Case 1.06 (the strongest weak-field case) and the strong-field Case 4.00. The CZ of Case 1.06 [Figure 6 (a)] is effectively hydrodynamic in its torque balance. It represents the “standard" by which current global models (e.g., [57]–[59]) maintain a solar-like differential rotation with fast equator and slow pole. That is, the rotational influence on the convection lead to Taylor columns with correlations in the components of \(\boldsymbol{u}^\prime\) (i.e., Reynolds stresses), which transport angular momentum away from the rotation axis and produce mostly positive torques at low latitudes (\(\lesssim45^\circ\); further from the rotation axis) and negative torques at high latitudes (\(\gtrsim45^\circ\); closer to the rotation axis). Meridional circulation also plays a role (a complicated one due to the presence of multiple circulation cells), especially at low latitudes. Viscosity always tries to eliminate gradients in \(\Omega\), in this case the latitudinal gradients, by spinning the equator down and the polar regions up. Note that this downward viscous spread of differential rotation is distinct from spread along poloidal field lines according to Ferraro’s law (e.g., [60], [61]) and in general, the RZ’s isorotation contours in the weak-field cases do not fall along poloidal field lines. The RZ of Case 1.06 [Figure 6 (b)] has a torque balance that is roughly an imprint of the balance in the CZ, but is overall much weaker and concentrated at high latitudes (and the magnetic torques are negligible). The torque balance in the CZ of Case 4.00 [Figure 6 (c)] still has a positive Reynolds-stress torque, but this positive torque is confined to significantly lower latitudes (\(\lesssim15^\circ\)). Consequently, most differential rotation is confined to a narrow prograde jet at the equator. The Maxwell-stress torque opposes the Reynolds-stress torque and effectively acts as an additional source of viscous torque. The meridional-circulation torque is significantly altered from its weak-field counterpart as well. Evidently, strong-field magnetism not only provides an additional torque, but also changes the structure of the convection and circulation so as to alter the hydrodynamic torques from their weak-field forms.

Finally, the torque balance in the RZ of Case 4.00 [Figure 6 (d)] was studied in and clearly this is the balance responsible for tachocline confinement. The profile of viscous torque has changed sign compared to the torque profiles in the other panels: it is now positive at low latitudes (\(\lesssim15^\circ\)) and negative at high latitudes (\(\gtrsim15^\circ\)), thus trying to imprint the equatorial jet and weak high-latitude retrograde differential rotation downward.⁴ The viscous torque is countered by the magnetic torque, which must come from the large-scale, non-axisymmetric (\(m=1,2\)) field components shown in Figure 2 (that this is true, at least for Case 4.00, was shown explicitly in ).

All our weak- and strong-field cases have torque balances like those in Figure 6. The medium-field cases have balances essentially similar to the strong-field cases, but the torques become more complicated due to the intermittent changes in field strength and differential rotation that was noted in connection with Figure 1. Regardless, the answer to how our simulated tachoclines are confined reduces to explaining the maintenance of large-scale, non-axisymmetric magnetism in the RZ. The following sections show how this maintenance can be understood in terms of the cycling dynamo and skin effect.

5 Cycling Behavior↩︎

5.1 Dynamo Cycles in the Weak- and Strong-Field Regimes↩︎

Figure 7 (left-hand panels) shows time-latitude diagrams of \(\left\langle B_\phi\right\rangle_{\phi}\) for the weak-field Case 1.00 and \(\text{real}(B_{\phi,1})\) for the strong-field Case 4.00 at two depths, one in the CZ and one in the RZ. Both cases cycle, although the polarity reversals in the weak-field case occur significantly more regularly than the reversals in the strong-field case. In each case, the cycle “imprints" from the base of the CZ onto the RZ with a phase lag (i.e., for every reversal in the CZ, there is a corresponding reversal in the RZ some time later). There is also significantly more rapid variation in the large-scale field in the CZ (seen as graininess in the time-latitude plots) than in the RZ. This again suggests that the RZ acts as a low-pass filter, in time as well as in space.

To describe these cycles more precisely, we define the frequency components of each \(\boldsymbol{B}_m\): \[\begin{align} \boldsymbol{B}_{m\omega} &\equiv\left\langle\boldsymbol{B}_m e^{i\omega t}W(t)\right\rangle_{t} = \left\langle\boldsymbol{B}e^{-i(m\phi-\omega t)}W(t)\right\rangle_{\phi,t},\label{eq:bmw} \end{align}\tag{22}\]

where \(W(t)\) is the Hanning window function and \(\omega\) is the discrete angular frequency. From the convention in the exponential (for nonzero \(m\) only), the \(\boldsymbol{B}_{m\omega}\) components with positive \(\omega/m\) move prograde in longitude and the components with negative \(\omega/m\) move retrograde.

We sample the spherical-slice magnetic-field data during the equilibrated state (\({t_{\rm{eq}}}\) to \({t_{\rm{max}}}\)). The sampling intervals are not uniform within a given simulation, but they are typically close to the mean interval \(\delta t\sim3\)–\(4 {P_{\rm{rot}}}\), with a typical standard deviation of \(\sigma_t\sim0.1\)–\(0.4{P_{\rm{rot}}}\) (see Table 4). We thus interpolate the non-uniform time series onto a uniform time series spaced by \(\delta t\) before computing the (windowed) discrete Fourier transform represented by Equation 22 .

Figure 7 (right-hand panels) shows the power in the large-scale toroidal field (\(|B_{\phi,0\omega}|^2\) for the weak-field case and \(|B_{\phi,1\omega}|^2\) for the strong-field case) corresponding to the time-latitude diagrams. The regularity of the weak-field cycle causes most of the power to be concentrated in the primary central frequency. By contrast, for the irregular strong-field cycle, there is a wide dispersion of power around a negative central frequency. This preference for negative frequencies suggests retrograde propagation of \(\boldsymbol{B}_1\), broadly consistent with transport by the negative background rotation rate in the RZ. Furthermore, the high-\(|\omega|\) “tail" in the strong-field case is significantly less pronounced in the RZ than in the CZ, again reinforcing the idea that the RZ acts as a low-pass filter in time.

Figure 7 (right-hand panels) shows that in each case, there is a central frequency (the “primary" cycle frequency \({\omega_{\rm{cyc}}}\)) and a dispersion (of width \(\sigma_\omega\)) in power about this central frequency. More precisely, for a given powerspectrum \(P(\omega)\), we define \({\omega_{\rm{cyc}}}\) as the median frequency associated with \(P(\omega)\) and \(\sigma_\omega\) as \(P(\omega)\)’s half-integral width: \[\begin{align} \sum_{\omega\leq{\omega_{\rm{cyc}}}}P(\omega)= \sum_{\omega={\omega_{\rm{cyc}}}-\sigma_\omega/2}^{{\omega_{\rm{cyc}}}+\sigma_\omega/2} P(\omega)\equiv\frac{1}{2} \sum_\omega P(\omega).\label{eq:omcyc} \end{align}\tag{23}\] The cycle period is \({P_{\rm{cyc}}}\equiv 2\pi/{\omega_{\rm{cyc}}}\) (since \({P_{\rm{rot}}}=2\pi\), note that \({P_{\rm{cyc}}}/{P_{\rm{rot}}}=1/{\omega_{\rm{cyc}}}\)). The quantity \(q\equiv{\omega_{\rm{cyc}}}/\sigma_\omega\) defines the regularity of the cycle, with higher \(q\) indicating a more regular cycle.

Table 4 contains values of \({\omega_{\rm{cyc}}}\), \({P_{\rm{cyc}}}\), \(\sigma_\omega\), and \(q\), along with the sampling parameters \(\delta t\), Nyquist frequency \({\omega_{\rm{nyq}}}\), and frequency resolution \(\delta\omega\). For the weak-field cases, we take \(P(\omega)= \left\langle|\boldsymbol{B}_{0\omega}|^2\right\rangle_{\rm{full}}\) (considering positive \(\omega\) only, since \(\boldsymbol{B}_0=\left\langle\boldsymbol{B}\right\rangle_{\phi}\) is real) and for the non-weak-field cases, we take \(P(\omega)= \left\langle|\boldsymbol{B}_{1\omega}|^2\right\rangle_{\rm{full}}\) (considering both positive and negative \(\omega\)). The weak-field solutions all have similar cycle periods (\({P_{\rm{cyc}}}\sim1400\)–\(2000{P_{\rm{rot}}}\)), with relatively high values of \(q\). This confirms the visual appearance of regular cycles in the weak-field cases (Figure 7).

The medium- and strong-field cases have more irregular cycles (with \(q\lesssim1\)) and the cycle period (with the exception of either one of cases 1.67 or 2.00) monotonically increases with increasing field strength. Since field strength increases with \({\rm{Pr_m}}\) and therefore with magnetic diffusion time \(P_\eta\) (see Table 5), this suggests that the cycle period for the non-weak cases is at least partly determined by the level of diffusion (i.e., \({P_{\rm{cyc}}}\) scales more or less monotonically with \(P_\eta\)).

6 Skin-Depth Interpretation↩︎

As mentioned at the conclusion of Section 4.2, explaining the presence (or not) of tachoclines in these simulations boils down to the maintenance of large-scale, non-axisymmetric (\(m=1,2\)) \(\boldsymbol{B}_{\rm{pol}}\) in the RZ.⁵ In , we showed that two effects were responsible for this maintenance: induction (possibly from inertial oscillations; see also [63]) and diffusion of CZ-produced field to roughly a (then ill-defined) skin-depth below the CZ. In this section, we precisely define the relevant skin effect and we show how the amplitude of \(\boldsymbol{B}_{\rm{pol}}\) in the RZ can be extremely well-predicted considering only diffusive skin effects.

As a first approximation, we assume that fluid motions produce no electromotive force (e.m.f.) below \(r_0\) (or a radius slightly below \(r_0\) for the weak-field cases; see Figure 8‘s caption). Then the evolution of \(\boldsymbol{B}_{\rm{pol}}\) in the RZ is governed by diffusion alone, with the upper boundary condition (at \(r=r_0\)) that \(\boldsymbol{B}_{\rm{pol}}\) matches what the CZ produces and the lower boundary condition (at \(r=r_{\rm{in}}\)) that the field decays with depth instead of grows. For axisymmetric weak-field dynamos, the regular polarity reversals provide an oscillating boundary condition at a single frequency. This is the classic form of Stokes’ problem of an oscillating boundary. In its solution, the field amplitude is contained in an envelope that decays exponentially downwards with a scale height (in this context, called the skin depth) that depends on the frequency of oscillation. This is the formalism expounded in the original fast magnetic confinement scenario of [21]. Note that in this axisymmetric case, the rotation rate of the frame in which the equations are solved does not matter.

However, for the non-axisymmetric medium- and strong-field dynamos, the choice of rotating frame does matter. Since advection in \(\phi\) of a non-axisymmetric \(\boldsymbol{B}_{\rm{pol}}\) constitutes an e.m.f., diffusion-only evolution is possible only if the RZ rotates approximately like a solid body. Then, to examine purely diffusive solutions, the induction equation must be written in the frame rotating at the solid-body rate \(\Omega_{\rm{RZ}}\) (see Table 3 for the simulated and solar values of \(\Omega_{\rm{RZ}}\)). Because the field at \(r=r_0\) is cycling with multiple frequencies (see the previous Section 5), this setup still corresponds to Stokes’ problem, but there is now a different skin-depth for each component \(\boldsymbol{B}_{{\rm pol},m\omega}\). Furthermore, since the equations must be solved in the frame of the RZ, the frequency determining the skin-depth is not \(\omega\), but the Doppler-shifted value \(\omega-m\Omega_{\rm{RZ}}\).⁶

Assuming that the spatial variation of \(\boldsymbol{B}_{\rm{pol}}\) is predominantly radial, Equation 10 leads to separate boundary-value problems for each \(\boldsymbol{B}_{{\rm pol},m\omega}\): \[\begin{align} \label{eq:indskin} -i(\omega-m\Omega_{\rm{RZ}})\boldsymbol{B}_{{\rm pol},m\omega} &\approx \frac{{\rm{Ek}}}{{\rm{Pr_m}}}\overline{\eta}(r)\frac{\partial^2\boldsymbol{B}_{{\rm pol},m\omega}}{\partial r^2} \end{align}\tag{24}\] for \(r\leq r_0\). Rapid variation in \(r\) allows us to neglect the terms in \(\nabla^2\) other than \((\partial/\partial r)^2\), sphericity terms, and the term from \(\nabla\overline{\eta}\). Note that Equation 24 is valid for all \(m\).

Because \(\overline{\eta}(r)\) varies with radius, we follow [64] and define \[\begin{align} \label{eq:reta} r_\eta \equiv r_{\rm{in}}+ \frac{\int_{r_{\rm{in}}}^{r}\overline{\eta}(r^\prime)^{-1/2}dr^\prime}{\int_{r_{\rm{in}}}^{r_0}\overline{\eta}(r^\prime)^{-1/2}dr^\prime}. \end{align}\tag{25}\] Note that \(r_\eta\) is a monotonically increasing function of \(r\) and is equal to \(r\) at \(r=r_{\rm{in}}\) and \(r=r_0\).⁷ Again assuming rapid radial variation, Equation 24 becomes \[\label{eq:indskin2} \begin{align} -i(\omega-m\Omega_{\rm{RZ}})\boldsymbol{B}_{{\rm pol},m\omega} \approx \frac{{\rm{Ek}}}{{\rm{Pr_m}}}\overline{\eta}_{\rm{const}}\frac{\partial^2\boldsymbol{B}_{{\rm pol},m\omega}}{\partial r_\eta^2},\\ \text{where}\ \ \ \ \ \overline{\eta}_{\rm{const}}\equiv\left[ \frac{r_0-r_{\rm{in}}}{\int_{r_{\rm{in}}}^{r_0}\overline{\eta}(r^\prime)^{-1/2}dr^\prime} \right]^2 \end{align}\tag{26}\] is an intermediate value of \(\overline{\eta}(r)\) in the RZ. For our chosen reference state, \(\overline{\eta}_{\rm{const}}=0.292\) and \(\overline{\eta}(r)\) achieves this value at \(r/R_\odot=0.599\).

Equation 26 is of course Stokes’ problem again and its exact solution yields \[\tag{27} \begin{align} \left\langle|\boldsymbol{B}_{{\rm pol},m\omega}|^2\right\rangle_{\rm{sph}}(r) = &\left\langle|\boldsymbol{B}_{{\rm pol},m\omega}|^2\right\rangle_{\rm{sph}}(r_0)\times\nonumber\\ & \exp{\left[-2\left(\frac{r_0-r_\eta}{\delta_{m\omega}} \right)\right]},\\ \text{where}\ \ \ \ \ \delta_{m\omega} \equiv& \;\sqrt{\frac{2{\rm{Ek}}\overline{\eta}_{\rm{const}}}{{\rm{Pr_m}}|\omega-m\Omega_{\rm{RZ}}|} }\tag{28} \end{align}\] is the \(m\)- and \(\omega\)-dependent skin-depth.

The “skin-predicted" amplitude of large-scale \(\left\langle|\boldsymbol{B}_{\rm{pol}}|^2\right\rangle_{ {\rm sph}, t}\) is then found by summing Equation 27 over all \(\omega\) and low \(m\). We choose \(m=0\) for the weak-field cases and \(m\in\{0,1,2\}\) for the medium- and strong-field cases. Figure 8 shows large-scale \(\left\langle|\boldsymbol{B}_{\rm{pol}}|^2\right\rangle_{ {\rm sph}, t}\) (both the skin-predicted and actually-realized values) for a weak-, medium-, and strong-field case. Equation 27 does an extremely good job of predicting the field strength for the weak- and strong-field cases and a reasonable job for the medium-field case. Overall, it thus seems highly likely that the magnetization of the RZ is determined primarily by the dynamo cycle of the CZ imprinting diffusively downward.

In , the strong \(\boldsymbol{B}_{\rm{pol}}\) in the RZ of Case 4.00 was attributed partially to deep dynamo action. For all the magnetic cases considered here, we have verified that the deep dynamo is still present, that is, the production of \(|\boldsymbol{B}_{\rm{pol}}|^2\) by diffusion (\(D_{\rm pol}\)) is negative in the RZ, while the production by e.m.f. (\(I_{\rm pol}\)) is positive.⁸ It was emphasized in that this implies (by definition; e.g., [65], p. 146) the presence of dynamo action deep in the RZ, and we argued in that this deep dynamo (possibly driven by Rossby waves) may have been responsible for tachocline confinement in Case 4.00. However, the results of this section indicate that the strength of \(\boldsymbol{B}_{\rm{pol}}\) in the medium- and strong-field cases can be almost fully accounted for by diffusive skin effects. It thus seems likely we would have tachocline confinement (in the simulations considered here) regardless of whether there was a deep dynamo or not. How the deep dynamo is driven—and whether it can confine the tachocline in the absence of large diffusion—remains an intriguing open question.

7 Polarity Reversals for Non-Axisymmetric Magnetic Fields↩︎

Polarity reversals in non-axisymmetric magnetic fields [e.g., Figures 7 (e,g)] can be accomplished in two distinct ways. For definiteness, consider \(B_{\phi,1}\) (i.e., the colatitudinal field associated with a single partial-wreath pair). At a given radius and latitude, we have \[\begin{align} \label{eq:tworeversals} W(t)B_{\phi,1}(t) = \sum_\omega B_{\phi,1\omega}e^{-i\omega t} = A(t)e^{i\varphi(t)} \end{align}\tag{29}\] The first equality comes directly from Equation 22 and the second equality is simply the mathematical statement that any complex number can be written as an amplitude [here \(A(t)\)] and a complex phase [here \(e^{i\varphi(t)}\)]. Note that the presence or not of the window function \(W(t)\) is immaterial to the following arguments.

The first type of non-axisymmetric polarity reversal is due to modulation of the amplitude \(A(t)\). These reversals contain cycle minima [for which \(A(t)=0\)] and are analagous to the reversals of full-wreath (i.e., axisymmetric) polarity in the weak-field cases [e.g., Figures 7 (a,c)], or equivalently to what we believe happens to the solar interior magnetic field to cause the observed butterfly diagram (e.g., [66]). The second type of non-axisymmetric polarity reversal is due to changes in the phase \(\varphi(t)\), which simply occur from advection of the whole structure in longitude (there are no cycle minima in this case). Equation 29 shows that in general, there is no straightforward way to separate which frequency components \(B_{\phi,1\omega}\) are due to each type of reversal. Indeed, in [32] (see Figures 11 and 12 from that paper), we showed that both types of reversal occur simultaneously in the CZ-only partial-wreath cycles, with the frequency of amplitude modulation similar to that of longitudinal advection.

Postponing for now the important investigation of how amplitude modulation occurs (it must be caused by non-axisymmetric dynamo processes; e.g., [67]–[70]), we discuss in this section which frequencies are consistent with longitudinal advection. Figure 9 shows Case 4.00’s poloidal powerspectra as functions of latitude and frequency for both \(m=1\) and \(m=2\) at the CZ–RZ interface \(r=r_0\). For \(m=1\) (panel a), the shape of the powerspectrum is nearly latitude-independent, with roughly constant values of the latitudinally dependent cycle frequency \({\omega_{\rm{cyc}}}(\theta)\) and dispersion \(\sigma_\omega(\theta)\). The central frequency \({\omega_{\rm{cyc}}}(\theta)\) overlaps with \(m\Omega(r_0,\theta)\) at low latitudes (about \(15^\circ\) north and south), which correspond to the retrograde (\(\Omega<0\)) jet of Figure 3 (c). For \(m=2\) (panel b), \({\omega_{\rm{cyc}}}(\theta)=m\Omega(r_0,\theta)\) at a slightly lower latitude (about \(10^\circ\) north and south) and the dispersion of power \(\sigma_\omega(\theta)\) is “stretched" by roughly a factor of two compared to the dispersion for \(m=1\). This stretching factor is consistent with the advective rate being proportional to \(m\)-value.

One valid interpretation of Figure 9 is that the large-scale (\(m=1,2\)) structures move as a cohesive structure (i.e., with rotation rates more or less independent of latitude), but at a time-varying rotation rate, which corresponds to a range of frequencies \({\omega_{\rm{cyc}}}\pm\sigma_\omega/2\) centered about the advective rate \(m\Omega(r_0,\theta)\) at low latitudes. Since these low latitudes are also the location of the retrograde jet, we may interpret the partial wreaths as being “anchored" to the jet. Another valid interpretation is that the partial wreaths have an intrinsic rate of rotation caused by the dynamo mechanism, which is apparently independent of latitude. This dynamo-intrinsic rotation rate then determines the rotation rate of the retrograde jet via the magnetic torques (since the magnetic torque should force the fluid to move with \(\boldsymbol{B}_{\rm{pol}}\)).

For the skin depth [Equation 28 ], it does not matter which physical process produces a given value of \(\omega\). As long as \(\omega\) is different from \(m\Omega_{\rm{RZ}}\) (see the orange lines in Figure 9), the amplitude of \(\boldsymbol{B}_{\rm{pol}}\) should decay downward with a finite skin-depth. Considering the first interpretation of Figure 9, we thus argue for the relevance of an important new type of skin effect that could operate in stars with both rigidly-rotating RZs and non-axisymmetric magnetic fields. This skin effect would arise from non-axisymmetric field at the CZ–RZ interface being advected by a background rotation rate that is different from the rotation rate of the RZ.

One particularly interesting consideration are the latitudes at the base of the CZ that co-rotate with the RZ. At those latitudes, \(\omega-m\Omega_{\rm{RZ}}=0\) and the skin depth in Equation 28 becomes infinite. What this really means is that any frozen-in non-axisymmetric field appears completely stationary to the RZ and spreads downward indefinitely on a diffusive time-scale. We explore this idea in the solar context in the following section.

8 Discussion: Non-Axisymmetric Dynamo Confinement of the Solar Tachocline↩︎

These results suggest that the fast magnetic confinement scenario—which was originally proposed, in 1D only, for axisymmetric \(\boldsymbol{B}_{\rm{pol}}\) cycling at a single frequency [21], [22], [71], [72]—should be expanded (into 3D) to include both non-axisymmetric \(\boldsymbol{B}_{\rm{pol}}\) and a spread in cycle frequencies. Whether this more general scenario is actually capable of confining the solar tachocline depends on several major differences between simulations and the Sun, which we briefly discuss here. Note that for this section (which is concerned with a real astrophysical object, namely, the Sun), we regard all physical quantities as dimensional.

Prior work on the fast magnetic confinement scenario has always assumed a turbulently enhanced magnetic diffusivity. For example, [22] (see Figure 5 from that paper) nicely show that for dynamo poloidal field strengths of \(\sim\) \(10^3\) G (and a cycle period of \(\sim\)​22 yr), the magnetic diffusivity must be larger than its molecular value by a factor of at least \(10^5\)–\(10^6\). However, as discussed in Section 1, how much turbulent enhancement of the viscosity occurs in the hydrodynamic scenario (and if the enhancement is primarily horizontal or vertical) is a subject of ongoing research with no firm conclusions at present. For magnetic diffusive enhancement, even less is known. Thus for simplicity, we assume here that the magnetic diffusion is not turbulently enhanced. We also leave aside for now ’s proposition that deep dynamo action may generate significant \(\boldsymbol{B}_{\rm{pol}}\).

8.1 Diffusive Equilibration in Simulations↩︎

In the Sun, all diffusive time-scales are significantly greater than the current solar age (\(t_\odot=4.6\) Gyr; see Table 6). By contrast, simulations that seek to address the tachocline confinement problem do so by evolving the MHD equations over significant fractions of the relevant diffusion times. If a statistically steady state is achieved, it thus likely contains significant diffusive effects in the dynamical balances.

We believe this may be one of the main reasons our tachocline cases have most of the differential rotation confined to a narrow equatorial jet near the outer boundary. The viscous and magnetic diffusion time-scales are similar (we have order-unity \({\rm{Pr_m}}\) values) and in most cases we run for several of each time-scale. The steady state thus necessarily has similar magnitudes for the viscous and magnetic torques in the CZ and RZ (compare the left-hand and right-hand columns of Figure 6). This means that any magnetic torque strong enough to prevent viscous tachocline spread is also strong enough to eliminate much of the differential rotation in the CZ.

8.2 Viscous versus Radiative Spread: General Torque Balance↩︎

Even barring the open question of whether circulation burrowing is hyperdiffusive in the Sun, it seems likely that radiative spread dominates viscous spread. This dominance is expressed via the “\(\sigma\)-parameter" (e.g., [11], [12], [73]–[75]): \[\begin{align} \label{eq:sigma} \sigma\equiv\sqrt{\frac{{P_{\rm{ES}}}}{{P_{\nu}}_{,\rm RZ}}}=\frac{\sqrt{{\rm{Pr}}\rm Bu}}{2} \end{align}\tag{30}\] For the Sun, \(\sigma_\odot=0.17\ll1\) (see Table 6). Since the Reynolds number in the solar CZ is extremely high, the viscous torque should drop out of the torque balance in the CZ as well. Global simulations seem to indicate that large-scale magnetic field (when strong enough) significantly reduces the differential rotation in the CZ (e.g., [32], [33], [45]–[47], [76]–[79]). For a fast magnetic confinement scenario to work (i.e., a scenario in which the dynamo-produced magnetic field diffusively penetrates into the upper RZ), we thus might require that the total magnetic torque (\(\tau_{\rm{mag}}\equiv\tau_{\rm{ms}}+ \tau_{\rm{mm}}\)) be both large enough in the RZ to counter radiative spread and small enough to drop out of the torque balance in the CZ. In that case, Equation 21 (its dimensional counterpart) becomes \[\begin{align} \label{eq:torquegeneral} 0=\begin{cases} &\underbrace{-\frac{4\Omega_\odot^2}{\overline{N^2}}r_0^2\overline{\rho}\;\overline{\kappa}\frac{\partial^4\left\langle\mathcal{L}\right\rangle_{t}}{\partial r^4}}_{\tau_{\rm{rad}}\text{ (radiative spread)}} + \tau_{\rm{mag}}\ \ \ \ \ \text{in the RZ}\\ &\tau_{\rm{rs}}+ \tau_{\rm{mc}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \;\;\; \text{in the CZ,} \end{cases} \end{align}\tag{31}\] where the form of \(\tau_{\rm{rad}}\) is derived in [6] [their Equation (4.9)] and we have assumed a thin tachocline (so that we retain only highest derivatives in \(r\)). We estimate \(\partial^4\left\langle\mathcal{L}\right\rangle_{t}/\partial r^4\sim (r_0/\sqrt{2})^2\Delta\Omega_{\rm{CZ}}/\Gamma_\odot^4\). If we take \(\Gamma_\odot=0.05R_\odot\) and take the Model S values in Table 6 (averaged over the upper solar RZ) for \(\overline{\rho}\), \(\overline{\kappa}\), and \(\overline{N^2}\), and take \(\Delta\Omega_{\rm{CZ}}=0.20\Omega_\odot\) from Table 3, we find \[\begin{align} \label{eq:tauradest} \tau_{\rm{mag}}\sim \tau_{\rm{rad}}\sim 0.84\;\rm{dyn}\ {\rm{cm}}^{-2}\ \ \ \ \ \text{in the RZ}. \end{align}\tag{32}\]

Meanwhile in the CZ, the Reynolds-stress torque has not been measured helioseismically (although could be in the future via ring analysis; e.g., [80]–[82]). Nonetheless, the meridional flow’s amplitude \(|\boldsymbol{u}_{\rm{pol}}|\sim10\;{\rm m\;s^{-1}}\) is fairly well known, at least in the upper half of the CZ (e.g., [83]–[85]) and so we estimate \(\tau_{\rm{mc}}\sim (3/2\pi)\tilde{\rho}|\boldsymbol{u}_{\rm{pol}}|R_\odot\Delta\Omega_{\rm{CZ}}\), or \[\begin{align} \label{eq:taumcest} \tau_{\rm{rs}}\sim \tau_{\rm{mc}}\sim 1.2\times 10^{6}\;\rm{dyn}\ {\rm{cm}}^{-2}\ \ \ \ \ \text{in the CZ}. \end{align}\tag{33}\] Equations 32 and 33 suggest that a diffusively coupled solar CZ and RZ (in which the magnitude of \(\tau_{\rm{mag}}\) is similar in both zones) can support the fast magnetic confinement scenario, i.e., \(\tau_{\rm{rad}}\sim\tau_{\rm{mag}}\ll \tau_{\rm{rs}}\sim\tau_{\rm{mc}}\). We can further express \(\tau_{\rm{mag}}\) from the large-scale non-axisymmetric field in terms of field strength: \(\tau_{\rm{mag}}\sim [1/(2\sqrt{2}\pi^2)]|B_{\phi}|^2\). Here, we have (crudely) assumed that \(|\boldsymbol{B}_{{\rm pol}}|\sim|B_{\phi}|\), that \(r\sin\theta\sim r/\sqrt{2}\), and that the typical length-scale for large-scale field variation is \(\sim\) \(\pi r/2\). Equations 32 and 33 then yield \[\begin{align} \label{eq:fieldstrength} 4.8 {\rm{G}}\;\lesssim |B_{\phi}| \ll 5800\;{\rm{G}}. \end{align}\tag{34}\] Equation 34 states that if the fast magnetic confinement scenario operates in the Sun, we expect the zonal field strength to be significantly less than \(5800\) G in the CZ (so as not to disturb the torque balance there) and to diffusively decay to a lower bound of at least \(4.8\) G in the tachocline region (to counter radiative spread). The value of the lower bound depends strongly on the actual value of \(\Gamma_\odot\) and the value of the upper bound on the reliability of the simulations’ prediction that strong field quenches differential rotation.⁹

8.3 Small Skin Depths and Spread of a Permanent Dynamo Field↩︎

Equation 28 shows that, except for \(\omega=m\Omega_{\rm{RZ}}\), any oscillatory component of the solar dynamo has a very small skin depth and thus cannot significantly penetrate into the RZ. This is the reason why prior 1D models like [22] required an \(\overline{\eta}\) greatly enhanced from its molecular value. Explicitly, we rewrite Equation 28 in dimensional form as \[\begin{align} \delta_{m\omega}=\left( \frac{2\left\langle\eta\right\rangle_{\rm{RZ}}}{|\omega-m\Omega_{\rm{RZ}}|}\right)^{1/2}= (0.027R_\odot) {P_{\rm{cyc}}}^{1/2}, \end{align}\] where here \({P_{\rm{cyc}}}\equiv 2\pi/|\omega-m\Omega_{\rm{RZ}}|\) and is measured in Gyr. If we require diffusive spread over (say) \(\Gamma_\odot=0.05R_\odot\), we need \({P_{\rm{cyc}}}\sim1.4\) Gyr. With the solar age at \(t_\odot=4.6\) Gyr, such a high \({P_{\rm{cyc}}}\) cannot unambiguously constitute a “cycle" and instead better corresponds to the permanent component of \(\boldsymbol{B}_{\rm{pol}}\) (as viewed in the frame rotating with the RZ), here denoted by \(\boldsymbol{B}_{\rm pol,perm}\).¹⁰ There are few, if any, constraints on the solar \(|\boldsymbol{B}_{\rm pol,perm}|\), only that it is significantly less than \(|\boldsymbol{B}_{\rm{pol}}|\) (e.g., [86]). It is not obvious, however, how much less than \(|\boldsymbol{B}_{\rm{pol}}|\) it really is, and thus whether we can rule out a dynamo confinement scenario entirely if \(\overline{\eta}\) is not turbulently enhanced.

For example, even if the solar dynamo were purely axisymmetric and perfectly cyclic with a period of 22 yr, we would expect at most \(N_{\rm cyc} = (4.6\;{\rm Gyr})/(22\;{\rm yr})=2.1\times 10^{8}\) cycles since the dynamo turned on. If we assume there have always been random modulations of the cycle amplitude (as are observed throughout recorded history), then we estimate \(|\left\langle\boldsymbol{B}_{\rm pol,perm}\right\rangle_{\phi}|=|\left\langle\boldsymbol{B}_{\rm{pol}}\right\rangle_{\phi}|/\sqrt{N_{\rm{cyc}}} = 6.9\times 10^{-5}|\left\langle\boldsymbol{B}_{\rm{pol}}\right\rangle_{\phi}|\). Given Equation 34 , this reduced field strength would only be a factor of \(\sim\)​10 too small to confine the tachocline.¹¹

As noted at the end of the last Section 7, any non-axisymmetric \(\boldsymbol{B}_{\rm{pol}}\) that co-rotates with the RZ is effectively non-oscillatory and thus contributes to \(\boldsymbol{B}_{\rm pol,perm}\). If the fast magnetic confinement scenario is generalized to include non-axisymmetric fields, it thus seems possible that \(\boldsymbol{B}_{\rm pol,perm}\) (including all \(m\)’s) could be significantly larger and more topologically complex than prior estimates like [64]. It is also worth noting that many prior simulations using a variety of codes (e.g., [23], [29], [47], [87]) have all suggested that large-scale magnetic field accumulates preferentially in the tachocline region. Furthermore, the presence of a tachocline was argued to significantly stabilize the large-scale magnetic fields, sometimes lengthening the dynamo cycle period or even producing time-steady dynamos.

Active longitudes (preferential solar longitudes at which sunspots emerge; e.g., [88]–[91]) are particularly striking as a possible contributor to non-axisymmetric \(\boldsymbol{B}_{\rm pol,perm}\). Although it would be a major leap to claim that active longitudes imply a permanent interior partial-wreath structure co-rotating with the RZ (authors have done so nonetheless; e.g., [92]), it is intriguing that: (1) they often come in opposite-polarity pairs separated in longitude by \(180^\circ\) (e.g., [93], [94]) and (2) they seem to persist, in a properly chosen rotating frame (or in a frame with time-dependent rotation rate), for long time-scales: 20 years [95] or even \(\sim\)​100 years [96].

Whatever the source of \(\boldsymbol{B}_{\rm pol,perm}\), it should penetrate into the RZ much deeper than any skin-depth. Considering the Rayleigh problem (i.e., Stokes’ first problem, of a boundary plate suddenly jerked from rest), we estimate (for \(r\leq r_0\)): \[\label{eq:bpolperm} \begin{align} |\boldsymbol{B}_{\rm pol,perm}|(r) &= |\boldsymbol{B}_{\rm pol,perm}|(r_0) {\rm erfc}\left( \frac{r_0-r}{\delta_{\rm perm}}\right),\\ \text{where}\ \ \ \ \ \delta_{\rm perm} &= \sqrt{4 \left\langle\eta\right\rangle_{\rm{RZ}}t_\odot} = 0.21R_\odot. \end{align}\tag{35}\] For \(r_0-r=\Gamma_\odot=0.05R_\odot\), we find \({\rm erfc}(0.05/0.21) = 0.73\), i.e., there should be only a \(\sim\)​27% reduction in \(|\boldsymbol{B}_{\rm pol,perm}|\) over the depth of the tachocline.

8.4 Conclusion↩︎

In summary, we have performed a suite of dynamo simulations in which tachocline confinement is achieved if the large-scale non-axisymmetric fields (partial wreaths) are strong enough. These partial-wreath structures cycle with frequencies consistent with advection by a low-latitude retrograde jet. The structures thus appear to cycle from the perspective of the rigidly-rotating RZ and penetrate diffusively downward, with the amplitude of the confining \(\boldsymbol{B}_{\rm{pol}}\) very well-predicted by the skin-depth Equation 27 .

As a whole, the simulations presented here effectively achieve a fast magnetic confinement scenario [21], which is now generalized to include non-axisymmetric fields and a spread in cycle frequencies. Our work thus offers a significantly wider range of applicability to the fastmagnetic confinement scenario. To further constrain if such a scenario is consistent with observations, we might recommend that future work explore in greater detail the processes giving rise to non-axisymmetric magnetic field (such as active longitudes), and determine observationally how fast active-longitude pairs rotate with respect to the RZ.

In this discussion section, we have argued that if the magnetic diffusivity is not enhanced, then only an effectively permanent component of the solar dynamo can play a role in tachocline confinement. This component can include both the axisymmetric dynamo field (averaged in time since the birth of the Sun) and, possibly more importantly, any non-axisymmetric field that co-rotates with the RZ. In addition to the fast magnetic confinement scenario, we thus might also recommend exploring a more general (possibly slow) “dynamo confinement scenario." This would be similar to the model of [20], but with the permanent dynamo field taking the place of the primordial field. lt would differ from [20] mainly in that no primordial field would need to be confined to the RZ.

Finally, in order to make further progress on the numerical side, future simulations need to accomplish several computationally challenging tasks. First, they need to achieve equilibrium that is not diffusively controlled.¹² Second, they must be run in the \(\sigma\lesssim1\) regime [see Equation 30 ]; only then can we assess whether a dynamo confinement scenario can operate in the solar regime of little viscous torque. Finally, simulations must be run with small skin-depths \(\delta_{m\omega}\) (i.e., low \({\rm{Ek}}\) or high \({\rm{Pr_m}}\)). Skin-depths as small as in the Sun would not be possible, but we may at least achieve \(\delta_{m\omega}<\Gamma\), which would help confirm if the permanent dynamo field could penetrate deeply enough [possibly according to Equations 35 ] to confine the tachocline.

9 Background State↩︎

Note that in this section, we discuss both the non-dimensional and dimensional versions of various quantities. To explicitly distinguish, we denote the dimensional version of a quantity with a “dim" subscript (quantities like \(c_{\rm{p}}\) and \(\tilde{\nu}\), which are always dimensional, do not require a subscript).

In terms of the dimensional background state, the perfect-gas law is \[\begin{align} \label{eq:idgasdim} \overline{P}_{\rm{dim}}= \left[\frac{(\gamma-1)c_{\rm{p}}}{\gamma}\right]\overline{\rho}_{\rm{dim}}\overline{T}_{\rm{dim}}, \end{align}\tag{36}\] hydrostatic balance is \[\begin{align} \label{eq:hydrdim} \frac{d\overline{P}_{\rm{dim}}}{dr_{\rm{dim}}}=-\overline{\rho}_{\rm{dim}}\overline{g}_{\rm{dim}}, \end{align}\tag{37}\] and the first law of thermodynamics is \[\begin{align} \label{eq:firstlawdim} \frac{1}{c_{\rm{p}}}\left(\frac{dS_{\rm{dim}}}{dr_{\rm{dim}}}\right)=\frac{1}{\gamma}\frac{d\ln\overline{T}_{\rm{dim}}}{dr_{\rm{dim}}}-\left(\frac{\gamma-1}{\gamma}\right)\frac{d\ln\overline{\rho}_{\rm{dim}}}{dr_{\rm{dim}}}. \end{align}\tag{38}\]

After non-dimensionalizing, Equations 36 –38 take the form \[\begin{align} \label{eq:idgas} \overline{P}&= \overline{\rho}\overline{T}, \end{align}\tag{39}\] \[\begin{align} \label{eq:hydr} \frac{d\overline{P}}{dr}=-{\rm{Di}}\left(\frac{\gamma}{\gamma-1}\right)\overline{\rho}\,\overline{g}, \end{align}\tag{40}\] and \[\begin{align} \label{eq:firstlaw} \frac{d\overline{S}}{dr}=\frac{1}{\gamma}\frac{d\ln\overline{T}}{dr}-\left(\frac{\gamma-1}{\gamma}\right)\frac{d\ln\overline{\rho}}{dr}, \end{align}\tag{41}\] where \(\overline{S}\equiv\overline{S}_{\rm{dim}}/c_{\rm{p}}\) and we recall that \({\rm{Di}}\equiv\tilde{g}H/(c_{\rm{p}}\tilde{T})\). We combine Equations 39 –41 to yield \[\begin{align} \frac{d\overline{T}}{dr}-\left(\frac{d\overline{S}}{dr}\right)\overline{T}&= -{\rm{Di}}\overline{g}, \end{align}\] which has the exact solution [after choosing, without loss of generality, \(\overline{S}(r_0)=0\)], \[\begin{align} \label{eq:tmphat} \overline{T}&= e^{\overline{S}}\left[\overline{T}(r_0) - {\rm{Di}}\int_{r_0}^r \overline{g}(x) e^{-\overline{S}(x)}dx\right]. \end{align}\tag{42}\] We then eliminate \(\overline{P}\) from Equations 39 and 40 to yield \[\begin{align} \label{eq:rhohat} \overline{\rho}&= \overline{\rho}(r_0) \exp{\left[-\left(\frac{\gamma}{\gamma-1}\right)\overline{S}\right]}\left[\frac{\overline{T}(r)}{\overline{T}(r_0)}\right]^{1/(\gamma-1)}. \end{align}\tag{43}\] There are three equations relating \(\overline{\rho}(r_0)\), \(\overline{T}(r_0)\), \({\rm{Di}}\), \(\gamma\), \(\beta\), and \(N_\rho\): two from our choice of non-dimensionalization—\((4\pi/V_{\rm{CZ}})\int_{r_0}^{r_{\rm{out}}}\overline{\rho}(r)r^2dr=1\) and \((4\pi/V_{\rm{CZ}})\int_{r_0}^{r_{\rm{out}}}\overline{T}(r)r^2dr=1\), where \(V_{\rm{CZ}}\equiv(4\pi/3)(r_{\rm{out}}^3-r_0^3)=(4\pi/3)[(1-\beta^3)/(1-\beta)^3]\)—and one from the definition \(N_\rho\equiv\ln[\overline{\rho}(r_0)/\overline{\rho}(r_{\rm{out}})]\). Thus, \(\overline{\rho}(r_0)\), \(\overline{T}(r_0)\), and \({\rm{Di}}\) may be regarded as functions of \(\gamma\), \(\beta\), and \(N_\rho\). For the values given in Table 1 (and our choices for \(\overline{g}\) and \(d\overline{S}/dr\) given below), we explicitly find \(\overline{\rho}(r_0)=2.67\), \(\overline{T}(r_0)=2.04\), and \({\rm{Di}}=1.72\).

For \(\overline{g}(r)\propto1/r^2\) and the condition \((4\pi/V_{\rm{CZ}})\int_{r_0}^{r_{\rm{out}}}\overline{g}(r)r^2dr=1\), we require \[\begin{align} \overline{g}(r) = \left[ \frac{1-\beta^3}{3(1-\beta)^3} \right]\frac{1}{r^2}. \end{align}\] To model the transition from convective stability to instability at the base of the CZ, we choose \(d\overline{S}/dr\) to be zero in the CZ, a constant positive (near-unity) value in the RZ, and continuously matched in between: \[\frac{d\overline{S}}{dr} = \begin{cases} \Sigma & r\leq r_0 - \delta\\ \Sigma\bigg{\{}1 - \left(1 - \Big{(}\frac{r-r_0}{\delta}\Big{)}^2\right]^2\bigg{\}} & r_0 - \delta < r < r_0\\ 0 & r\geq r_0, \end{cases} \label{eq:dsdr}\tag{44}\] where \(\Sigma=0.453\) (note that \(\Sigma\) is not really a free parameter, since it can always be absorbed into the fluid control parameter \(\rm Bu\)). The choice of quartic matching ensures that the ultimate stability transition (determined by the total entropy gradient \(d\overline{S}/dr+d\langle S\rangle_{\rm{sph}}/dr\) in the equilibrated state) is never too far from \(r_0\). By contrast, for a tanh matching [\(d\overline{S}/dr=(\Sigma/2)(1-\tanh{[(r-r_0)/\delta]})\); e.g., [43]], the stability transition can occur significantly above \(r_0\). In our cases, it could occur as high up as \(r_0+5\delta\), since \(d\left\langle S\right\rangle_{\rm{sph}}/dr\) is generally \(10^4\)–\(10^5\) times smaller than \(d\overline{S}/dr\) and \((1/2)[1-\tanh{(5)}]\approx5\times 10^{-5}\).

With \(d\overline{S}/dr\) and \(\overline{g}\) chosen, we numerically integrate Equations 42 and 43 to find \(\overline{\rho}\) and \(\overline{T}\). This approach to defining the background state (also used by [43]) has the main advantage that hydrostatic balance is satisfied everywhere, even in the transition region. This stands in contrast to polytropic matching (e.g., [59], [76]).

Note that equation 44 also defines the buoyancy frequency through \[\begin{align} \frac{\overline{N^2}}{\overline{g}} = \frac{d\overline{S}/dr}{\langle\overline{g}d\overline{S}/dr\rangle_{\rm{RZ}}}, \end{align}\] where \(\langle\overline{g}d\overline{S}/dr\rangle_{\rm{RZ}}=0.597\). We choose \(\overline{Q}(r)\) to occupy primarily the CZ: \[\overline{Q}= \frac{c}{2}\left[1 + \tanh{\left(\frac{r-r_0}{\delta_{\rm{heat}}}\right)}\overline{\rho}\overline{T}\right], \label{nxjqsflv}\tag{45}\] where \(c=0.944\). This value of \(c\) is required because \(\overline{Q}=Q_{\rm{dim}}H/\widetilde{F_{\rm{nr}}}\). The definition \((\overline{F_{\rm nr}})_{\rm{dim}}\equiv(1/Hr^2)\int_r^{r_{\rm{out}}}Q_{\rm{dim}}(x)x^2dx\) then yields \(1/c = (2\pi/V_{\rm{CZ}})\int_{r_0}^{r_{\rm{out}}}(1/r^2)\int_r^{r_{\rm{out}}}f(x)x^2dxdr\), where \(f(x)\equiv 1+\tanh[(x-r_0)/\delta_{\rm{heat}}]\).

Table 5 gives some additional (derivative) input parameters that can be computed from the parameters of Table 1 and the form of the reference state just described.

We compare this reference state to the standard solar Model S [98] in Figure 10. Note that Model S profiles (denoted by an “S" subscript) are originally in dimensional form.

For the dimensional molecular diffusivities associated with Model S, we define \[\begin{align} \overline{\nu}_S &\equiv1.2\times 10^{-16}\;\frac{(\overline{T}_S/{\rm{K}})^{5/2}}{\overline{\rho}_S/({\rm{g}}\ {\rm{cm}}^{-3})}\;\rm{cm^2\ s^{-1}},\\ \overline{\kappa}_S &\equiv\frac{16\sigma_{\rm SB}\overline{T}_S^3}{3\overline{\chi}_S\overline{\rho}_S^2(c_{\rm{p}})_S},\\ \overline{\eta}_S &\equiv5.2\times 10^{11}\;\frac{\ln\Lambda}{(\overline{T}_S/{\rm{K}})^{3/2}}\;\rm{cm^2\ s^{-1}}, \end{align}\] where \(\sigma_{SB}\equiv5.67\times 10^{-5}\;\rm{erg\;cm^{-2}\;s^{-1}\;K^{-4}}\) is the Stefan-Boltzmann constant and \(\overline{\chi}_S\) is the opacity from Model S. The forms of the molecular viscosity \(\overline{\nu}_S\) and the radiative thermal diffusivity \(\kappa_S\) are given in (e.g.) [99] via [100]. The form of \(\overline{\eta}_S\) is given in (e.g.) [101]. The Coulomb logarithm \(\ln\Lambda\) is tabulated by (e.g.) [102], and we approximate \(\ln\Lambda\approx2.5 + r/r_0\) (see [64]).

To non-dimensionalize Model S, we take \((R_\odot)_{\rm{dim}}=6.96\times 10^{10}\) cm and set \((r_{\rm{in}})_{\rm{dim}}=0.491(R_\odot)_{\rm{dim}}\), \((r_0)_{\rm{dim}}=0.719(R_\odot)_{\rm{dim}}\), and \((r_{\rm{out}})_{\rm{dim}}=0.947(R_\odot)_{\rm{dim}}\) [and thus \(H=0.228(R_\odot)_{\rm{dim}}\)]. This choice means we compare to the bottom three density scale-heights of Model S’s CZ, i.e., \(\ln{\{ \overline{\rho}_S[(r_0)_{\rm{dim}}]/\overline{\rho}_S[(r_{\rm{out}})_{\rm{dim}}]\}}=3\). For a given reference-state quantity \(\psi\), we then define \(\left\langle\psi\right\rangle_{\rm{CZ}}\) (or \(\tilde{\psi}\)) as a volume average of \(\psi\) over \(((r_0)_{\rm{dim}}, (r_{\rm{out}})_{\rm{dim}})\) and \(\left\langle\psi\right\rangle_{\rm{RZ}}\) as a volume average of \(\psi\) over \([(r_{\rm{in}})_{\rm{dim}}, (r_0)_{\rm{dim}}]\). We then scale \(\overline{\rho}_S\), \(\overline{T}_S\), and \(\overline{g}_S\) by \(\tilde{\rho}_S\), \(\tilde{T}_S\), and \(\tilde{g}_S\) (respectively), \(\overline{Q}_S\) by \((\widetilde{F_{\rm nr}})_S/H\), and \((\overline{N^2})_S\) by \(\left\langle(N^2)_S\right\rangle_{\rm{RZ}}\). Then for the rest of this section, the Model S profiles denote their non-dimensional forms.

Figure 10 shows that Rayleigh’s non-dimensional reference state is fairly solar-like, being equivalent to the adiabatic polytrope from our prior work (e.g., [34], [57], [103], [104]). The biggest discrepancies occur near \(r=r_0\), where our reference state has relatively wide and smooth transitions in \(\overline{N^2}\) and \(\overline{Q}\) compared to the narrow and sharp transitions from Model S.

Note that the background diffusivities from the simulations’ reference state, \(\overline{\nu}\), \(\overline{\kappa}\), and \(\overline{\eta}\), are specified independently from the thermodynamic profiles. We choose all simulation diffusivities to increase with radius like \(\overline{\rho}^{-1/2}\) (and of course they are normalized to have a volume-average over the CZ of 1). Note that this choice does not in any sense correspond to the non-dimensional Model S profiles, \(\overline{\nu}_S\), \(\overline{\kappa}_S\), and \(\overline{\eta}_S\).

10 Dimensional Solar Analog↩︎

Here, we “re-dimensionalize" the models considered in the current paper to match the presentation of cases H and 4.00 in . A non-dimensional simulation can be re-dimensionalized by assuming dimensional values for quantities like \(H\) and \(\Omega_0\) and then computing the associated scales for the fluid variables (e.g., \([S]\) and \([\boldsymbol{u}]\)) and reference-state profiles (e.g., \(\tilde{\rho}\) and \(\tilde{T}\)), as described in Section 2. To list the input dimensional quantities in the conventional way, we define the luminosity \(L\equiv\tilde{Q}V_{\rm{dim}}\), where \(V_{\rm{dim}}\equiv(4\pi H^3/3)(r_{\rm{out}}^3-r_{\rm{in}}^3)\). \(L\) is typically the control parameter that sets \(\widetilde{F_{\rm{nr}}}\). We also define the stellar mass \(M\equiv[(1-\beta^3)/3(1-\beta)^3](H^2/G)\tilde{g}\), so \(\overline{g}_{\rm{dim}}=GM/H^2r^2\). \(M\) is typically the parameter that sets \(\tilde{g}\). The full set of dimensional input parameters is then: \(H\), \(\Omega_0\), \(L\), \(M\), \(c_{\rm{p}}\), \(\tilde{\rho}\), \(\tilde{T}\), \(\tilde{\nu}\), \(\tilde{\kappa}\), \(\tilde{\eta}\), and \(\langle N^2\rangle_{\rm{RZ}}\).

Some of the input dimensional parameters are obviously redundant, so there are infinitely many ways to re-dimensionalize. The only requirement is that the chosen dimensional values be consistent with the input non-dimensional numbers. Historically, we in the solar and stellar communities have chosen stellar-like dimensional values for as many parameters as possible except for the diffusivities, which are chosen to be unrealistically high. This choice is exemplified in Table 6, which contains the scaling employed in . Most of the chosen parameters are solar-like, except for the diffusivities. The rotation rate \(\Omega_0\) is chosen to be about three times higher than the solar Carrington value.

The inherent non-uniqueness associated with dimensional simulations is one of the main reasons we report only the non-dimensional versions of the simulations in this work. For example, comparing the simulated \(\boldsymbol{B}_{\rm{dim}}\) (measured in G) to an observed \(\boldsymbol{B}\) at the solar surface (also measured in G) is fundamentally ambiguous. For, we could have re-dimensionalized using the values in Table 6, but instead chosen \(\Omega_0\rightarrow \Omega_\odot\), \(L\rightarrow L_\odot/27\), \(\langle N^2\rangle_{\rm{RZ}}\rightarrow \langle N^2\rangle_{\rm{RZ}}/3\), \(\tilde{\nu}\rightarrow\tilde{\nu}/3\), \(\tilde{\kappa}\rightarrow\tilde{\kappa}/3\), and \(\tilde{\eta}\rightarrow\tilde{\eta}/3\). This would have yielded dynamically identical simulations, but with all values of \(\boldsymbol{B}_{\rm{dim}}\) three times smaller.

11 Output Non-Dimensional Numbers↩︎

We define the Reynolds (Re), Rossby (Ro), and magnetic Reynolds (\({\rm{Re_m}}\)) numbers, separately for the mean and fluctuating flows: \[\begin{align} {\rm{Re}}_{\rm{mean}} \equiv\frac{(\left\langle\boldsymbol{u}\right\rangle_{\phi})_{\rm{rms}}}{{\rm{Ek}}}, \ \ \ \ \ {\rm{Re}}_{\rm{fluc}} \equiv\frac{(\boldsymbol{u}^\prime)_{\rm{rms}}}{{\rm{Ek}}},\tag{46}\\ {\rm{Ro}}_{\rm{mean}}\equiv\frac{(\left\langle\boldsymbol{\omega}\right\rangle_{\phi})_{\rm{rms}}}{2}, \ \ \ \ \ {\rm{Ro}}_{\rm{fluc}} \equiv\frac{(\boldsymbol{\omega}^\prime)_{\rm{rms}}}{2},\tag{47}\\ {\rm{Re_m}}_{,\rm{mean}}\equiv{\rm{Re}}_{\rm{mean}}{\rm{Pr_m}}, \ \ \ \ \ {\rm{Re_m}}_{,\rm{fluc}} \equiv{\rm{Re}}_{\rm{fluc}}{\rm{Pr_m}}\tag{48}, \end{align}\] where \(\boldsymbol{\omega}\equiv\nabla\times\boldsymbol{u}\) is the vorticity and the mean in the rms is taken in volume (over the CZ or RZ) and in time over the equilibrated state. Table 7 contains the values of these numbers for each simulation considered in this work.

We thank the following (undoubtedly incomplete) list of individuals for helpful discussions of the solar tachocline: Pascale Garaud, Mark Miesch, Lydia Korre, Nicholas Featherstone, Gustavo Guerrero, Connor Bice, Sacha Brun, Antoine Strugarek, Catherine Blume, Matthew Browning, Steven Tobias, David Hughes, and Jørgen Christensen-Dalsgaard. This work was partly done in collaboration with the COFFIES DRIVE Science Center (NASA grant 80NSSC22M0162). L.I.M. and N.H.B. thank the Isaac Newton Institute for Mathematical Sciences (Cambridge, UK) for support and hospitality during the 2022 Programme “Frontiers in Dynamo Theory: from the Earth to the Sun and Stars", where part of the work on this paper was undertaken. L.I.M. was primarily supported during this work by a National Science Foundation Astronomy & Astrophysics Postdoctoral Fellowship under award AST-2202253, as well as a Future Investigators in NASA Earth and Space Sciences Technology (FINESST) award 80NSSC19K1428. The computations integral to this work were supported by NASA grant 80NSSC18K1127. This research was further supported by NASA grants 80NSSC18K1125, 80NSSC19K0267, and 80NSSC20K0193. Computational resources were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. Rayleigh is hosted and receives support from the Computational Infrastructure for Geodynamics (CIG), which is supported by the National Science Foundation awards NSF-0949446, NSF-1550901, and NSF-2149126. The input files, final checkpoints, and some data-analysis products (averaged data and basic plots) for all simulations are publicly accessible via Zenodo [105].

References↩︎