###### Abstract

We characterize and analyze rotational torsional oscillations developing in a large-eddy magnetohydrodynamical simulation of solar convection (Ghizaru, Charbonneau, and Smolarkiewicz, Astrophys. J. Lett. 715, L133 (2010); Racine et al., Astrophys. J. 735, 46 (2011)) producing an axisymmetric large-scale magnetic field undergoing periodic polarity reversals. Motivated by the many solar-like features exhibited by these oscillations, we carry out an analysis of the large-scale zonal dynamics. We demonstrate that simulated torsional oscillations are not driven primarily by the periodically-varying large-scale magnetic torque, as one might have expected, but rather via the magnetic modulation of angular-momentum transport by the large-scale meridional flow. This result is confirmed by a straightforward energy analysis. We also detect a fairly sharp transition in rotational dynamics taking place as one moves from the base of the convecting layers to the base of the thin tachocline-like shear layer formed in the stably stratified fluid layers immediately below. We conclude by discussing the implications of our analyses with regards to the mechanism of amplitude saturation in the global dynamo operating in the simulation, and speculate on the possible precursor value of torsional oscillations for the forecast of solar cycle characteristics.

## Section 1 Introduction

Differential rotation of the solar surface was first noted in the seventeenth century by Christoph Scheiner on the basis of his extensive sunspot observations. In his 1632 Rosa Ursina, Scheiner states that sunspots move more slowly the farther away they are from the solar Equator, and even concludes that “…From this phenomenon is drawn the strongest argument for a fluid surface of the Sun.” (Rosa Ursina, p. 559; as cited and translated by ?, ?, p. 440). Rediscovered in the mid-nineteenth century by R.C. Carrington and G. Spörer and soon thereafter extended to high latitudes by Doppler measurements, solar differential rotation has now been mapped deep into the Sun by helioseismology [Christensen-Dalsgaard, 2002, Howe, 2009]. From the frequency splitting of acoustic eigenmodes of varying azimuthal orders, it has now been shown that the surface latitudinal differential-rotation pattern, with the solar Equator rotating approximately 30% faster than the Poles, persists throughout the bulk of the solar convective envelope, down to (with the Sun’s radius), where it abruptly vanishes across a thin spherical shear layer, known as the tachocline, located immediately beneath the core–envelope interface. The underlying stably stratified core appears to be rotating rigidly (or nearly so) down to , at a rate equal to that of the surface mid-latitudes (e.g. ?, ?, and references therein).

This internal differential-rotation pattern has remained generally steady since the first helioseismic rotational inversions carried out in the late 1980s; but not exactly steady. Rotational torsional oscillations were first noted in surface Doppler measurements [Howard and Labonte, 1980], and later shown by helioseismology to extend all the way to the base of the Sun’s convective envelope. The torsional oscillation signal reaches only a few nHz in amplitude (about 0.5% of the rotational frequency), and peaks at high latitudes and in surface and subsurface layers. The oscillations develop at twice the magnetic-cycle frequency and, at mid- to high latitudes (where the signal is the strongest), retain the same phase at all depths (e.g. Figure 26 in ?, ?). More elaborate phasing patterns occur with latitudes, with two diverging “branches” of faster rotating fluid appearing at mid-latitudes around solar-activity minimum: one migrating all the way down to the Equator in the span of two full activity cycles, the other migrating poleward to cause a marked spin-up of the polar region peaking at around the time of activity maximum (see Figure 25 in ?, ?, and accompanying discussion).

Numerous models have been suggested to explain the observed behavior of solar torsional oscillations, the vast majority relying directly or indirectly on the Lorentz force associated with the Sun’s magnetic field. ? (?) in Section 9.5 gives a succinct overview of these various theoretical explanations, which turn out to be difficult to confirm or refute on the basis of extant observations. Torsional oscillations having higher amplitudes near the surface and at high latitudes are certainly to be expected; subjected to a torque of a given magnitude (and of whatever origin), a ring of fluid centered on the solar rotation axis will experience greater angular acceleration if located high up in the envelope since its moment of inertia will be reduced through the lower density; and at a given depth, that same moment will also be smaller at higher latitude because of the shorter moment arm, yielding again greater angular acceleration. At any rate, both the good phase locking with the magnetic-activity cycle as measured, e.g. through the sunspot number, and the close tracking of the equatorially migrating band of rotational acceleration with the activity belts, point toward a close dynamical link between torsional oscillations and the cyclically varying large-scale solar magnetic field. This has in fact remained the favored explanation ever since the discovery of torsional oscillations [Schüssler, 1981, Yoshimura, 1981]. Interest in this possible dynamical linkage has in fact recently ramped up, due primarily to the curious observation that the poleward branch of the torsional oscillations, due to appear in the final years of Cycle 23, has failed to show up as anticipated. Taken together with other peculiar features of the extended activity minimum having followed Cycle 23, this has prompted speculations regarding the possibility that the Sun is about to enter a phase of strongly suppressed magnetic activity, perhaps akin to the Maunder Minimum [Hill et al., 2011].

From a dynamical point of view, the simplest hypothesis would be to assume that torsional oscillations are directly driven by the Lorentz force associated with the cyclic large-scale magnetic component that we associate with the magnetic-activity cycle, acting on the zonal flow as a time-varying perturbation of the global hydrodynamical (HD) balance setting the form of solar internal differential rotation.
Such a balance would involve Reynolds stresses, Maxwell stresses, and angular-momentum advection by the meridional flow within the convection zone^{1}^{1}1Some level of dynamical coupling to the underlying stably stratified radiative core
may also play an important role; on the possible participation of the outer radiative core in setting the angular-momentum balance within the convection zone; see the prescient analysis presented in ? (?).. As will become apparent in what follows, the situation may well be far more complex.

In this article, we present an analysis of the dynamics and energetics of torsional oscillations arising in a global magnetohydrodynamical (MHD) simulation of solar convection producing solar-like cycles in its dynamo-generated large-scale magnetic field. We first (Section 2) give an overview of the simulation itself, together with a description of torsional oscillations arising therein. We then recast the azimuthal component of the momentum equation in conservative form, which allows the study of azimuthal force balance in terms of fluxes of angular momentum and their temporal variations (Section 3). We also examine the energetics of torsional oscillations, and conclude (Section 4) by elaborating on some consequences of our analysis for dynamo saturation, and for the possible use of torsional oscillations as precursors of cycle-amplitude fluctuations.

## Section 2 Numerical Data

### 2.1 The Global Simulation

We use numerical data produced by one of the global implicit large-eddy simulations (ILES) of MHD solar convection of the type presented by ? (?), and ? (?). These remain unique so far in producing an axisymmetric large-scale magnetic-field component undergoing cyclic polarity reversals on a multi-decadal timescale, in a manner similar in many ways to what is observed on the Sun. The underlying mathematical and computational frameworks are described by ? (?), and represent a MHD generalization of the well-proven general-purpose geophysical flow simulation code EULAG (see ?, ?, and references therein). A unique feature of both EULAG and EULAG–MHD is the possibility to delegate all dissipation to the underlying advection scheme, which makes it possible to reach a maximally turbulent state at a given grid size, with stability persisting even when field gradients reach spatial scales commensurate with the computational cell size. The interested reader is referred to ? (?) for further details on algorithmic implementation.

