© ESO, 2017

1. Introduction

The solar magnetic cycle and the magnetohydrodynamic (MHD) dynamo that powers it are the drivers of space weather, the latter impacting technological assets ranging from power line networks to satellite communications. The solar dynamo is a MHD system in which the flow of plasma across a pre-existing magnetic field induces more magnetic field, offsetting Ohmic dissipation. This process, in the case of the Sun, also leads to cyclic polarity reversals of the large-scale magnetic component. The bulk of this field induction takes place in the Sun’s convection zone (CZ), where convective turbulence and differential rotation act as energy reservoirs and primary inductive flows. A detailed understanding of the dynamics of these flows, including the non-linear backreaction of the Lorentz force associated with the dynamo-generated magnetic field, is therefore crucial for understanding the dynamo process at a fully dynamical level.

The total velocity field can be decomposed into turbulent motions (from granulation to the larger convective scales) and large-scale, globally coherent motions organized on the scale of the Sun itself. The latter can be subdivided into (differential) rotation, fluid motion in the azimuthal direction, ê_φ, and meridional circulation (MC) motions in the {ê_r, ê_θ} plane.

Since the initial proposals of mean-field solar dynamo theory (Parker 1955; Steenbeck et al. 1966), differential rotation has been identified as a critical ingredient for magnetic field amplification, through what is known as the Ω effect. Thanks to the advances in helioseismology over the last decades, this flow has now been mapped for the interior of the Sun with a high degree of confidence for a wide range of latitudes and depths (Howe 2009). Analytical profiles that closely match the observed differential rotation profile are currently being used in several kinematic mean-field dynamo models (see review by Charbonneau 2010).

Increasing attention has been given to the meridional circulation following the development of surface magnetic flux transport models (Wang et al. 1989; Wang & Sheeley 1990) and flux transport (FT) dynamos (e.g., van Ballegooijen & Choudhuri 1988; Choudhuri et al. 1995; Dikpati & Charbonneau 1999). It also gained particular visibility as a key ingredient in the Babcock-Leighton (BL) dynamo framework (Babcock 1961; Leighton 1969; Wang & Sheeley 1991) since it carries the product of decayed active regions from low latitudes to the poles in these models. It is believed that this magnetic flux transport is what eventually reverses the polar field and switches the polarity of the solar magnetic dipole.

Observational data of the deep meridional flow is challenging to obtain because of its relatively low amplitude when compared to the omnipresent background of convective turbulence. This implies that the inclusion of the MC profile in dynamo models is prone to more uncertainties (when compared to the differential rotation) given the unavailability of a complete mapping throughout the convection zone. Until recently the MC was only known with confidence in the near surface layers where a mean poleward flow of around 15 ~ 20 m s^-1 is observed (Duvall 1979; Giles et al. 1997; Gizon et al. 2003; Ulrich 2010). In the past, this led mean-field modelers to adopt simple MC profiles (one cell per hemisphere) using extrapolations based on mass conservation in the solar interior (e.g., van Ballegooijen & Choudhuri 1988). Nowadays, improved helioseismic techniques allow the measurement of the MC between ±60° in latitude and down to 0.75 R_⊙. Nevertheless, these extended measurements are not yet conclusive as different groups obtain different radial profiles for the MC (cf. Zhao et al. 2013; Jackiewicz et al. 2015; Rajaguru & Antia 2015). One thing is almost certain: the meridional circulation profile is more complex than the commonly used one-cell-per-hemisphere configuration, with the possibility of having multiple cells stacked in radius and in latitude.

The first evidence of variations in the amplitude of the MC came from the correlation tracking measurements of surface features presented by Komm et al. (1993). Inspired by these early results, it was later proposed that coherent variations (on the cycle timescale) in the MC strength could directly account for the variations of the amplitude of the solar cycle (e.g., Nandy 2004; Passos & Lopes 2008; Lopes & Passos 2009; Karak 2010). Evidence that the MC actually exhibits this type of moderate variations came from the magnetic feature tracking measurements of Hathaway & Rightmire (2010) and Upton & Hathaway (2014). These authors showed that the surface poleward meridional flow varies in anti-phase with the solar cycle, decreasing its amplitude around solar maximum and speeding up again towards solar minimum. This type of solar-cycle related large-scale flow perturbation is also visible in the rotation through the torsional oscillations (Howard & Labonte 1980). This suggests that the source mechanism for torsional oscillations and MC cyclic variations has a magnetic origin. Recent observations of the equatorward drift velocity of the sunspot butterfly wings impose strong constraints on the variation of the meridional circulation at the storage depth of the toroidal flux tubes which, upon destabilization and buoyant rise across CZ, will emerge as active regions. According to Hathaway (2011) the drift velocity is independent of cycle strength, which means that either the MC is constant at the depth where the toroidal field accumulates or it might not be directly responsible for the equatorward migration of the activity belt.

The MC role in dynamo theory is that of a transport mechanism of magnetic flux between different regions. In BL flux-transport dynamo models, it transports magnetic flux towards the equator at the base of the CZ, and towards the poles at the surface. In models operating in the advection dominated regime, it also couples the bottom of the CZ where the strong radial shear is located with the photospheric layers where the Babcock-Leighton mechanism operates. However, the relevance of this coupling mechanism can be deemed less relevant if faster transport processes are considered (e.g., buoyancy, turbulent pumping, or high magnetic diffusivity). For this reason, and because of different parameterizations, the impact that MC amplitude variations have on the cycle strength varies between models operating in the advection or diffusion dominated regimes (Yeates et al. 2008; Lopes & Passos 2009; Nandy et al. 2011; Karak et al. 2014). Most of the time the authors loosely attribute these changes to the Lorentz force feedback. These models typically run in the kinematic regime, i.e., the background velocity field (differential rotation and MC) induces and organizes the magnetic field, but the latter does not feed back into the velocity field. A few exceptions in this respect are the non-kinematic mean-field simulations of Rempel (2006) and the low-order dynamo simulations of Passos et al. (2012), which both include magnetic feedback onto the MC.

An alternative to the Lorentz force feedback to explain why the MC varies, had its origins in the magnetic feature tracking measurements of Meunier (1999). The author identified variations in the surface meridional flow that can be associated with local inflows into active regions. Based on earlier ideas by H. C. Spruit and on these measurements, Cameron & Schüssler (2010, 2012) showed that, in the context of the surface flux transport (SFT) model framework, inflows into active regions could account for the global variations measured in the surface meridional flow. The cyclic surface MC variations are characterized by a weakening of the meridional flow on the poleward sides of the active region belt. This can be interpreted as the joint action of inflows towards the sunspot areas superimposed on a mean poleward meridional flow (Hathaway & Upton 2014). Surface flux transport simulations have also shown that variations in the amplitude of the MC can have an important impact on the amplitude of the next solar cycle because it determines the amount of polar field that is available for the next cycle’s production (Jiang et al. 2013; Upton & Hathaway 2014; Martin-Belda & Cameron 2016). The bottom line is that for both types of models, axisymmetric FT and SFT, the dynamics of the MC is a key element that needs to be better understood.

From a fluid dynamics point of view, the MC arises from angular momentum redistribution in the presence of rotation, convection and thermal gradients (e.g., Tassoul & Tassoul 1982; Rüdiger 1989; Tassoul 2000). Therefore, one of the most appropriate ways to study the MC is by means of global 3D MHD simulations. Such simulations allow researchers to disentangle the complex chain of interactions that the various physical processes create. Several authors have investigated the physical origin of differential rotation and MC profiles obtained in global simulations (for a detailed review, see Miesch 2005). From analysis based on the balances of the angular momentum and the thermal wind balance, there seems to be an agreement to the conjecture that the meridional motions appear as the consequence of the turbulent transport of angular momentum (e.g., Guerrero et al. 2013; Gastine et al. 2014; Featherstone & Miesch 2015, hereafter FM15). However, these studies were conducted in purely hydrodynamic simulations without considering the effects of magnetism.

Here we present an analysis that follows a similar methodology to that used in Brun et al. (2011), Miesch & Hindman (2011), and FM15, but for a 3D MHD global simulation that develops a large-scale magnetic field cycle (Passos & Charbonneau 2014, hereafter PC14). The MC profile obtained in this model was presented in Passos et al. (2015, hereafter PCM15) and exhibits interesting parallels to the helioseismic profile measured by Zhao et al. (2013). Another motivation for this work comes from the early realization that this meridional flow varies in intensity along the magnetic cycle, as noticed in Passos et al. (2012).

In this work we describe how the large-scale magnetic field interacts and modifies the MC profile in this simulation, and extrapolate our findings to the solar case. The paper is organized as follows. In Sects. 2 and 3 we present the 3D model used and describe the MC profile obtained. In Sect. 4 we conduct an angular momentum balance analysis and study the role of the large-scale magnetic field in inducing variations in the MC through gyroscopic pumping (GP). In Sect. 5 we develop an equation for the meridional force balance in the presence of magnetic fields and study the contribution from its various terms to the variations in the MC cell structure. We conclude in Sect. 6 with a discussion about the relevance of our findings to the actual case of the Sun.

2. The model

The results presented here are obtained through a dynamics analysis to the EULAG millennium simulation described in PC14, an Implicit Large-Eddy Simulation (ILES) of global solar convection produced with the EULAG-MHD code (Smolarkiewicz & Charbonneau 2013). The model solves the MHD extension of the anelastic equations of Lipps & Hemler (1982) in a spherical shell defined between 0.61 <r/R_⊙< 0.96. The convection is driven via a volumetric heating/cooling term in the energy equation. The domain is gravitationally stratified, rotates at the solar rate, and includes a convectively stable fluid layer underlying the convection zone. The governing equations read where D / Dt ≡ ∂/∂t + u·∇ is the convective (Lagrangian) derivative, u is the flow velocity, B is the magnetic field, Θ is the potential temperature, Ω = Ω₀(sinθ,cosθ,0) is the angular velocity with θ being the latitude, μ the magnetic permeability, ℋ the radiative diffusion, g the gravitational acceleration, and ρ_o the density stratification; α = 1 /τ defines the timescale of the Newtonian cooling term that drives convection and π′ is a density-normalized pressure perturbation that subsumes centrifugal forces and magnetic pressure. Prime quantities represent perturbations with respect to an arbitrarily selected ambient state (denoted by the subscript “e”). Quantities related to the basic state of the anelastic asymptotic expansion are denoted by the subscript “o”. A more detailed description of the model and its parameters can be found in Ghizaru et al. (2010); Racine et al. (2011); Smolarkiewicz & Charbonneau (2013); Cossette (2015); Guerrero et al. (2016a).

The numerical scheme used in the EULAG-MHD simulations does not consider any explicit dissipative terms. The MPDATA algorithm introduces a subgrid scale numerical viscosity at the minimal level required to maintain stability (Margolin et al. 2006). This effectively maximizes the Reynolds numbers and turbulence levels of the model for a given grid size. For the grid resolution used in this simulation, Strugarek et al. (2016) estimates an effective viscosity on the order of 1 × 10¹² cm² s^-1 for the convection zone. These authors show that even in an average sense this dissipation is strongly scale-dependent, varying by more than a factor of 10 between global and convective scales (see in particular their Fig. 5 and accompanying discussion). Using this information and the rms velocity in the bulk of the CZ we can say, as a conservative estimate, that this simulation is operating in a Reynolds number between 40 and 50. Strugarek et al. (2016) also show that in this simulation the Prandtl number is slightly higher than 1, which is consistent with the values assumed for ILES simulations. Taking an average value for the vorticity in the middle of the convection zone and the solar rotation rate, we can estimate a Rossby number of approximately 0.02. Comparison to equivalent dimensionless numbers computed from simulations using explicit dissipation is meaningful only up to a point. The MPDATA algorithm at the core of EULAG introduces dissipation in a spatiotemporally intermittent manner in response to the development of strong local gradient in advected variables.

The low viscosity regime allows us to reproduce a tachocline by placing a convectively stable layer at the bottom of the domain. The shear layer develops naturally and persists on a timescale much longer than the dynamo cycle period. Including the tachocline in the model allows some of the features of the solar magnetism to be reproduced, the most important being the generation of a deep-seated strong toroidal magnetic field which seems to govern the field dynamics in the entire domain. For instance, it drives the development of magneto-shear instabilities and the generation of non-axisymmetric turbulent modes in the stable layer that might influence the dynamo period (Lawson et al. 2015; Guerrero et al. 2016a). Furthermore, Guerrero et al. (2016b) argue that the magnetic tension due to the large-scale field at the tachocline induces the speedup and slowdown of the axial motions, ultimately establishing a torsional oscillations pattern (see also Beaudoin et al. 2013).

The millennium simulation was performed on a mesh of N_φ × N_θ × N_r = 128 × 64 × 47. This relatively low resolution allows for a long integration time, generating a solution extending over more than 1600 yr which spans about 41 polarity reversals (half magnetic cycle with a period of ~40 yr). The numerical experiments carried out by Strugarek et al. (2016) show that for this resolution the EULAG-MHD solutions agree with those obtained with the pseudo-spectral code ASH, which in turn has been accurately benchmarked against four other pseudo-spectral codes (Jones et al. 2011). While this simulation generates rotational torsional oscillations of solar-like amplitude (see Beaudoin et al. 2013), it does not reproduce the near-surface shear layer (Miesch & Hindman 2011; Hotta et al. 2015; Guerrero et al. 2016b), nor the inflows associated with active regions (e.g., Cameron & Schüssler 2012; Martin-Belda & Cameron 2016). Therefore, the relevance of these effects for the global meridional circulation is not addressed here.

3. Flows in the meridional plane

To highlight the influence that the magnetic field has on the dynamics of the large-scale meridional flows, we compare first the morphology of the MC of the millennium simulation (MHD) with a purely hydrodynamic (HD) analog simulation spanning nearly 330 yr for the same resolution and physical parameters. To ensure that the time intervals in which we will study the flow dynamics is not influenced by the initial conditions, we excluded the first 82 yr of the HD and 873 yr of the MHD simulated data from this analysis.

As usual, the MC radial and latitudinal components, respectively ⟨ u_r ⟩ and ⟨ u_θ ⟩, are obtained by zonally averaging these components of the velocity field. In Fig. 1 we present these quantities and the corresponding stream function averaged over a time interval of ~246 yr (for the HD simulation) and ~406 yr (for the MHD).

Fig. 1 Meridional plane displays of time averaged ⟨ ur ⟩, ⟨ uθ ⟩ and stream function. The top row is for a HD simulation, while the bottom row shows the same quantities for a MHD counterpart. For the HD (MHD) case, the temporal average corresponds to 246 (406) yr. In panels A and D, red/yellow (blue/green) denotes flows towards the north (south). In panels B and E, red (blue) denotes rising (sinking) flows. In panels C and F, we present the stream lines of the MC. The red (blue) tones denotes circulation in the counterclockwise (clockwise) direction. The black dashed line indicates the depth of the base of the convection zone. Below it we have the stable layers. The vertical dotted line in C) indicates the position of the tangent cylinder for reference. At the top layer, thin ticks mark 15° intervals. Panel C can be directly compared with Fig. 1g of FM15.

In both cases (HD and MHD), ⟨ u_r ⟩ and ⟨ u_θ ⟩ have a column-like morphology at low latitudes, with several thin cells parallel to the rotation axis (in agreement with the Taylor-Proudman constraint imposed by rotation). This behavior is typical of other similar global MHD simulations (see Miesch 2005, for examples). The small-scale pattern that appears near the top layers at low latitudes is an artifact from the low numerical viscosity of the code and from the time average of multiple asymmetric larger-scale flows over time.

Things become much more distinct between the two cases in the area contained inside a cylinder tangent to the base of the convection zone (in the equator) and parallel to the rotation axis (highlighted in Fig. 1C) by the vertical dotted line and the horizontal arrow).

For the HD case, inside the tangent cylinder (TC) we observe two circulation cells. One of them rotates counterclockwise (CCW) in the north, and spreads from the cylinder border to higher latitudes. This big cell is easier to identify in the southern hemisphere of Fig. 1C. There is another small cell, rotating clockwise (CW) in the north near the pole (see Figs. 1A–C).

Since outside the TC the flows (axis-oriented cells) have much higher velocities we purposefully used a low color saturation threshold in these figures to reveal the MC components everywhere. Otherwise, only the cells outside the tangent cylinder would be visible.