The MPDATA algorithm (?, ?; ?, ?) at the core of both EULAG and EULAG–MHD belongs to the class of advection algorithms known as non-oscillatory forward-in-time (NFT), which also includes a number of advection schemes of wide usage in computational fluid dynamics and numerical astrophysics, e.g. the Flux Corrected Transport algorithm of ? (?) (and its many subsequent variants and elaborations), and the PPM scheme of ? (?). In such schemes, dissipation is introduced at the algorithmic level, rather than as explicit dissipation terms in the governing equations, as in conventional Large Eddy Simulation approaches. The crux in the design of ILES scheme is to keep this numerical dissipation at a minimum, i.e. to ensure it activates only when and where it is needed to maintain locally smooth (i.e. spurious oscillation free) solution. In such schemes the level of implicit diffusivities typically decreases with increasing mesh size, but the absolute level of numerical dissipation also depends on the specific algorithm used. Numerous examples are provided in the volume edited by ? (?), including applications to local and global solar–stellar convection.

One drawback of ILES schemes is the difficulty in calculating, even a posteriori, conventional dimensionless numbers such as the Reynolds or Rayleigh numbers, which complicates comparison with simulations carried out including explicit dissipation.

From the practical point of view, however, their main advantage is that nonlinearly stable turbulent simulations can be produced on relatively coarse grids (?, ?; ?, ?; ?, ?; ?, ?), which, in particular, allows long temporal integration, as required in the study of behaviors such as magnetic cycles, developping on timescales very much longer than the turbulent turnover time.

The specific simulation segment analyzed in what follows spans 180 years, and is executed on a relatively small mesh of size in radius latitude longitude. The spatial domain is a thermally forced thick spherical shell () rotating initially rigidly at the solar rate, convectively unstable in its outer two thirds (). Small perturbations in radial flow speed and magnetic field are introduced in the convectively unstable portion of the domain at .

The foregoing analysis begins with the four-dimensional datasets (three spatial dimensions plus time) returned by the simulation. The first step is to extract the axisymmetric components of the total flow and magnetic field. As shown by ? (?) through modal decomposition, these axisymmetric components evolve on a timescale much longer than their non-axisymmetric counterparts, and can thus be legitimately considered as a distinct dynamical entity. We therefore express the total flow [] and magnetic field [] as

(\theequation) | |||||

(\theequation) |

where

(\theequation) | |||||

(\theequation) |

are the axisymmetric large-scale components, calculated by zonal averaging:

(\theequation) |

Note that under these definitions, , so that the non-axisymmetric contributions of the flow and field play the role of the “small scales” in mean-field theory.

Figure 1 offers four views of the large-scale (axisymmetric) toroidal magnetic-field component evolving over the timespan of the simulation. The top two panels (a) and (b) show time–latitude cuts, the first extracted at the depth coinciding with the base of the convecting layers (), and the other near its top (). The bottom two panels show time–radius cuts extracted at (c) low and (d) mid-latitudes in the southern hemisphere. Regular polarity reversals of the large-scale magnetic field stand out prominently in these diagrams, here on a period of about 36 years for each half cycle (equivalent to a sunspot cycle); thus the magnetic cycle period in this simulation is a little over three times longer than the years observed on the Sun. The large-scale toroidal component is antisymmetric about the equatorial plane, in agreement with Hale’s polarity laws, and peaks at mid-latitudes (panel (a)) and immediately beneath the core–envelope interface (panel (d)); this latter property is in line with the need to form and store in the convectively stable layer the toroidal magnetic flux ropes that, upon buoyancy-driven destabilization and emergence, will give rise to sunspots (see ?, ?, and references therein). The subsurface time–latitude diagram on panel (b) and time–radius diagram on panel (c) also show a hint of a secondary dynamo mode, of much shorter period and lower amplitude than the primary mode, producing what looks like an oscillation superimposed on the more slowly evolving magnetic component pervading the bulk of the domain. Interestingly, a similar combination of long- and short-period dynamo modes was also observed in the spherical wedge simulations of ? (?). This intriguing dynamo feature will be revisited in what follows.

Other features of this simulation are discrepant with respect to the solar cycle, besides the period. Most notably, the toroidal magnetic component at the core–envelope interface, where sunspots are presumed to originate, peaks at too high latitudes compared to the sunspot butterfly diagram, and only shows a hint of equatorward migration. Moreover, analysis of the poloidal large-scale component reveals that the latter oscillates in phase with the deep-seated toroidal component, whereas in the Sun a phase lag of is inferred. The specific simulation we are using for the foregoing analysis develops a slow phase drift between hemispheres, which eventually leads, after some 300 years, to a switch to a non-axisymmetric large-scale dynamo mode, a fascinating dynamo behavior in and of itself. Despite these departures with respect to observed solar behavior, the presence of a well-defined cyclic behavior in the large-scale magnetic field offers a unique opportunity to investigate quantitatively the magnetic back reaction on large-scale flows building up in the simulation and observed in the Sun, in particular differential rotation.

### 2.2 The Mean Differential Rotation

Mechanically speaking, solar differential rotation is driven primarily by Reynolds stresses arising through rotation-driven anisotropies in convective turbulence and angular-momentum transport by meridional flows. Helioseismology has now mapped, with good accuracy, differential rotation throughout the bulk of the solar convection zone and upper radiative core (e.g. Figure 18 in ?, ?). If one excludes subphotospheric layers, the primary rotational gradient in the solar convection zone is latitudinal, with the rotational frequency of equatorial regions exceeding that of polar regions by about 30%. This latitudinal gradient vanishes at the interface between the convection zone and the underlying radiative core, across a thin shear layer known as the tachocline.

With simulations computed in a reference frame rotating at the mean solar rate (), differential rotation can be directly computed from the zonally averaged -component of the flow velocity. We first compute the mean differential rotation profile by temporally averaging over the full temporal extent of the simulation, which is in a statistically equilibrated state. We have carried out this averaging exercise for the MHD simulation of Figure 1, as well as a purely HD convection simulation operating under the same forcing regime and rotation rate, and computed using the same mesh size. The results are shown on Figure 2, expressed as angular velocity , with the latitude. The left panels show isocontour maps, with corresponding radial cuts plotted on the right panels, on the same scales to allow quick visual comparison of the two simulations. Both are characterized by equatorial acceleration, but with isocontours too closely aligned with the rotation axis, and too concentrated towards the middle of the convection zone at low latitudes, as compared to the helioseismically-inferred solar internal differential rotation. These features are in fact typical of these types of simulations (e.g. Figure 9 in ?, ?; Figure 1 in ?, ?; Figure 3 in ?, ?), unless a latitudinal gradient in the heat flux is artificially imposed at the base of the domain (see ?, ?). Nonetheless, the differential rotation characterizing the HD simulation (top row) shows some remarkably solar-like features, notably the magnitude of the Pole-to-Equator angular velocity contrast, and, in particular, a thin tachocline-like rotational shear layer located immediately beneath the core–envelope interface. The thinness of this layer ( here) is a direct reflection of the very low dissipation levels characterizing this simulation (cf. Figure 1 in ?, ?).