For the MHD case, on the other hand, Figs. 1D and E show a different scenario. We note that inside the TC the MC pattern is much more complex than in the HD case. In this region, the flow morphology shows two main cells of opposite circulation, with a prominent upflow between them at around ± 48° latitude in the bottom half of the CZ. Near the poles we have several smaller cellular structures (see Fig. 1F).

The morphology of the MC evolves as the magnetic cycle progresses and these changes are significant especially inside the TC. This region coincides with the location where the large-scale magnetic cycle develops and where other solar-like dynamic phenomena such as torsional oscillations are observed (Beaudoin et al. 2013). In order to investigate how the MC evolves with the magnetic cycle, we consider a time interval of ten sunspot-like cycles (half magnetic cycles), spanning from cycle 21 to cycle 30 of the millennium simulation and covering a period of 406 yr (see Fig. 2).

Fig. 2 Magnetic cycle proxy (for the northern hemisphere). The vertical gray areas represent the maxima and minima phases considered. For reference we labeled these phases for simulated cycle number 24.

The proxy for the cycles presented here is the same as defined in PC14, i.e., the normalized amplitude of ⟨ B_φ ⟩ integrated over an extended region (in r and θ), centered around 48° near the base of the CZ, where this field component accumulates (here just for the northern hemisphere). In order to emphasize the differences along the cycle, we compare quantities averaged out around times of minimum and maximum. These two epochs are labeled for cycle 24 in Fig. 2. We generically define the minima and maxima phases as intervals of ± 3 yr around the time of cycle minimum and maximum respectively. Figure 3 presents snapshots of the magnetic field components and MC morphology taken at these two phases of the simulated cycle 25 (starting around t = 1045 yr).

Fig. 3 Zonally averaged magnetic field components sampled at the minimum (top row) and maximum (bottom row) of cycle 25. Panels A and D show ⟨ Bφ ⟩ with the MC stream function overplotted, where the dashed (solid) contour lines represent cells circulating in the counterclockwise (clockwise) direction. Panels B and E show ⟨ Bθ ⟩, and panels C and F show ⟨ Br ⟩.

Outside the TC, the behavior of the MC seems to be driven mainly by the same mechanism as in the HD case since it still maintains the well-defined column patterns. In what follows we focus on the region inside the TC.

At maximum (Fig. 3D we have a cell rotating CCW between 48° and the pole in the bottom half of the CZ, and another smaller cell rotating CW between 70° and the pole in the upper half of the CZ. Between 25° and 46° (still inside the TC area), there is a CW rotating cell that spreads through most depths (and with a tendency to become oriented almost parallel to the rotation axis). These mid-latitude cells share an upflow section at around 48°, the latitude at which the zonally averaged toroidal field achieves stronger values. During the minimum phase the two mentioned cells above 48° get mixed and almost vanish, while the lower latitude cell also decreases in strength. Studies by Passos et al. (2012, 2016) indicate that the magnetic field has an important role in modulating the horizontal component of the MC within the CZ (specially inside the TC region). In the next section, we explore the nature of the interaction between magnetic field and flows that lead to cyclic variations in the MC.

4. Mechanisms of MC spatiotemporal variations

The acceleration of meridional flows is described by the zonal component of the vorticity equation, which we consider in Sect. 5. A variety of global convection simulations, mean-field models, and theoretical arguments consistently suggest that the dominant terms in this equation are the Coriolis and baroclinic terms, which give rise to thermal wind balance (TWB); see Kitchatinov & Rüdiger (1995), Durney (1999), Brun & Toomre (2002), Miesch et al. (2006), Balbus et al. (2009), Miesch & Toomre (2009), Gastine et al. (2014), Warnecke et al. (2016), Kitchatinov (2016), and references therein.

However, the TWB equation itself is degenerate in the meridional flow. It cannot be solved to yield the steady-state MC profile. Global circulations can be driven by departures from TWB, which are thought to occur in the upper and lower boundary layers of the convection zone, due to turbulent stresses (Balbus & Latter 2010; Miesch & Hindman 2011; Warnecke et al. 2016). This is commonly found in mean-field models where the meridional momentum transport by the convective Reynolds stress is modeled as a turbulent diffusion (Dikpati 2014; Kitchatinov 2016). Here the poleward circulation in the solar surface layers is attributed to the axial gradient of the angular velocity, ∂Ω /∂z, which appears in the zonal vorticity equation through the Coriolis force (see Sect. 5).

A related possibility is the phenomenon of gyroscopic pumping (GP). This occurs when a steady source of zonal momentum (an axial torque) induces a meridional flow by means of inertia. This operates by means of the Coriolis term in the zonal vorticity equation (∝∂Ω /∂z) as noted in the previous paragraph, but it does not require a sustained departure from TWB. In a steady state, the amplitude and structure of the MC is largely determined by the nature of the axial torque, as described in detail by Miesch & Hindman (2011). Gyroscopic pumping (GP) has been definitively demonstrated by Haynes et al. (1991) and is thought to play an important role in the exchange of air masses between the stratosphere and troposphere (Holton et al. 1995). It has been reproduced in the laboratory in the classic Plumb-McEwan experiment (see discussion by McIntyre 1998). GP has also been invoked to account for the circulation in planetary atmospheres (Read 1986) and stellar radiative zones (Spiegel & Zahn 1992; Fritts et al. 1998; Garaud & Bodenheimer 2010; Wood et al. 2011).

The maintenance of meridional flows by GP has been found in previous mean-field solar convection models by Rempel (2005) and in previous global convection simulations by Brun et al. (2011), FM15, Gastine et al. (2014), Hotta et al. (2015), and Guerrero et al. (2016b). In non-magnetic models, the source of the zonal torque is the angular momentum transport by the Reynolds stress. However, magnetic torques can also induce a meridional flow in an analogous way. In this section we investigate this mechanism in detail for our cyclic convective dynamo.

The angular velocity Ω in our model is defined as (6)where Ω₀ = 2.42405 × 10^-6 s^-1 is the angular velocity of the coordinate system, and λ = rcos(θ) is the momentum arm with θ being the latitude. The differences between the HD and MHD differential rotation profiles are shown in detail by Beaudoin et al. (2013, see their Fig. 2).

The redistribution of specific angular momentum, ℒ = λ²Ω, plays a central role in the establishment and maintenance of mean flows. Angular momentum transport by the convective Reynolds stress not only governs the magnitude of the differential rotation ΔΩ, but it also regulates the structure and amplitude of the meridional circulation by means of GP (Miesch & Hindman 2011). In a steady state, the meridional acceleration induced by the inertia of the differential rotation is offset mainly by horizontal pressure gradients. This can be expressed as a balance between the Coriolis and baroclinic terms in the zonal vorticity equation and is known as the thermal wind balance (TWB). Any torque that disrupts this balance through a local acceleration or deceleration of Ω will induce a meridional flow that will act to restore the equilibrium profile of Ω that is consistent with TWB.

Thus, there are several ways in which magnetism can influence the MC. The first is through the direct acceleration of the meridional flow due to the mean¹ meridional Lorentz force. The second is by exerting a torque through the zonal component of the mean Lorentz force that alters the rotation profile, Ω. This is the mechanism of GP. A third way for magnetism to influence the MC is by altering the convective momentum and energy transport by means of the non-axisymmetric components of the Lorentz force, namely the Maxwell stress. In this section we will demonstrate, as in the non-magnetic convection simulations of FM15, that it is the second mechanism, GP, that largely accounts for the structure and variability of the meridional flow that we see (Sect. 3).

The equation that describes the conservation of angular momentum in an anelastic system (7) (see Appendix of Miesch & Hindman 2011) gives us some information about the physical mechanisms involved. This equation is obtained by multiplying the zonal component of the momentum equation by λ, and then averaging over longitude (indicated by angular brackets, ⟨⟩). For an inviscid simulation like ours (neglecting numerical diffusion), this yields (7)where , and the terms on the right-hand side include angular momentum fluxes due to the Reynolds stresses, Maxwell stresses, and the large-scale magnetic fields (magnetic tension). These fluxes are defined as where we use the classical Reynolds decomposition of the field and flow components in zonal means (associated with large-scale flows) and fluctuations (associated with small-scale turbulence), e.g., . Equation (7) tells us how angular momentum variations due to local zonal forcings (net axial torque ℱ) can induce variations in the meridional flow. The same approach was used by Miesch & Hindman (2011) to study how the MC is established in the solar near surface shear layer. It was also used by Brun et al. (2011) and FM15 to study how the MC is established in hydrodynamic (non-magnetic) simulations of global convection (with and without a tachocline). These authors found that the Reynolds stress term has a preponderant role in establishing the overall amplitude and morphology of the MC. More recently, Guerrero et al. (2016b) used Eq. (7) to identify the main agent driving torsional oscillations in a global dynamo simulation. Their results suggest that the magnetic tension at the bottom of the convection zone induces axial torques that periodically speedup and slowdown the angular velocity.

4.1. Angular momentum balance in the HD case

For comparison purposes we begin by performing this analysis to our HD simulation. Since there are no magnetic fields involved, Eq. (7) reduces to (11)Figure 4 shows the individual terms of Eq. (11) averaged over a time interval of 246 yr. When ⟨ ρ₀ ⟩ ∂ℒ /∂t is smaller than the other terms, we can assume that ℱ is due to the sole action of the Reynolds stresses. When ℱ > 0 (red lines and shades) the net torque is prograde inducing a meridional flow away from the rotation axis. While ℱ < 0 (blue lines and shades), the net torque is retrograde and induces a flow towards the rotation axis.

To interpret this figure we first focus on the Reynolds stress shown in Fig. 4C. Outside the TC, below 45°, it exhibits a divergence in the mid-lower CZ (blue) and a convergence in the upper CZ (red), with contours approximately aligned with the rotation axis. This is consistent with the work of FM15 and signifies the transport of angular momentum by sheared banana cells².

Fig. 4 Panel A: ρ0∂ℒ /∂t, B) ⟨ ρ0um ⟩ ·∇ℒ and C) in kg m-1 s-2, averaged over 246 yr. Panels B and C have the same scale (105) but A) is two orders of magnitude smaller (103), indicating that a steady state has been achieved.