The resemblance to solar differential rotation degrades, however, upon moving to the MHD simulation (bottom row in Figure 2). The Pole-to-Equator angular velocity contrast is now reduced by a factor of three as compared to the HD simulation, and the latitudinal gradient has all but vanished at mid to high latitudes. Such strong magnetic backreaction on the mean differential rotation is in fact typical of these types of global MHD convection simulations, as shown already by ? (?) and ? (?). A residual tachocline remains, in the sense that the weak convection zone latitudinal differential rotation again vanishes across a thin shear layer beneath the core–envelope interface. Although quite weak, this differential rotation remains important for the operation of the dynamo, as the analysis of a similar simulation carried out in ? (?) shows that it contributes approximately equally with the turbulent electromotive force to the regeneration of the large-scale toroidal magnetic field near the core–envelope interface.

Another important difference between the differential rotation profiles characterizing the HD and MHD simulations is that the former is temporally steady once the simulation has attained a statistically stationary state. The dotted lines bracketing each radial cut on the right panels of Figure 2 indicate the deviations of the zonal averages about their temporal mean plotted on the left panels. These are quite small in the HD simulation, except near the poles. The MHD simulation, on the other hand, exhibits somewhat larger one-sigma deviations at most latitudes, but these now reflect the presence of spatiotemporally coherent cyclic variations superimposing themselves on the mean rotational profile. We now turn to the characterization of these torsional oscillations.

### 2.3 The Torsional Oscillations

Torsional oscillations are best visualized by subtracting the temporally averaged rotational frequency profile of Figure 2 from the corresponding zonally-averaged rotational frequency at each time step:

(\theequation) |

where the overbar denotes temporal averaging over the timespan of the simulation:

(\theequation) |

amounting to the zonal and temporal average of the full azimuthal velocity [] numerical data set, with the length of the simulation. The result of this procedure is shown on Figures 3 and 4, which display respectively time–latitude diagrams at four fixed depths, and radius–time diagrams at four fixed latitudes. Several features visible in these plots are noteworthy: i) A cyclic signal is clearly present at all depths and latitudes so sampled, at twice the frequency characterizing the magnetic cycle (cf. Figure 1); ii) The torsional oscillations peak in amplitude at high latitudes and in the surface and subsurface layers, reaching there nHz; iii) At mid- to high latitudes, the oscillations show a phase approximately independent of depth. The oscillations reach their peak prograde phase (i.e. peaking at positive values) at about the peak of the magnetic cycle. All of these features are remarkably solar-like, as can be inferred from comparison with a similarly formatted diagram in Howe (2009; cf. Figure 26 to Figure 4 herein).

At first glance, the surface spatiotemporal pattern of the torsional oscillations is not particularly solar-like (cf. Figure 25 of ?, ?, and Figure 3(d) herein). Near the surface, a strong and rather complex oscillatory signal is present at low latitudes, arising from the superimposition of an oscillation associated with the large-scale magnetic cycle with a second, characterized by higher frequency and restricted to the subsurface equatorial regions. Examination of the simulation reveals that this is associated with a secondary dynamo mode feeding on the strong latitudinal shear present in the outer half of the convection zone at low latitudes (see Figure 2; also Figure 1(b) and (c)).

However, the dissimilarities greatly diminish if one focuses on the pattern present at latitudes higher than the equatorial magnetic “activity belts”, as defined by the time–latitude distribution of the large-scale zonal magnetic component at the core–envelope interface (Figure 1(a)), where sunspots are presumed to originate. Figure 5 illustrates the idea. It is essentially a closeup of the northern hemisphere portion of Figure 3(c), on which have been superimposed a few isocontours of mean toroidal magnetic field strengths at the core–envelope interface (cf. Figure 1(a)) for the first three half-cycles in the simulation. This stretching procedure “renormalizes” the activity belts to low latitudes, and allows a comparison that is arguably more relevant to the patterns observed on the Sun. If one accepts this stretch at face value, then the comparison becomes actually quite good (cf. e.g. Figure 25 in ?, ?). In particular, the torsional acceleration () is seen to begin at high latitudes (here ) at about the time of magnetic-polarity reversal (akin to solar minimum here), and develops in two diverging branches, one propagating poleward and the second, of lower amplitude, propagating equatorward. In the Sun, this second branch requires two activity cycles to reach the Equator, while here it does so in only one cycle; this may be a reflection of the fact that our activity belts are located at too high latitudes, but this remains to be demonstrated through further simulations. It is also quite remarkable that this latitudinal double-branch pattern in the torsional oscillations persists all the way to the base of the convecting layers in the simulation. Careful examination of Figure 4(d) also reveals that an oscillatory signal even penetrates all the way through the underlying stably stratified layer to the base of the computational domain (), although this signal may reflect the excitation of gravity waves.

At any rate, the results presented here demonstrate clearly that this simulation generates a global torsional oscillation pattern that is solar-like in a number of ways. Remaining discrepancies notwithstanding, we have in hand a unique “virtual laboratory” allowing a quantitative and fully dynamical investigation of the mechanism(s) driving these torsional ocillations. This is the topic to which we now turn.

## Section 3 The Dynamical Drivers of Torsional Oscillations

### 3.1 The Zonal Momentum Equation

Torsional oscillations are generated because the (magneto)dynamics of the convecting layers lead to a systematic, time-dependent redistribution of angular momentum that is driven, directly or indirectly, by the magnetic cycle. In parallel to the MHD equations solved in the simulation, the foregoing analysis focuses on the azimuthal component of the Navier–Stokes equations, including the Lorentz and Coriolis forces, and written under the anelastic approximation. The latter implies, in particular, that mass conservation reduces to , with .

The starting point of our analysis is to recast the azimuthal momentum equation in conservative form involving fluxes of angular momentum. We begin by applying the azimuthal averaging operator previously defined in Equation (2.1):

(\theequation) |

with the pressure profile, , where is the velocity field of the plasma in the rotating frame of the Sun, is the mean angular velocity of the Sun, is a radial vector locating a given fluid element in a spherical coordinate system with origin at the Sun’s center, and once again the variable denotes the latitude. Now, the azimuthal averaging operator commutes with derivatives acting on large-scale quantities, and the Lorentz force term can be rewritten as . Moreover, under the axisymmetry assumption, all azimuthal derivatives of averaged quantities automatically vanish, so that the above expression reduces to

(\theequation) |

Use of judicious vector identities also leads to

(\theequation) | |||||

(\theequation) |

Using the scale separation introduced in Equation (2.1) – (2.1), the terms and can be separated between large-scale and small-scale contributions, e.g.:

(\theequation) |

and likewise for . Inserting Equation (3.1) and Equation (3.1) into Equation (3.1) yields

(\theequation) |

One can separate the rigid rotation component by writing , noting that because there is no small-scale contribution to . Since the latter term is also divergence free, we can write

(\theequation) |

The end result of all this is to recast Equation (3.1) in conservative form:

(\theequation) | |||||

The four terms adding up under the divergence on the LHS act as volumetric force densities in the longitudinal direction, and define here the four contributors to zonal dynamics:

(\theequation) | |||||

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

These are, respectively, turbulent Reynolds stresses, Coriolis force acting on the meridional flow, turbulent Maxwell stresses (small-scale Lorentz force) and the magnetic torque (Lorentz force associated with the large-scale magnetic component). Note that viscous stresses do not appear explicitly here, as our simulation is of the implicit large-eddy type, i.e. it does not include explicit dissipative terms in the momentum equation.

### 3.2 Angular-Momentum Fluxes

As per Equation (\theequation), the mean radial and latitudinal angular-momentum fluxes are given by

(\theequation) | |||||

(\theequation) | |||||

Following ? (?), we first examine the global rotational dynamics by computing from the simulation output the fluxes of angular momentum integrated across spherical shells or conical wedges centered on the rotation axis:

(\theequation) | |||||

(\theequation) |

so that is the net angular-momentum transport rate through shells of different radii, and through cones tangent to different latitudes.

In order to disentangle the various physical contributions to angular-momentum transport, these integrals are computed separately for the four distinct contributions to the total angular-momentum fluxes. The result of this procedure is shown in Figure 6, where the fluxes have also been temporally averaged over the extent of the simulation. The procedure was also carried out for the same parent HD simulation whose mean rotational frequency is plotted on Figure 2 (top). In this HD simulation, the only contributors to zonal dynamics are the Reynolds stresses and Coriolis force acting on the meridional flow.

In the HD simulation, the Reynolds stresses and Coriolis force are seen to act
in opposition at all latitudes and depths. This is precisely what one would expect
for a stationary rotational state to ensue. With all fluxes vanishing
at domain boundaries, the fact that these two contributions
do not add up to zero reflects the presence, in the simulation, of a dissipative
force associated with the numerical scheme, which here contributes
to the zonal dynamics, especially at low latitudes in the convection zone.
This
can be traced in part to the very strong rotational shear which builds up there
in the HD simulation
(cf. top panels in Figure 2)^{2}^{2}2Such a dissipative force,
of purely numerical origin, is the hallmark of ILES simulations.
In the specific context of EULAG, the hydrodynamical parent of EULAG–MHD,
it has been shown by ? (?)
to mimic bona fide subgrid scale parametrization in the context
of the planetary boundary layer,
and, more importantly, to effectively vanish
if an appropriate explicit subgrid model is introduced in the model.

Turning to the MHD simulation, the most obvious feature is perhaps the fact that all four potential fluxes of angular momentum contribute more or less equally to the global rotational dynamics. Moreover, the presence of magnetic fields has greatly altered the non-magnetic fluxes of angular momentum. This is particularly the case for the radial flux of angular momentum associated with the Coriolis force (cf. green curves on panels (a) and (b)), and latitudinal fluxes by Reynolds stresses (cf. blue curves on panels (c) and (d)), which both undergo reversals in direction over a substantial portion of their spatial range when going from HD to MHD. Also noteworthy in this MHD simulation, the Maxwell stresses associated with the small-scale magnetic field contribute more or less equally to the Lorentz force associated with the large-scale, cyclic magnetic field. Finally, the presence of a magnetic field leads to a significant rotational coupling between the convection zone and underlying stably stratified core, extending deep into the latter. Such a coupling is almost entirely absent in the HD simulation. The high degree of (anti)symmetry about the equator apparent on panels (c) and (d) is a true feature of these simulations, as no averaging of hemispheres has been carried out here.

It is particularly interesting to compare Figures 6(b) and (d) to the corresponding diagrams presented by ? (?), namely their Figure 10. Note that in ? (?), their equivalent of our integrated radial fluxes are divided by and they are using CGS units instead of SI units. Moreover, we use latitudes rather than polar angles, so that on Figures 6(c) and (d) herein a positive latitudinal flux implies northward transport, while in ? (?) northward transport correspond to negative flux values. Their simulation differs from ours in three important ways : i) it covers only the convecting layers; ii) it includes substantial explicit viscosities and magnetic diffusivities throughout the simulation domain; and iii) it does not generate a large-scale magnetic component. Even though their mean differential rotation is qualitatively similar to ours (equatorial acceleration, polar deceleration, tendency towards cylindrical isocontours in equatorial regions), major differences exist in the underlying rotational dynamics. Viscous forces play a major role in the radial transport of angular momentum in their simulation, teaming with Maxwell stresses to offset the Coriolis force and Reynolds stresses throughout the whole convection zone. In our simulation, this dynamical balance (minus explicit viscous force) materializes only in the lower third of the convecting layers, with Reynolds and Maxwell stresses acting in opposition to the Coriolis force higher up. In both simulations, the latter leads to a net upward transport of angular momentum, and Maxwell stresses drive a downward transport.

For the simulation analyzed on Figure 10 of ? (?), viscous diffusion makes a lesser contribution to angular momentum transport in the latitudinal direction than it does in the radial direction. There are now more similarities with our latitudinal fluxes, Reynolds and Maxwell stresses now opposing each other at most latitudes. The primary difference, besides the presence of a significant large-scale magnetic torque in our simulation, is found with the net latitudinal angular-momentum transport by the meridional flows, which is equatorward at low to mid-latitudes in our simulation, but poleward in the ? (?) simulations. In our MHD simulation, except very near the Equator, the meridional flow and Reynolds stresses team up to sustain equatorial acceleration. These two are resisted by magnetic forces, consistent with the Pole-to-Equator angular velocity contrast being some three times smaller in the MHD simulation than its purely HD parent simulation (cf. top and bottom panels on Figure 2).

The differences between these two sets of simulations most likely do not arise exclusively from the presence of a large-scale cyclic magnetic field in our simulations, as angular-momentum fluxes calculated in purely HD versions of the ? (?) simulations (see ?, ? and ?, ?) also differ markedly from the HD balance depicted on Figures 6(a) and (c) herein. Looking at case (c) in Figure 11 from ? (?), which is the most turbulent, and comparing it with our Figures 6(a) and (c), we observe very different patterns for both the Reynolds stresses and the Coriolis force acting on the meridional flow. The radial fluxes show distinct depth variations, particularly in the middle of the convection zone, but the difference is most striking in the latitudinal fluxes distributions. At low latitude [] in Figure 6(c) and Figure 11(c) from ? (?), both the Reynolds stresses and the meridional-circulation contributions are of opposite sign when comparing both simulations. Additionally, there is a sign change around in our simulation that has no counterpart in theirs. ? (?) present another set of HD simulations, all strongly turbulent but where they also vary the rotation rate. More specifically, in their Figure 9 they show angular-momentum fluxes for two cases: a star which has a solar-like rotation rate in panels (a) and (b), and one with a rotation rate five times greater than the Sun in (c) and (d). While panels (a) and (b) remain similar to the aforecited equivalent plots in ? (?), the more rapidly rotating case (panels (c) and (d)) reveals yet again a distinct dynamical balance. The major players in the radial transport of angular momentum are the Reynolds stresses and the viscous transport terms, with the Coriolis force exerted on the meridional circulation playing a lesser role. While latitudinal transport has a very complex profile, the Reynolds stresses and meridional circulation terms usually have the same sign and are counterbalanced by the viscous transport.