As described by FM15, the divergence of F^RS in the lower CZ near the equator induces a prominent CW circulation cell in the northern hemisphere (NH), immediately outside the TC. This effect is also observed in our simulation (Fig. 1C).

In FM15 (and in the MHD simulation considered here (Fig. 1F), this CW (blue) cell lies between two CCW (red) cells above and below. However, in the HD simulation considered here, the upper cell at low latitudes (outside the TC) is absent, with the circulation in this region dominated by a series of smaller-scale cells. This difference can be most likely attributed to different effective viscosities between the two models.

Inside the TC, the CCW cell that pervades most of the CZ in the NH is somewhat less evident here, but it is present as well. In Fig. 1C is easier to identify its antisymmetric counterpart in the southern hemisphere (SH). The different morphologies in this CCW cell found in FM15 and our HD simulation can be partially attributed to the presence of the stable zone and overshoot region, absent in FM15. The inward angular momentum transport by downflow plumes leads to a convergence of the angular momentum flux in the overshoot region that persists to very low latitudes, establishing a strong equatorward flow. The high density in the overshoot region makes this a substantial contribution to the mass flux and mass conservation largely accounts for the poleward flow in the mid CZ. This establishes the strong CCW circulation cells near the base of the CZ. Similar results were seen in the penetrative convection simulations by Miesch et al. (2000, see their Fig. 16a).

The close correspondence between panels B and C of Fig. 4 and the small amplitude of A indicate a statistically steady state. Thus, the advection of angular momentum by the MC balances the transport of angular momentum by the convective Reynolds stress. Furthermore, readjustments in ℒ occur on the timescale of a several days, which means that any averaging over longer periods will result only in a small residual. The small differences between frames B and C can be attributed to the contribution of numerical viscous fluxes (especially in the θ direction, see Guerrero et al. 2016b), to a residual ∂ℒ /∂t, to the finite duration of the temporal averaging, and to the different numerical methodologies used while running the model and the a posteriori analysis. We elaborate on this in the next section.

4.2. Angular momentum balance in the MHD case

Next, we apply the same analysis procedure to the MHD simulation where a dynamo generated large-scale magnetic cycle contributes to the transport of angular momentum. In Fig. 5 we show a comparison between the left-hand side and the right-hand side of Eq. (7) averaged over the ten cycles. Here again, the first term of Eq. (7) shows only a small residual time dependence two orders of magnitude smaller than the terms on the right-hand side. As above, the numerical diffusion can account for most of the small differences between panels 5B and C. The three components that make up the right-hand side in Fig. 5C are shown individually in Fig. 6.

Fig. 5 Panels A and B: same terms as in Fig. 4 for the MHD simulation also in kg m-1 s-2. Panel C: entire right-hand side of Eq. (7), which include the contributions from FRS, FMS, and FMT. All quantities are averaged over ten cycles (406 yr).

Fig. 6 Components of the right-hand side of Eq. (10) and Fig. 5C, plotted individually: A) , B) , and C) .

The influence of magnetism on the net axial torque can be seen by comparing Figs. 4C and 5C. The general pattern outside the TC is similar in the two figures, with a region of divergence (blue) in the lower CZ at low latitudes, straddled above and below by regions of convergence (red). This is somewhat expected since the strong columnar flow behavior in this region is mainly maintained by Reynolds stresses, and the large-scale magnetic field is weak. The Reynolds stresses induce a series of three stacked circulation cells at low latitudes. The small differences with the HD case in this region can be attributed mainly to the contribution of Maxwell stresses; see Figs. 1C and F. As in the HD case, the angular momentum transport by these circulation cells largely balances the net torque; compare Figs. 5B and C. Still outside the TC, but in the upper part of the CZ, around 35° there are some small differences that can be attributed to local magnetic fields. We address this further below.

As noted in Sect. 3, the main difference between the HD and MHD cases is the paired set of MC cells inside the TC at mid-latitudes, seen for example in Fig. 1F in conjunction with an upwelling in the lower CZ at a latitude of about 48°. We attribute these circulation cells to GP induced by the magnetic part of the axial torques seen in Fig. 5C. At the base of the CZ, these torques are accelerating the rotation rate at latitudes higher than 48° (red), and decelerating lower latitudes (blue), inducing a convergent flow which produces the mid-latitude upwelling. At slightly larger radii in the lower CZ the pattern of torques reverses, with deceleration (blue) and acceleration (red) at latitudes poleward and equatorward of 48°. This establishes the horizontally diverging flow in the mid CZ that closes off the pair of mid-latitude circulation cells. Though such a closure is not required by GP (closed circulation cells are always ensured by mass conservation), the quadrupolar pattern of red-blue-red-blue serves to enhance the mid-latitude circulation cells and to keep them localized in the lower CZ (see Fig. 3D). This interpretation is confirmed by studying the angular momentum transport by the MC in Fig. 5B, which shows a similar quadrupolar pattern at mid-latitudes.

It is clear from Fig. 6 that this quadrupolar pattern in the net torque arises from the large-scale Lorentz force (panel C). However, in order to interpret this, we must begin with the Reynolds stress component in panel A. As noted above, F^RS is dominated by banana cells, which transport angular momentum cylindrically outward (away from the rotation axis) and equatorward. The equatorward component of the transport leads to a divergence of F^RS at mid-latitudes in the upper CZ (blue), i.e., the Reynolds stress is extracting angular momentum from the top half of the CZ at latitudes 30°–55° in order to establish the solar-like differential rotation. At the same time, in the bottom half of the CZ, angular momentum is also being transported downwards. At mid-latitudes, downflow plumes carry angular momentum to the stable zone and their deceleration in the overshoot region gives rise to a convergence of the angular momentum flux. This is visible in Fig. 6A as a red stripe near the base of the CZ at latitudes between ± 60°, which becomes particularly prominent between latitudes of ± 25°–50°. This acts to accelerate the rotation rate near the base of the CZ and decelerate the rotation rate in the upper CZ. The magnetic torques respond to this Reynolds stress. In particular, the quadrupolar pattern of F^MT at mid-latitudes arises when the torques exerted by the Reynolds stress are extended poleward by magnetic tension. For example, in the upper CZ at a latitude of 40°, the Reynolds stress is acting to decelerate Ω (blue). Magnetic tension opposes this local deceleration (red) and spreads it to higher latitudes (blue). Similarly, when the Reynolds stress acts to accelerate the fluid near the base of the convection zone, the rigidity imparted by magnetic tension serves to “drag” higher latitudes along. This acts to decelerate mid-latitudes (blue) and accelerate higher latitudes (red).The Maxwell stress also opposes the Reynolds stress, both in the overshoot region and in the upper CZ (Fig. 6B). However, in contrast to the large-scale magnetic tension, F^MS is more diffusive in nature, and thus more localized.

The quadrupolar mid-latitude pattern of angular momentum divergence and convergence bears an interesting resemblance to the angular momentum cycle described by Gilman et al. (1989; see their Fig. 1a). They postulated a poleward angular momentum transport by some unspecified process in the solar tachocline that offset the equatorward transport by convective Reynolds stresses in the CZ. They even identified magnetic stresses as a possible candidate process. Poleward angular momentum transport by magnetic stresses in the tachocline is an important component of several tachocline confinement models (reviewed by Miesch 2005). Our global MHD convection simulation clearly demonstrates this. Though the magnetic stresses in our model vary over the course of a magnetic cycle (see below), they do induce a net angular momentum transport toward the poles near the base of the CZ, as seen in Fig. 6C.

Another minor difference in the axial torques distribution between the HD and MHD cases is that the latter has a more prominent prograde (positive) torque in the upper CZ at the equator (compare Figs. 4C and 5C). This is reflected in the MC profile of the MHD model where an upflow seems to be induced (see in Fig. 1F the paired set of circulation cells in the upper CZ near the equator, red in the NH, blue in the SH). This does not appear to be due directly to the Lorentz force, but rather to a modification of the Reynolds stress by magnetism (Fig. 6A). This may be attributed to the diffusive nature of the Maxwell stress, which tends to make the flow more laminar, enhancing the cylindrically outward angular momentum transport by banana cells (Brun et al. 2004; Nelson et al. 2013; Fan & Fang 2014; Karak et al. 2015; Hotta et al. 2016).

Furthermore, in the upper CZ near a latitude of about 28°, the MHD case exhibits a torque pattern that is not present in the HD case (compare Figs. 4C and 5C). This is due to the presence of a secondary dynamo that exists in that region. A complete characterization of this secondary (weaker and short period) dynamo mode is presented in Beaudoin et al. (2016). The authors study its interaction with the primary mode and show how it is maintained by a latitudinal shear. Although the magnetic field created in this region its weaker than that produced by the main dynamo mode that operates near the bottom of the convection zone around 50°, it still contributes to MC variability.

The dominant role of the large-scale magnetic torque in regulating the MC inside the TC suggests that the cyclic variation of the large-scale field should account for the cycle dependence of the MC, as discussed in Sect. 3. We find that this is indeed the case, as demonstrated in Fig. 7. The most apparent difference between cycle minimum and cycle maximum is in the mean magnetic torque, represented in panels D and H. Although the transport of angular momentum by the large-scale Lorentz force relies on the zonal poloidal field, the mid-latitude poloidal fields in the lower CZ are generally strongest when the toroidal bands are strongest. Thus, the establishment of meridional flows by magnetic torques via GP is most efficient at cycle maximum. This is when the mid-latitude upwelling at ± 48° in the lower CZ is strongest. Signatures from the secondary dynamo mode operating in the upper CZ near ± 28° latitude are also apparent in panels C, D, G, and F of Fig. 7.

As expected, F^RS shows little difference between cycle minima and maxima (Figs. 7B and F). The Maxwell stress (Figs. 7C and G) acts as a diffusive component (opposing the Reynolds stress below 48° in the CZ) with a time-varying component (opposing the large scale magnetic torque in the stable layers and above 48° in the lower CZ). This likely reflects non-axisymmetric structure in the bands, as opposed to turbulent diffusion by smaller-scale, more chaotic motions. The transport of angular momentum by the MC depicted by panels A and E of Fig. 7 shows the spatial combination of the action of the right-hand side torques.

Fig. 7 Several terms of Eq. (7) averaged over ten cycle minima (A–D) and 10 cycle maxima (E–H). We focus only on the NH for simplification. Panels A and D: left-hand side of Eq. (7); B and F: ; C and G: ; and D and H: . All the contours have the same color scale as in Fig. 6. The vertical dotted lines indicate the TC for reference.

In Fig. 8 we take a closer look at the time evolution of the various torque components and the response of the MC. Panels 8A and B show the net axial torque overlaid with the direction and amplitude of the meridional flow, shown as arrows. This confirms our interpretation in terms of GP; regions of negative torque (flux divergence: blue) generally exhibit a flow component towards the rotation axis (in addition to an axial flow component sustained by mass conservation) and regions of positive torque (flux convergence: red) generally exhibit a flow component away from the rotation axis. In panel 8B, it is clear that the mid-latitude upwelling in the lower CZ during cycle max, and the associated torques, are induced by the presence from the toroidal bands (white contours).

Fig. 8 Panels A and B: right-hand side of Eq. (7) in the NH averaged over cycle minima and cycle maxima. The vector field represents the direction of the meridional flow and the vertical dotted line represents the TC. The white contour lines in B) show the area where the toroidal field accumulates at cycle maximum. The bottom panels show the time evolution of the left-hand side of Eq. (7) and the individual terms on the right-hand side (as labeled), sampled at the numbered locations. The data in the lower three panels were smoothed with a one-year average filter.

The bottom panels of Fig. 8 show the evolution of the left-hand side of Eq. (7) together with the individual components of the net axial torque sampled at three different locations:1 ◯the area of higher correlation between the amplitudes of the toroidal field and the horizontal flow,2 ◯where the secondary (weaker) dynamo mode is operating and3 ◯an area outside the TC and away of the main influence of magnetic torques.

The variability in region1 ◯is clearly driven by the large-scale Lorentz force, F^MT (orange line), as discussed above. In the heart of the toroidal bands, magnetic tension is accelerating the rotation rate by extracting angular momentum from lower latitudes. Maxwell stresses resist this acceleration (blue line). The sum of these two components (not shown here) closely matches the black line that represents the left-hand side of Eq. (7). The Reynolds stress is almost negligible here, showing signs of cycle modulation but at 1 to 2 orders of magnitude below the amplitude of the other signals. The small amplitude of the Reynolds stress reflects the location where the toroidal bands accumulate: inside the TC and close to the base of the CZ where the convective amplitude is weak. In Fig. 9A we present the individual contributions of the two terms of the left-hand side of Eq. (7) and the joint contribution of the torques of the right-hand side sampled at region1 ◯for an interval of 5 yr taken around a cycle maximum³. This figure clearly shows that the angular momentum response due to variations in the zonal flows happens at a timescale on the order of one month (black line) and that the long-term variations in the MC (blue line) follow the GP forcing exerted by the right-hand side torques.

Fig. 9 Panel A: individual terms of the left-hand side of Eq. (7) (in black and blue) and the total contribution of the right-hand side (orange) sampled at region1 ◯for a five-year interval, taken at cycle maximum. Panel B: individual terms of the left-hand side (black and purple) and right-hand side (red) of the meridional force balance (Eq. (13) in Sect. 5) for the same location and for the same time interval.

In terms of dynamics, near the bottom of the CZ, the acceleration of the fluid within the heart of the toroidal bands induces an equatorward flow, while the extraction of angular momentum from lower latitudes induces a poleward flow. This establishes a meridional flow that converges horizontally into the magnetic toroidal bands and then turns upwards (see Figs. 8B and 10). The location of this deep convergence and upwelling is shifted slightly towards the equatorward edge of the bands, as shown in Fig. 10B. It is also interesting to notice from panel 10C that for latitudes higher then 50° (above region1 ◯) we can observe an apparent migration of equatorward flows (blue in the north) from the middle of the CZ to the upper layers as the cycle unfolds (marked with an arrow). This is not actually a migration. It is a decrease in the radius of the near surface CCW MC cell that can be found at those latitudes (see Fig. 3) caused by the decrease in the torque near the surface. As the magnetic torque becomes weaker from the cycle maximum to minimum, the near surface CCW cell becomes thinner and its equatorward flow (blue) approaches the top layers. At the same time there is a CW MC cell near the pole (see Fig. 3) that gradually expands to lower latitudes. The surface equatorward section of this higher latitude cell associated with the previous behavior is what finally establishes the observed superficial dynamical pattern above 50°.