Another major difference lies of course with the fact that our dynamical balance is not perfectly stationary, showing instead periodic variations which drives the torsional oscillations visible on Figures 3 and 4. It is therefore also interesting to examine the temporal evolution of those angular-momentum fluxes over a magnetic half-cycle. This is carried out on Figure 7, which shows the evolution of each individual flux component, from the beginning of the second half-cycle on Figure 1 to its end, i.e. from one minimum to the next, at a nine-year temporal cadence, as indicated by the vertical line segments on Figures 1(a) and (d). Each panel also reproduces the temporal average of the corresponding contribution to the angular-momentum transport rate over the full simulation duration (in black), taken directly from Figure 6(b).

The magnetic-torque contribution (panel (d)) shows large variations about its temporal average, which is of course expected in view of the cyclic evolution characterizing the large-scale magnetic field. Far less expected a priori, however, is the fact that all other flux contributions also undergo similar variations in the bulk of the convecting layers. This is even the case for the nominally non-magnetic contributions, namely the turbulent Reynolds stresses and Coriolis force acting on the meridional flow. The latter’s temporal variations show little spatial coherence, whereas Reynolds stresses vary largely in unison at all depths. A similar situation arises with the magnetic contributions, with the magnetic torque showing a complex spatio–temporal behavior, while temporal variations in Maxwell stresses are very well-correlated spatially. Also noteworthy, flux contributions associated with the small-scale flow and magnetic field — Reynolds and Maxwell stresses — fall rapidly with depth at and below the core–envelope interface, where they show almost no temporal variability. Consequently, the time-dependence of the rotational coupling between the convection zone and underlying stably stratified fluid layers is driven by the interplay between the strongly time-varying contributions of large-scale magnetic torques and angular-momentum advection by the meridional flow.

Another interesting aspect of our results relates to the rotational coupling between the convection zone, where differential rotation is generated, and the underlying stably stratified fluid layers. Figure 8(a) shows time series of the radial fluxes of angular momentum, computed via Equation (\theequation) separately for each of its four contributions, as labeled. At the core–envelope interface proper (), the net transport of angular momentum across the corresponding spherical shell is directed upward, and it is driven primarily by the Coriolis force acting on the meridional flow, and resisted by the magnetic forces. Although the magnitude of the large-scale magnetic torque varies cyclically in phase with the magnetic cycle, as one would have expected, the Coriolis term does also, which results in a net upward flux of angular momentum (in black) that does not show a well-defined cyclic signal. This general pattern is maintained down to in the stable layers, but with the disappearance of the HD forces further below, the magnetic terms take over completely, as shown on Figure 8(b). The large-scale torque now drives an upward angular-momentum flux, and is opposed by the Maxwell stresses. Cyclic variations on the magnetic cycle frequency can still be detected, but quasi-cyclic modulations on shorter periods are also apparent at these depths.

A similar dynamical transition as a function of depth in the stable layer is also present in the latitudinal angular-momentum fluxes, as illustrated on Figures 8(c) and (d). These time series result from the calculation of the four individual contributions on the RHS of Equation (\theequation) at latitude and two different depths, as labeled. At the mid-latitude core–envelope interface (panel (c)), the latitudinal flux is directed equatorward and dominated by the meridional flow contribution, which shows a strong cyclic signal in phase with the magnetic cycle. This is associated with a strong driving of the meridional flow by the large-scale magnetic field (on this point see also ?, ?), but with regards to the zonal dynamics, the large-scale magnetic torque here drives angular momentum poleward, as expected from the shearing of a latitudinally oriented poloidal large-scale magnetic field by a latitudinal differential rotation characterized by equatorial acceleration. Moving inward, already at (not shown) the magnetic torque has reversed and now drives angular momentum equatorward, and by (panel (d)) the contribution of the meridional flow has vanished and the large-scale magnetic torque is the sole driver of the equatorward angular-momentum flux, which remains strongly modulated by the magnetic cycle. The Reynolds and Maxwell stresses are minor contributors to the latitudinal flux of angular momentum at all depths in the stable layer.

The downward penetration of the large-scale magnetic fields within the stable layer allows a rough estimate of the average numerical diffusivity characterizing the simulation in that part of the domain. The large-scale toroidal magnetic field peaks around , oscillates on a -year full magnetic period, and penetrates inward by (see Figure 1(d)). Equating this observed penetration distance to the electromagnetic skin depth leads to an estimate of the effective magnetic diffusivity m s. This is quite small, as expected given that large velocity and magnetic field gradients, which determine the level of implicit diffusivities within our computational framework, diminish rapidly upon moving downward into the stable layers. One must bear in mind, however, that the effective diffusivities are likely much higher in the turbulent environment of the convecting layers.

### 3.3 Volumetric Force Densities

We now turn to explicit calculations of volumetric force densities, by taking the divergence of the various contributions to the total angular-momentum flux, as appearing on the LHS of Equation (\theequation), at each grid point in the meridional plane. In both HD and MHD simulations, temporally averaging each set of grid point values over four magnetic half-cycles again produces patterns with a high degree of symmetry with respect to the Equator. Comparing HD and MHD simulations reveals the most pronounced differences at mid- to high latitudes, where even the nominally HD forces — Reynolds stresses and Coriolis force acting on the meridional flow — show large differences in their spatial distributions. This is particularly striking in the Coriolis term, which tends to accelerate (decelerate) the zonal flow in the outer (inner) half of the convecting layers of the HD simulation, and shows the opposite pattern in the MHD simulation. In the MHD simulation, Maxwell and Reynolds stresses tend to oppose each other at most locations in the meridional plane, and the two magnetic contributions are the sole significant contributors within the underlying stable fluid layer, consistent with Figure 6.

Figure 9 offers a more focused look at the torsional-oscillation dynamics at high latitudes, in the form of time series of the various azimuthal-force components extracted at a specific grid point located at depth and latitude (southern hemisphere). At first glance, the Reynolds and Maxwell stresses contribute very little to the zonal dynamics at this location, which is primarily driven by the large-scale Lorentz force (red) and Coriolis force (green). However, these two contributions are of similar magnitudes but strongly anticorrelated, and so nearly cancel each other. The resulting total force (black) is then of much smaller magnitude than either of these two contributions, and comparable again to the lower amplitude Maxwell stresses (orange).