Fig. 10 Mean latitudinal flow ⟨ uθ ⟩ as a function of latitude and time at A) the top layers (r = 0.95 R⊙) and B) near the base of the CZ (r = 0.72 R⊙). For a better contrast, in panel B the displayed quantity is actually 2⟨ uθ ⟩. The black contour lines represents ⟨ Bφ ⟩ (solid and dashed for positive and negative polarities respectively). The horizontal solid line represents the latitude at which panel C is sampled. The time interval covers cycles 23 to 28. The color scale saturates at ± 3 m s-1 (red northward, blue southward). Panel C is a radius vs. latitude plot of ⟨ uθ ⟩ taken at 65° north. Both ⟨ uθ ⟩ and ⟨ Bφ ⟩ are smoothed over 6 months. The vertical dotted lines represent the maximum and minimum of cycle 25. The dashed line marks the tachocline depth, the dotted line below it is the depth where ⟨ Bφ ⟩ contours are computed and the solid line is the depth where panel B is sampled.

It is clear from Figs. 8–10 that the accumulation of strong toroidal fields induces cyclic variations in the meridional flow. Though MC variability induced by Lorentz force feedbacks has been studied within the context of non-kinematic mean field models (e.g., Rempel 2006), it is worth noting that this line of causality is in stark contrast to the kinematic assumption that is often adopted in mean-field solar dynamo modeling.

The time evolution of the axial torques and meridional flow in region2 ◯is more erratic (see Fig. 8). In this small border-like region the large-scale magnetic torque and the Maxwell stress have the same sign. The meridional flow compensates this changes in the axial torque and a small CCW circulation cell appears during the maximum. In this region, the magnetic torque peaks (at negative values) around the time when the deep toroidal field in region1 ◯starts to rise. Meanwhile, there is also a pronounced shorter-term variability that is associated with the periodicity of the secondary dynamo mode.

Finally, we also highlight region3 ◯in Fig. 8. Here, where the large-scale magnetic field has a small influence, we observe almost no variation with the magnetic cycle. As in the HD case, the angular momentum transport is dominated by banana cells (F^RS).

5. Meridional force balance in the presence of large-scale magnetic fields

In addition to the GP mechanism presented in the previous section, the MC might also be influenced by the presence of entropy gradients throughout the CZ. These thermal gradients imply that surfaces of constant mean pressure and density do not overlap completely, which gives rise to a baroclinicity contribution in the vorticity equation. This contribution influences the whole system causing the MC to readjust in order to achieve the TWB. In the context of purely HD simulations, e.g. Brun & Toomre (2002); Miesch et al. (2006); Brun et al. (2010, 2011) have shown that the action of baroclinicity in TWB, together with the Reynolds stresses are the key factors for the system to achieve equilibrium. In addition, by studying fully MHD models of solar-like stars, Varela et al. (2016) unveiled the relevant role of a large-scale magnetic field in the TWB influence in differential rotation. In this section, we assess the influence of the magnetic field on the MC by comparing the balance conditions in HD and MHD models.

We find it useful to start by defining an equation for the evolution of the vorticity, ω = ∇ × u. Applying ∇ × to the momentum Eq. (1) yields (12)where ω_a = (∇ × u) + 2Ω₀ is the absolute vorticity.

As mentioned before, TWB studies have been carried out mainly for HD simulations and generalized under the argument that the magnetic field influence can be neglected. However, as we show in the previous section, the magnetic field has an important role in the GP forcing mechanism. Therefore, in order to gauge how important this magnetic contribution is for TWB, we develop an equation for the meridional force balance (MFB) by computing the zonally averaged component of the vorticity evolution Eq. (12). A similar equation is also presented Strugarek et al. (2011) and in the recent work of Varela et al. (2016)⁴ but it is used in a different context. Further details are presented in the appendix. The MFB equation assumes the form (13)where we can consider a stationary state (system in equilibrium) for ∂ ⟨ ω_φ ⟩ /∂t = 0. The second term on the left-hand side is the ê_φ component of (2Ω₀·∇)u, more commonly represented by 2Ω₀∂ ⟨ u_φ ⟩ /∂z, where ∂/∂z is the derivative in the direction parallel to the rotation axis. We follow Brun et al. (2011) for the naming for the several terms in Eq. (13). The main difference between the two MFB Eqs. (11) is that their baroclinic term is written in terms of the entropy while ours is written in terms of potential temperature; the viscous term is absent in our case because our simulation has no explicit viscosity and we take into account the contribution of the magnetic field. We would like to note that ∂ω_φ/∂t is not strictly zero but it oscillates around zero on a timescale commensurate with rotation, as shown in Fig. 9B. As expected, this term varies with the same frequency of the accelerations and decelerations of the zonal flows depicted by ∂ℒ /∂t. Therefore, using the previous arguments, we can consider that this system is in a quasi-static equilibrium.

5.1. MFB analysis for the HD case

All the quantities needed to compute Eq. (13) can be extracted directly from our numerical simulations. For comparison purposes we start by showing the meridional force balance computed for the HD case. In Fig. 11 we show meridional plots for the left- and right-hand sides (and the first four individual terms) of Eq. (13) averaged over the 246 yr interval and for the NH. For simplicity we restrict our discussion to the NH and divided all terms by 2Ω₀.

Fig. 11 Panels A and D: left- and right-hand sides of the MFB equation for the HD simulation. The terms were computed for the NH only, and with color saturation at ± 3 × 107 s-1. The vertical dotted line represents the TC. In the other four panels we present the individual contributions on the right-hand side: B) stretching, C) advection, E) compressibility and F) baroclinicity (with color saturation at ± 2 × 10-7 s-1).

The good agreement between the left- and right-hand sides of (13) represented by panels A and D of Fig. 11 is a clear indication that the system is actually very close to stationarity. Although panel D is the sum of the individual contributions of the right-hand side of Eq. (13), we can see that for most of the CZ the baroclinic term (panel F) is dominant. It has a large positive contribution in regions of strong rotational shear, such as the base of the CZ and at low latitudes outside the TC. This is in agreement with the results of Brun et al. (2011) who associate this behavior with the presence of strong thermal gradients. The baroclinic term is only weakly negative in the inner border of the TC, a region where rotation has almost no spatial variation. The contributions of stretching and advection are only relevant near the pole, while vorticity compressibility is negligible.

5.2. MFB analysis for the MHD case

The same analysis is now applied to the MHD simulation, but this time taking into consideration the complete form of Eq. (13). We use the same time interval of ten cycles as in the previous section for the averaging. In Fig. 12 we opted not to show the compressibility contribution because, as in the HD case, it is negligible when compared to the other terms.

Fig. 12 Panels A and D: left- and right-hand sides of the MFB equation in the NH for the MHD case. The individual terms of the right-hand side are shown in panels B, C, E, and F. The magnetic contribution is the sum of the two individual magnetic terms. The vertical dotted line represents the TC. The color scale saturates at ± 1.5 × 10-7 s-1.

Outside the TC, the MHD simulation behaves similarly to its HD counterpart, with baroclinicity being the dominant contribution. In terms of balance between left- and right-hand sides of (13), panels A and D show a very good agreement in most of the CZ. Differences are only evident inside the TC in the top layers and slightly near the poles. We attribute these differences to the upper and polar boundary conditions for the magnetic field. The difference observed between panels 12A and D in the top boundary might be due to the radial magnetic field we enforce at the surface. This means that when we have a strong poloidal field located near the surface, it will be forced over a couple grid points by the boundary condition into the radial direction. This is exactly the case of the poloidal field configuration during cycle maximum (see Fig. 3E). This problem might be alleviated in future simulations by introducing more realistic top boundary conditions like that used in Warnecke et al. (2016). The differences between the left- and right-hand sides during cycle minimum (not shown here) are much smaller. It is the radial derivative of the poloidal field present in the second magnetic term of Eq. (13) that is responsible for this “artificial” contribution. There are two other possible sources of error that can explain the minute differences we find in this balance calculation (and in the previous section as well). The first is numerical diffusivity, which we cannot measure directly. The other issue is related with the different numerical methods used to compute derivatives and other composite quantities in the main code during the simulation and a posteriori. During the simulation run, EULAG numerics computes central cell values and fluxes across the cell borders, while the type of analysis that we perform a posteriori assumes values computed in the cell corners using centered finite differences. Differentiation across the poles can also introduce some artifacts. Nevertheless, the very good match obtained in the HD case indicates that these two sources of error are in fact very small, and that the main issue here seems related to the magnetic field boundary conditions. Despite these possible sources of uncertainty, we consider that there is a general good agreement between left- and right-hand sides for most of the CZ.