Figure 10(a) shows time series of this volumetric force density (black) and zonal flow speed (blue), both now averaged radially over the convection zone () and latitudinally over a band of angular width centered on the South Pole. The red line is a proxy for the cycle amplitude, constructed by integrating the magnetic energy over a spherical shell () straddling the core–envelope interface. All three time series have been normalized to their peak unsigned amplitude, for plotting purposes. Despite strong temporal fluctuations, the plot reveals that the net zonal force peaks throughout the rising phases of each magnetic cycle, leading the rise in the zonal prograde flow. The latter peaks shortly before cycle maximum, and reaches its lowest speed at times of cycle minimum (cf. Figures 1 and 4). This pattern is illustrated differently on Figure 10(b), which shows a correlation plot of peak longitudinal velocity amplitude versus magnetic energy, calculated by averaging temporally over a two-year temporal window centered on either cycle maximum or minimum. The calculation is carried out independently in each hemisphere, and reveals a good correlation ().

Figure 11 offers yet a different look at the torsional-oscillation dynamics. The curves are trajectories in a two-dimensional phase space defined in terms of the zonally averaged zonal flow deviation about its temporal mean over the simulation timespan (horizontal), versus the zonally-averaged latitudinal flow deviation about its own temporal mean (vertical). Four such trajectories are shown, for meridional plane grid points located at the subsurface high latitudes in panels (a) and (b), and mid-high-latitude core–envelope interface in panels (c) and (d).

In the high-latitude subsurface layers (panels (a) and (b)), both flow perturbations are strongly correlated in time, with the poleward latitudinal flow varying in phase with the prograde rotational acceleration. That local rotational acceleration (deceleration) should correlate in this manner to the variation of the latitudinal flow component, consistent with conservation of angular momentum in an axisymmetric fluid ring symmetrically contracting (stretching) as it gets displaced towards (away) from the rotation axis by the latitudinal flow. This suggests that, in this high-latitude location, the zonal dynamics are “enslaved” to the meridional flow dynamics. Interestingly, a similar situation was observed by ? (?) when analyzing the torsional oscillations generated by his pioneering global MHD simulation of solar dynamo action, even though magnetic driving of the meridional flow perturbation could not be firmly established. This is all the more remarkable given that in the ? (?) simulations, angular-momentum transport by the meridional flow played little significant role in sustaining the average differential rotation profile.

Enslavement of zonal dynamics to meridional flow variations does not hold everywhere, however, as evidenced by panels (c) and (d) of Figure 11. These diagrams are now constructed from time series of zonal- and latitudinal-flow variations extracted at latitude at the core–envelope interface. The corresponding phase diagrams of zonal- and latitudinal-flow residuals are now markedly different from their high-latitude subsurface counterparts on panels (a) and (b), with the two flow residuals now varying cyclically but out of phase with one another, with a phase lag . This indicates that the zonal dynamics cannot be reduced to angular-momentum conservation in a contracting (expanding) fluid ring, and results from a more complex interplay of time-varying direct and indirect magnetically mediated forcing (see also the companion analysis presented in ?, ?).

### 3.4 Energetics of the Torsional Oscillations

The conclusion drawn above can be further substantiated through an analysis of the flux of energy to and from the various energy reservoirs defined by the flow and magnetic field. The evolution equation for the kinetic energy density of the flow [] takes here the form:

(\theequation) |

where the volumetric force densities associated with the flow [], magnetic field [], and buoyancy [] have been schematically grouped on the RHS. With , both terms on the RHS correspond to the volumetric work done by the various forces on, or against, the flow u. In what follows, we are interested in the energy flow in association with the zonal dynamics, so we will set , , , and , as defined by Equations (\theequation) — (\theequation). Table 1 lists the total power associated with these zonal volumetric force components (left column), averaged over the full 180-year time span of the simulation segment and integrated spatially over the southern polar caps (middle column), and over the full domain (right column).

Over the full simulation domain, Reynolds stresses and the Coriolis force inject energy into the zonal flow, each contributing more or less equally, and are resisted by both magnetic contributions, each again contributing approximately equally to the draining of zonal kinetic energy. Moreover, all four of these power contributions turn out to be fairly steady in time. This state of affairs reflects primarily the sustenance of the mean differential rotation (Figure 2).

A different picture emerges, however, if one focuses on high latitude regions, where torsional oscillations have their highest amplitudes (see Figure 4). Figure 12 shows time series of each power-density contribution from Equation (\theequation) integrated over a conical volume going from the South Pole to latitude . Energy injection into the torsional oscillations is now dominated by the Coriolis force exerted on the meridional circulation, with peak energy-transfer rates usually occurring at the peak phase of each magnetic half-cycle. Reynolds stresses still contribute significantly, but only around the times of magnetic-polarity reversals. The most important agent systematically extracting energy from the zonal flow is now the large-scale magnetic torque, doing so in a very well-defined cyclic fashion. This reflects the action of the azimuthal Lorentz force associated with the shearing, by the torsional oscillations, of the strong polar magnetic fields generated in this simulation.

The rotational energy extracted by the large-scale magnetic torque represents an energy input into the large-scale magnetic field. In the absence of explicit Ohmic dissipation, magnetic energy evolves according to:

(\theequation) |

where is the Poynting electromagnetic energy flux, and the second term on the RHS is the direct counterpart of the same term appearing on the RHS of Equation (\theequation), except of course for the sign. Torsional oscillations, far from being directly driven by large-scale magnetic torques, instead divert energy back into the large-scale magnetic field through the agency of magnetic torques. The conclusion is that direct magnetic driving of torsional oscillations does not represent a saturation mechanism for the global dynamo operating in this simulation. However, the magnetic cycle does drive large fluctuations in the meridional flow, and the Coriolis force acting on this cyclically forced flow turns out to be the main driver of torsional oscillations. Indeed, the energy analysis presented in Section 5 of ? (?) indicates that magnetic driving of the latitudinal flow is the primary sink of magnetic energy in this simulation. Schematically, the energy flow is thus of the form:

Magnetic energy meridional flow torsional oscillations (numerical) dissipation.

The most important “take-home” message of the above analyses is the following: torsional oscillations are not driven by a cyclic, magnetically mediated perturbation superimposing itself on an otherwise steady hydrodynamical zonal balance. Instead, all angular-momentum flux contributions, including those of a purely HD nature, are strongly modulated by the magnetic cycle. In these simulations, torsional oscillations are a fully nonlinear and truly MHD phenomenon.

## Section 4 Concluding Remarks

In this article, we have carried out a focused analysis of one of the implicit large-eddy MHD simulations of solar convection computed by ? (?) (see also ?, ?). To the best of our knowledge, these simulations remain unique in generating a spatially well-organized large-scale magnetic-field component undergoing regular polarity reversals in a manner resembling in many ways the solar magnetic cycle. We have shown that a well-defined rotational torsional oscillation signal is present in the simulation, showing a high degree of similarity with those observed in the Sun, including: i) frequency twice that of the magnetic cycle; ii) greater amplitudes in polar and subsurface regions, peaking at a few nHz; iii) peak prograde phase coinciding approximately with the peak in large-scale magnetic field; iv) diverging double-branch latitudinal structure; v) depth-independent phase at most latitudes.

We investigated rotational dynamics by first computing from the simulation output the various contributions to angular momentum fluxes, in the MHD simulation as well as a parent, unmagnetized simulation otherwise operating under the same numerical and physical parameter settings. We could show that in the MHD simulation, the presence of a large-scale cyclic magnetic field drives torsional oscillations not just directly through the associated large-scale magnetic torque, but also indirectly by modulating the other forces influencing zonal dynamics, most notably the transport of angular momentum by meridional flow. In fact, all force components driving the zonal flows undergo cyclic variations driven by the magnetic cycle, including the nominally “small-scale” Reynolds and Maxwell stresses.

We also examined the dynamical character of the rotational coupling between the convecting layers and underlying stably stratified fluid layers. Because of the low-dissipative properties of the numerical scheme underlying the simulations, a tachocline-like shear layer builds up immediately beneath the nominal base of the convective layers, within which significant radial and latitudinal cyclically-varing fluxes of angular momentum develop. The net radial flux is determined by competition between the meridional flow and Maxwell stresses within the tachocline () and by the magnetic torque and Maxwell stresses below. The latitudinal flux shows a strong cyclic signal, in phase with the magnetic cycle. It is dominated by the meridional flow within the tachocline, and by the large-scale magnetic torque below. The upper part of the stable layer is here an important player in setting the global cycle of angular-momentum redistribution — and thus torsional oscillations — within the convection zone (on these issues see also ?, ?).

Turning to a simple analysis of the energetics of torsional oscillations, we could also show that the primary direct power source for torsional oscillations arising in the simulation is the action of the Coriolis force on large-scale meridional fluid motion, with large-scale magnetic torques acting to oppose these oscillations at most phases of the cycle. This is in agreement with a parallel investigation by ? (?), who carried out a similar energy analysis for the meridional flow and could show that magnetic driving of this flow represented the primary sink of magnetic energy. This suggests that saturation of global dynamo action in this simulation occurs through the magnetic driving of flows on large spatial scales, magnetic energy being first diverted into the meridional components, and subsequently, through angular-momentum conservation, into torsional oscillations, where it eventually damps through turbulent stresses.

This state of affairs, should it carry over to the real Sun, has interesting consequences with regards to attempts to use fluctuations in large-scale flows as precursors to the solar-cycle amplitude. More specifically, taken jointly with the results presented by ? (?), our analysis suggests that fluctuations in the meridional flow may be better potential precursors than torsional oscillations, because the bulk of magnetic driving of large-scale flows occurs on and through this flow component. Re-analysis of extant surface Doppler measurements back to the beginning of Cycle 22 [Dikpati et al., 2010, Ulrich, 2010] have indeed revealed significant cycle-to-cycle variations in the surface meridional flow. We are currently pushing our simulations much further in time, as well as under different parameter regimes, which should allow a statistically sound investigation of the precursor potential of torsional oscillations and surface latitudinal flow variations.