By comparing Figs. 12D–F, we can see that inside the TC the baroclinic term also has an important role in establishing MFB, especially in the area close to the inner TC border (left side of the vertical dotted line in the panels). In the remainder of the CZ (roughly above 40°), the MFB is maintained by a combination of baroclinic and magnetic contributions (especially from the mean magnetic field). Miesch et al. (2006) argued that MFB (baroclinicity mainly) should have important effects in the tachocline and lower CZ except in regions where this equilibrium can be disrupted by strong magnetic fields. Our analysis also points in that direction. The stretching and advection terms are almost completely overshadowed by the contributions of the two previous terms.

Fig. 13 Panels A to C: left-hand side, baroclinic, and magnetic contributions of the MFB averaged over 10 cycle minima. Panels D to F: same quantities as the top row averaged over ten cycle maxima. The color scale saturates at ± 1.5 × 10-7 s-1.

We can highlight the influence of the magnetic field by comparing the MFB at times of cycle minima and maxima. In Fig. 13 we present an average of the left-hand side, baroclinic, and magnetic contributions over ten cycle minima (top row) and maxima (bottom row).

Inside the TC, panels 13A and D show that for latitudes between 0° and 41° (approximately where the TC intersects the top boundary) the variation between minima and maxima for the left-hand side is small and almost restricted to the stable layers. This is also applicable for the other two terms and is somewhat expected because the magnetic field does not have a significant presence at low latitudes in the bulk of the CZ. Inside the TC, during cycle maxima there is an enhanced negative (purple) baroclinic region that spreads from ~48° to 85° in latitude (panel E). During minimum the baroclinicity tends to “relax” to a profile closer to the HD case, with a positive enhancement in the shear region at the bottom of the CZ at high latitudes and in most of the high latitude range of the bulk of the CZ (panel B). The “quadrupolar cell pattern” that we find around the interface between stable and convective layers in the magnetic term, is similar (with opposite sign) to that found in the baroclinic term during both maxima and minima, indicating that these two quantities tend to balance each other.

The existence of this pattern in the stable layers in both terms as well as the cyclic variation of the baroclinicity between maximum and minimum are clear indications that the magnetic field is influencing the entropy (temperature) distribution. This is in close agreement with the physical mechanism of magnetic modulation of the thermal flux transport recently proposed by Cossette et al. (2017). We are currently investigating the evolutionary patterns of this thermal modulation.

The quadrupolar cell pattern in the stable layers is also visible in Fig. 7 and is associated with the large-scale magnetic torque. Its presence in Fig. 13 reflects the adjustments in the MC in response to GP.

6. Conclusions and final remarks

In this paper we have examined in detail the dynamical driving of meridional circulation in numerical simulations of solar convection, with and without a large-scale magnetic cycle. This work is motivated by the fact that the internal meridional flow in the Sun is a weak flow (a few m s^-1) compared to convection and differential rotation, and therefore quite difficult to measure helioseismically. Yet, this large-scale flow is believed to play a key role in the class of solar cycle models known as flux transport dynamos: it sets the cycle period in the advection-dominated regime, and its temporal variations on these decadal timescales are believed to greatly influence the amplitude of activity cycles. Moreover, the surface component of this flow contributes to the poleward transport of photospheric magnetic flux released by the decay of active regions, and thus has a direct – and observed – impact on the reversal and amplitude of the surface dipole moment.

In a thick, rotating, stratified turbulent fluid layer, the meridional flow dynamics are closely coupled to the balance of angular momentum. Consequently, we analyzed the evolution of the angular momentum as well as the meridional force balance in two analog simulations of global solar convection. The first is purely hydrodynamic (i.e., unmagnetized), while the second includes magnetic fields and self-consistently generates a large-scale magnetic cycle undergoing regular polarity reversals on a multi-decadal timescale. The comparison between the two simulations highlights the role of magnetism as a driver of meridional flow, and of its spatiotemporal variations.

The results obtained in the HD regime are in good agreement with previous analyses, and indicate that the convective angular momentum transport is responsible for the establishment of the large-scale meridional flow through the mechanism of gyroscopic pumping (GP). This mechanism reflects mainly the action of Reynolds stresses in areas where the rotation profile presents strong gradients, here primarily outside the TC.

The MHD simulation used in this study is the millennium simulation presented in PC14. It exhibits a number of solar-like features, including cyclic large-scale magnetic activity as well as cyclic variability patterns in the large-scale flows (torsional oscillations and MC variations) and convective heat flux. However, in this simulation, these characteristics appear in a range of latitudes that spans from 45° to 85°, i.e., inside the TC, away from the low latitude strong cylindrical differential rotation, while in the Sun they occur at lower latitudes (where cylindrical rotation is not observed). However, if we compare the active latitudinal range in the simulation to the Sun’s active latitudes, then the spatiotemporal patterns of rotational torsional oscillations become quite solar-like (see, e.g., Beaudoin et al. 2013), and the temporally averaged meridional flow pattern becomes remarkably similar to the helioseismic measurements of Zhao et al. (2013; see Passos et al. 2015).

Using the same methodology used to analyze the HD simulation we find, in the MHD case, that gyroscopic pumping is strongly influenced by the large-scale magnetic field. In Sect. 4 we showed that this influence materializes via a magnetic torque located at mid-latitudes at the bottom of the CZ. This magnetic torque accelerates and decelerates bands of rotation situated polewards and equatorwards of the toroidal field bands. Essentially, the magnetic field briefly changes the differential rotation, and the system re-establishes equilibrium by continuously altering the MC cell morphology in all of the CZ. The timescale associated with the recovery of this balance, i.e., the MC response, is on the order of one or two rotations (months). This is to be contrasted with the timescale associated with the variation of the large-scale magnetic field (several years). Since the system readjusts quickly, we can interpret the variation along the magnetic cycle as a quasi-static process where the system is always very close to equilibrium. We also note that the area where the magnetic torque is concentrated, and where the core of the GP mechanism takes place, is situated at the bottom of the CZ, away from the boundaries of the simulation domain, and therefore unlikely to be influenced by boundary conditions.

As the MC readjusts its structure to maintain angular momentum balance, it also has to satisfy the meridional force balance condition. In Sect. 5 we showed how this MFB is influenced not only by baroclinicity, but also by a magnetic contribution. Moreover, we find evidence that the baroclinic term itself is being modulated by magnetism as well. One of the ways this may happen is by the modulation of heat flux transport, as discussed in Cossette et al. (2017). The presence of magnetic field in certain areas alters the way heat is transported by convection and establishes thermal gradients. In the millennium simulation we observe a cyclic temperature variation pattern where the poles get cooler during cycle maxima and hotter than average at cycle minima. We are currently investigating whether the small CW rotation cell that appears near the poles at cycle maximum is a consequence of this MFB constraint in the adjustment of the MC.

One of the most prominent features resulting from these MC variations is a horizontal convergence of fluid into the equatorward edge of the toroidal magnetic field bands building up at the base of the CZ, and the associated mid-latitude upwelling it generates. This upwelling waxes and wanes in phase with the cyclic evolution of the magnetic field, becoming most prominent at cycle maximum. It may interfere with equatorward flux transport, promote flux emergence, and more generally affect the dynamical coupling between the convection zone and underlying radiative core over long timescales.

Another noteworthy specific spatiotemporal MC variation pattern (among the several) that can be extracted from this numerical experiment merits attention. In the top layers of our simulation domain, between 50° and the poles, we observe a poleward flow (see Fig. 10A). This surface flow exhibits a characteristic temporal modulation pattern, due to the appearance of an equatorward flow at high latitudes, peaking at cycle minima. This pattern is associated with a decrease in the magnetic torque at these latitudes near the surface and the appearance of a CW rotating MC cell (in the NH) near the poles (see Figs. 3B and C). A similar pattern is observed at the solar surface. The observational evidence (Haber et al. 2002; Ulrich & Boyden 2005; Ulrich 2010; Hathaway & Upton 2014; Bogart et al. 2015) indicates that this counter-cell tends to appear in the descending and minimum phase of the cycle. This cannot be explained by localized surface inflows because at those latitudes there are no active regions.