Acknowledgements : We are indebted to an anonymous referee for useful comments and suggestions. This work was supported by Canada’s Natural Sciences and Engineering Research Council, Research Chair Program, and by the Canadian Space Agency’s Space Science Enhancement Program (grant # 9SCIGRA-21). PB is also supported in part through a graduate fellowship from the Université de Montréal’s Physics department. PKS acknowledges partial support by the NSF grant OCI-0904599. NCAR is sponsored by the National Science Foundation.

## References

- Boris and Book, 1973 Boris, J.P., Book, D.L.: 1973, Flux-Corrected Transport. I. SHASTA, A Fluid Transport Algorithm That Works. J. Comp. Phys. 11, 38. doi:10.1016/0021-9991(73)90147-2.
- Brown et al., 2008 Brown, B.P., Browning, M.K., Brun, A.S., Miesch, M.S., Toomre, J.: 2008, Rapidly Rotating Suns and Active Nests of Convection. Astrophys. J. 689, 1354 – 1372. doi:10.1086/592397.
- Browning et al., 2006 Browning, M.K., Miesch, M.S., Brun, A.S., Toomre, J.: 2006, Dynamo Action in the Solar Convection Zone and Tachocline: Pumping and Organization of Toroidal Fields. Astrophys. J. Lett. 648, L157–L160. doi:10.1086/507869.
- Brun and Toomre, 2002 Brun, A.S., Toomre, J.: 2002, Turbulent Convection under the Influence of Rotation: Sustaining a Strong Differential Rotation. Astrophys. J. 570, 865 – 885. doi:10.1086/339228.
- Brun, Miesch, and Toomre, 2004 Brun, A.S., Miesch, M.S., Toomre, J.: 2004, Global-Scale Turbulent Convection and Magnetic Dynamo Action in the Solar Envelope. Astrophys. J. 614, 1073 – 1098. doi:10.1086/423835.
- Christensen-Dalsgaard, 2002 Christensen-Dalsgaard, J.: 2002, Helioseismology. Rev. Mod. Phys. 74, 1073 – 1129. doi:10.1103/RevModPhys.74.1073.
- Dikpati et al., 2010 Dikpati, M., Gilman, P.A., de Toma, G., Ulrich, R.K.: 2010, Impact of changes in the Sun’s conveyor-belt on recent solar cycles. Geophys. Res. Lett. 37, 14107. doi:10.1029/2010GL044143.
- Domaradzki, Xiao, and Smolarkiewicz, 2003 Domaradzki, J.A., Xiao, Z., Smolarkiewicz, P.K.: 2003, Effective eddy viscosities in implicit large eddy simulations of turbulent flows. Phys. Fluids 15, 3890 – 3893. doi:10.1063/1.1624610.
- Fan, 2009 Fan, Y.: 2009, Modeling the Subsurface Evolution of Active-Region Flux Tubes. In: Dikpati, M., Arentoft, T., González Hernández, I., Lindsey, C., Hill, F. (eds.) Solar-Stellar Dynamos as Revealed by Helio- and Asteroseismology: GONG 2008/SOHO 21, *CS-416*, Astron. Soc. Pacific, San Francisco, 489.
- Ghizaru, Charbonneau, and Smolarkiewicz, 2010 Ghizaru, M., Charbonneau, P., Smolarkiewicz, P.K.: 2010, Magnetic Cycles in Global Large-eddy Simulations of Solar Convection. Astrophys. J. Lett. 715, L133–L137. doi:10.1088/2041-8205/715/2/L133.
- Gilman, 1983 Gilman, P.A.: 1983, Dynamically consistent nonlinear dynamos driven by convection in a rotating spherical shell. II - Dynamos with cycles and strong feedbacks. Astrophys. J. Suppl. Series 53, 243 – 268. doi:10.1086/190891.
- Gilman and Miller, 1981 Gilman, P.A., Miller, J.: 1981, Dynamically consistent nonlinear dynamos driven by convection in a rotating spherical shell. Astrophys. J. Suppl. Series 46, 211 – 238. doi:10.1086/190743.
- Gilman, Morrow, and Deluca, 1989 Gilman, P.A., Morrow, C.A., Deluca, E.E.: 1989, Angular momentum transport and dynamo action in the sun - Implications of recent oscillation measurements. Astrophys. J. 338, 528 – 537. doi:10.1086/167215.
- Grinstein, Margolin, and Rider, 2007 Grinstein, F.F., Margolin, L.G., Rider, W.J.: 2007, In: Grinstein, F.F., Margolin, L.G., Rider, W.J. (eds.) Implicit Large Eddy Simulation : Computing Turbulent Fluid Dynamics, Cambridge University Press, UPH, Shaftesbury Road, Cambridge, UK. doi:10.1017/CBO9780511618604.
- Hill et al., 2011 Hill, F., Howe, R., Komm, R., Christensen-Dalsgaard, J., Larson, T.P., Schou, J., Thompson, M.J.: 2011, Large-scale Zonal Flows During the Solar Minimum – Where Is Cycle 25?, Bull. Am. Astron. Soc., 1610.
- Howard and Labonte, 1980 Howard, R., Labonte, B.J.: 1980, The sun is observed to be a torsional oscillator with a period of 11 years. Astrophys. J. Lett. 239, L33–L36. doi:10.1086/183286.
- Howe, 2009 Howe, R.: 2009, Solar Interior Rotation and its Variation. Liv. Rev. Sol. Phys. 6, 1. http://solarphysics.livingreviews.org/Articles/lrsp-2009-1/, consulted in May 2011.
- Käpylä et al., 2010 Käpylä, P.J., Korpi, M.J., Brandenburg, A., Mitra, D., Tavakol, R.: 2010, Convective dynamos in spherical wedge geometry. Astronom. Nach. 331, 73. doi:10.1002/asna.200911252.
- Margolin, Smolarkiewicz, and Sorbjan, 1999 Margolin, L.G., Smolarkiewicz, P.K., Sorbjan, Z.: 1999, Large-eddy simulations of convective boundary layers using nonoscillatory differencing. Phys. D Nonlinear Phenom. 133, 390 – 397. doi:10.1016/S0167-2789(99)00083-4.
- Margolin, Smolarkiewicz, and Wyszogradzki, 2006 Margolin, L.G., Smolarkiewicz, P.K., Wyszogradzki, A.A.: 2006, Dissipation in Implicit Turbulence Models: A Computational Study. J. App. Mech. - Trans. ASME 73, 469 – 473. doi:10.1115/1.2176749.
- Miesch, Brun, and Toomre, 2006 Miesch, M.S., Brun, A.S., Toomre, J.: 2006, Solar Differential Rotation Influenced by Latitudinal Entropy Variations in the Tachocline. Astrophys. J. 641, 618 – 625. doi:10.1086/499621.
- Miesch et al., 2000 Miesch, M.S., Elliott, J.R., Toomre, J., Clune, T.L., Glatzmaier, G.A., Gilman, P.A.: 2000, Three-dimensional Spherical Simulations of Solar Convection. I. Differential Rotation and Pattern Evolution Achieved with Laminar and Turbulent States. Astrophys. J. 532, 593 – 615. doi:10.1086/308555.
- Mitchell, 1916 Mitchell, W.M.: 1916, The history of the discovery of the solar spots. Pop. Astron. 24, 562.
- Passos, Charbonneau, and Beaudoin, 2012 Passos, D., Charbonneau, P., Beaudoin, P.: 2012, An Exploration of Non-kinematic Effects in Flux Transport Dynamos. Solar Phys. 279, 1 – 22. doi:10.1007/s11207-012-9971-2.
- Prusa, Smolarkiewicz, and Wyszogrodzki, 2008 Prusa, J.M., Smolarkiewicz, P.K., Wyszogrodzki, A.A.: 2008, EULAG, a Computational Model for Multiscale Flows. Comp. Fluids 37, 1193 – 1207. doi:10.1016/j.compfluid.2007.12.001.
- Racine et al., 2011 Racine, É., Charbonneau, P., Ghizaru, M., Bouchat, A., Smolarkiewicz, P.K.: 2011, On the Mode of Dynamo Action in a Global Large-eddy Simulation of Solar Convection. Astrophys. J. 735, 46. doi:10.1088/0004-637X/735/1/46.
- Schüssler, 1981 Schüssler, M.: 1981, The solar torsional oscillation and dynamo models of the solar cycle. Astron. Astrophys. 94, L17.
- Smolarkiewicz, 2006 Smolarkiewicz, P.K.: 2006, Multidimensional positive definite advection transport algorithm: An overview. Interna. J. Numerical Met. Fluids 50, 1123 – 1144. doi:10.1002/fld.1071.
- Smolarkiewicz and Charbonneau, 2012 Smolarkiewicz, P.K., Charbonneau, P.: 2012, EULAG, a Computational Model for Multiscale Flows: an MHD Extension. J. Comp. Phys., (submitted).
- Smolarkiewicz and Margolin, 1998 Smolarkiewicz, P.K., Margolin, L.G.: 1998, MPDATA: A Finite-Difference Solver for Geophysical Flows. J. Comp. Phys. 140, 459 – 480. doi:10.1006/jcph.1998.5901.
- Smolarkiewicz and Margolin, 2007 Smolarkiewicz, P.K., Margolin, L.G.: 2007, In: Grinstein, F.F., Margolin, L.G., Rider, W.J. (eds.) Implicit Large Eddy Simulation : Computing Turbulent Fluid Dynamics, Cambridge University Press, UPH, Shaftesbury Road, Cambridge, UK, 413 – 438. Chap. Studies in geophysics. doi:10.1017/CBO9780511618604.
- Smolarkiewicz and Prusa, 2002 Smolarkiewicz, P.K., Prusa, J.M.: 2002, In: Drikakis, D., Guertz, B.J. (eds.) Turbulent Flow Computation, Kluwer Academic Publishers, Dordrecht, 279 – 312. Chap. Forward-in-time differencing for fluids: Simulations of geophysical turbulence.
- Ulrich, 2010 Ulrich, R.K.: 2010, Solar Meridional Circulation from Doppler Shifts of the Fe I Line at 5250 Å as Measured by the 150-foot Solar Tower Telescope at the Mt. Wilson Observatory. Astrophys. J. 725, 658 – 669. doi:10.1088/0004-637X/725/1/658.
- Woodward and Colella, 1984 Woodward, P., Colella, P.: 1984, The numerical simulation of two-dimensional fluid flow with strong shocks. J. Comp. Phys. 54, 115 – 173. doi:10.1016/0021-9991(84)90142-6.
- Yoshimura, 1981 Yoshimura, H.: 1981, Solar cycle Lorentz force waves and the torsional oscillations of the sun. Astrophys. J. 247, 1102 – 1112. doi:10.1086/159120.