Generally speaking, the magnetic field alters the characteristics of convection and mean flows on both small and large-scales. What we observe in this simulation is a general magnetic modulation of convection and its associated dynamics during cycle rise and maximum, and a subsequent relaxation towards an HD-like profile when the cycle drops to a minimum. This raises concerns regarding the kinematic approach generally used in mean field and mean field-like axisymmetric dynamo models of the solar cycle. If the MHD effects that we see in this simulation scale up to solar conditions, then the kinematic approximation might be missing important physical effects.

Our analyses have led to a dynamically consistent scenario for the spatiotemporal evolution of the large-scale meridional flow in a simulated solar convection zone. However, this scenario was established on the basis of numerical simulation results carried out in a physical parameter regime far removed from solar internal conditions, so that one may legitimately question their relevance to the real Sun and stars. As a magnetized fluid system, our simulation does generate behaviors resembling solar observations, most notably decadal large-scale magnetic cycles, a reasonably solar-like internal differential rotation and pattern of torsional oscillations, and a high-latitude pattern of surface meridional flow variations. This suggests – certainly without proving – that the overall dynamical interactions between cyclic magnetism, angular momentum balance, and thermal wind balance taking place in the simulation do capture similar effects taking place in the solar interior. One potentially testable prediction emerging from our analysis is the buildup of a large-scale upwelling starting deep with the convection zone at active latitudes and peaking at cycle maximum. Because it is sustained over a time period commensurate with that of the magnetic cycle, such a spatiotemporally coherent upflow may actually be detectable helioseismically. We leave open the search for an associated helioseismic signature in extant data as an interesting observational challenge.

Acknowledgments

D. Passos is thankful to Sandra Braz for support, and acknowledges the financial support from the Fundação para a Ciência e Tecnologia grant SFRH/BPD/68409/2010 (POPH/FSE), CENTRA-IST, the GRPS-UdeM and the University of the Algarve for providing office space. P. Charbonneau is supported primarily by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada. All EULAG-MHD simulations reported upon in this paper were performed on the computing infrastructures of Calcul Québec, a member of the Compute Canada consortium. The National Center for Atmospheric Research is sponsored by the National Science Foundation.

References

Appendix A: Meridional force balance equation

The equation for the meridional force balance (13) is derived by computing the zonally averaged component of (12). In this appendix we detail the derivation of each individual term of (13) for our spherical latitudinal coordinate system.

i) Vorticity stretching

(A.1)

ii) Vorticity advection

(A.2)

iii) Vorticity "compressibility"

(A.3)

iv) Baroclinicity

(A.4)

v) Magnetic contribution 1

(A.5)

vi) Magnetic contribution 2

(A.6)

All Figures

Fig. 1 Meridional plane displays of time averaged ⟨ ur ⟩, ⟨ uθ ⟩ and stream function. The top row is for a HD simulation, while the bottom row shows the same quantities for a MHD counterpart. For the HD (MHD) case, the temporal average corresponds to 246 (406) yr. In panels A and D, red/yellow (blue/green) denotes flows towards the north (south). In panels B and E, red (blue) denotes rising (sinking) flows. In panels C and F, we present the stream lines of the MC. The red (blue) tones denotes circulation in the counterclockwise (clockwise) direction. The black dashed line indicates the depth of the base of the convection zone. Below it we have the stable layers. The vertical dotted line in C) indicates the position of the tangent cylinder for reference. At the top layer, thin ticks mark 15° intervals. Panel C can be directly compared with Fig. 1g of FM15.

In the text

	Fig. 2 Magnetic cycle proxy (for the northern hemisphere). The vertical gray areas represent the maxima and minima phases considered. For reference we labeled these phases for simulated cycle number 24.
In the text

	Fig. 3 Zonally averaged magnetic field components sampled at the minimum (top row) and maximum (bottom row) of cycle 25. Panels A and D show ⟨ Bφ ⟩ with the MC stream function overplotted, where the dashed (solid) contour lines represent cells circulating in the counterclockwise (clockwise) direction. Panels B and E show ⟨ Bθ ⟩, and panels C and F show ⟨ Br ⟩.
In the text

	Fig. 4 Panel A: ρ0∂ℒ /∂t, B) ⟨ ρ0um ⟩ ·∇ℒ and C) in kg m-1 s-2, averaged over 246 yr. Panels B and C have the same scale (105) but A) is two orders of magnitude smaller (103), indicating that a steady state has been achieved.
In the text

	Fig. 5 Panels A and B: same terms as in Fig. 4 for the MHD simulation also in kg m-1 s-2. Panel C: entire right-hand side of Eq. (7), which include the contributions from FRS, FMS, and FMT. All quantities are averaged over ten cycles (406 yr).
In the text

	Fig. 6 Components of the right-hand side of Eq. (10) and Fig. 5C, plotted individually: A) , B) , and C) .
In the text

	Fig. 7 Several terms of Eq. (7) averaged over ten cycle minima (A–D) and 10 cycle maxima (E–H). We focus only on the NH for simplification. Panels A and D: left-hand side of Eq. (7); B and F: ; C and G: ; and D and H: . All the contours have the same color scale as in Fig. 6. The vertical dotted lines indicate the TC for reference.
In the text

Fig. 8 Panels A and B: right-hand side of Eq. (7) in the NH averaged over cycle minima and cycle maxima. The vector field represents the direction of the meridional flow and the vertical dotted line represents the TC. The white contour lines in B) show the area where the toroidal field accumulates at cycle maximum. The bottom panels show the time evolution of the left-hand side of Eq. (7) and the individual terms on the right-hand side (as labeled), sampled at the numbered locations. The data in the lower three panels were smoothed with a one-year average filter.

In the text

Fig. 9 Panel A: individual terms of the left-hand side of Eq. (7) (in black and blue) and the total contribution of the right-hand side (orange) sampled at region1 ◯for a five-year interval, taken at cycle maximum. Panel B: individual terms of the left-hand side (black and purple) and right-hand side (red) of the meridional force balance (Eq. (13) in Sect. 5) for the same location and for the same time interval.

In the text

Fig. 10 Mean latitudinal flow ⟨ uθ ⟩ as a function of latitude and time at A) the top layers (r = 0.95 R⊙) and B) near the base of the CZ (r = 0.72 R⊙). For a better contrast, in panel B the displayed quantity is actually 2⟨ uθ ⟩. The black contour lines represents ⟨ Bφ ⟩ (solid and dashed for positive and negative polarities respectively). The horizontal solid line represents the latitude at which panel C is sampled. The time interval covers cycles 23 to 28. The color scale saturates at ± 3 m s-1 (red northward, blue southward). Panel C is a radius vs. latitude plot of ⟨ uθ ⟩ taken at 65° north. Both ⟨ uθ ⟩ and ⟨ Bφ ⟩ are smoothed over 6 months. The vertical dotted lines represent the maximum and minimum of cycle 25. The dashed line marks the tachocline depth, the dotted line below it is the depth where ⟨ Bφ ⟩ contours are computed and the solid line is the depth where panel B is sampled.

In the text

Fig. 11 Panels A and D: left- and right-hand sides of the MFB equation for the HD simulation. The terms were computed for the NH only, and with color saturation at ± 3 × 107 s-1. The vertical dotted line represents the TC. In the other four panels we present the individual contributions on the right-hand side: B) stretching, C) advection, E) compressibility and F) baroclinicity (with color saturation at ± 2 × 10-7 s-1).

In the text

	Fig. 12 Panels A and D: left- and right-hand sides of the MFB equation in the NH for the MHD case. The individual terms of the right-hand side are shown in panels B, C, E, and F. The magnetic contribution is the sum of the two individual magnetic terms. The vertical dotted line represents the TC. The color scale saturates at ± 1.5 × 10-7 s-1.
In the text

	Fig. 13 Panels A to C: left-hand side, baroclinic, and magnetic contributions of the MFB averaged over 10 cycle minima. Panels D to F: same quantities as the top row averaged over ten cycle maxima. The color scale saturates at ± 1.5 × 10-7 s-1.
In the text