© ESO, 2016

1. Introduction

The Sun has an activity cycle of about 11 yr, with an underlying magnetic field that oscillates with a period of around 22 yr. A dynamo operating in the convection zone below the solar surface is responsible for generating cyclic magnetic fields (see e.g., Brandenburg & Subramanian 2005; Charbonneau 2005, and reference therein). The occurrence of sunspots, which is the main surface manifestation of the solar cycle, varies regularly over the cycle. At the beginning of each cycle sunspots tend to appear at mid-latitudes, while toward the end of each cycle they tend to appear at low latitudes. It is believed that these sunspots and their occurrence are connected to an underlying toroidal magnetic field, which is migrating equatorward during the cycle. Theoretical models of solar magnetic field evolution have been studied for several decades. Mean-field models, where turbulence effects are parameterized through transport coefficients (see e.g., Krause & Rädler 1980), have been successful in producing some observed magnetic field properties (e.g., Käpylä et al. 2006; Kitchatinov & Olemskoy 2012). Another class of dynamo models relies on the Babcock-Leighton effect (Babcock 1961; Leighton 1964) and flux transport by meridional circulation (e.g., Choudhuri et al. 1995; Dikpati & Charbonneau 1999). The low Reynolds numbers compared with the Sun limit the usefulness of global simulations of self-consistent convection, where the equations of magnetohydrodynamics are solved directly. However, the increasing computing power has led to the successful reproduction of some observed features of the solar magnetic field by such models.

For a long time global simulations were only able to generate either poleward migrating fields (Gilman 1983; Brun et al. 2004; Käpylä et al. 2010; Nelson et al. 2013) or oscillatory ones with no clear migration pattern (Ghizaru et al. 2010; Racine et al. 2011). For the first time, Schrinner et al. (2011) and Käpylä et al. (2012) could produce clear equatorward migration of the toroidal magnetic field. Stratification and rotation rate had to be high enough for this to work (Käpylä et al. 2013). Recently, several groups have been able to produce grand minima-type events (Passos & Charbonneau 2014; Augustson et al. 2015; Käpylä et al. 2016a). In Käpylä et al. (2016a), a secondary dynamo mode disturbed the surface field of the primary mode. Warnecke et al. (2014) could explain the equatorward migration seen in simulations of Käpylä et al. (2012, 2013) and those of Augustson et al. (2015) as a propagating αΩ dynamo wave following the Parker-Yoshimura rule (Parker 1955; Yoshimura 1975). Here, α was estimated via the kinetic helicity and Ω is the local solar rotation rate. An equatorward migrating dynamo wave is possible if α is positive (negative) in the northern (southern) hemisphere (Steenbeck et al. 1966) and the radial gradient of Ω is negative. These interpretations have been confirmed independently by Duarte et al. (2016) and by the computation of turbulent transport coefficients obtained through the test-field method Warnecke et al. (2016b).

Given that this simple relation can describe the behavior of dynamos driven by self-consistent turbulent convection in simulations, the Parker-Yoshimura rule can also be a possibility to explaining the equatorward migration of the magnetic field of the Sun. Indeed, the differential rotation of the Sun (Schou et al. 1998) has a negative radial gradient in the near-surface shear layer (Thompson et al. 1996; Barekat et al. 2014, 2016). This makes it a possible location of the solar dynamo (Brandenburg 2005).

The generation of the solar differential rotation can well be described by mean-field models, where the productive parts of the off-diagonal Reynolds stress are parameterized by the Λ effect (Rüdiger 1989) and the turbulent heat transport in terms of anisotropic turbulent heat conductivity (see also Brandenburg et al. 1992). These models can reproduce the surface differential rotation, the spoke-like rotation profile and the near-surface shear layer (e.g., Kitchatinov & Rüdiger 1995, 2005; Rüdiger et al. 2013). With global models of turbulent convection, it is challenging to generate such rotation profiles (e.g., De Rosa et al. 2002; Brandenburg 2007). Miesch et al. (2006) were able to produce a spoke-like rotation profile by imposing a latitudinal entropy gradient at the bottom boundary, while Brun et al. (2011) used a stably stratified layer below the convection zone. Guerrero et al. (2013) could produce a near-surface shear layer at lower latitudes. Recently, Hotta et al. (2015) used a reduced sound speed technique (Hotta et al. 2012) achieving high stratification to produce a near-surface shear layer, the generation of which was suggested to be related with the meridional Reynolds stress component Q_rθ. However, they still had cylindrical rotation contours within the convection zone.

Warnecke et al. (2013a) applied a different approach and used an outer coronal envelope above the dynamo domain (Warnecke & Brandenburg 2010) to reproduce spoke-like differential rotation at low latitudes with a weak near-surface shear layer. This two-layer approach has also been used to successfully simulate coronal ejections driven by dynamos arising from forced turbulence (Warnecke et al. 2011, 2012a) as well as by convective dynamos (Warnecke et al. 2012b). Furthermore, the outcome of dynamo simulations suggests that the presence of a coronal envelope supports the dynamo and leads to a higher field strength (Warnecke & Brandenburg 2014).

To investigate and follow up on the findings of Warnecke et al. (2013a), we perform a detailed study of similar simulations with and without a coronal envelope to investigate the effect of a coronal envelope as a free boundary on a convectively driven dynamo. We vary the size of the envelope, as well as the cooling profile, the magnetic boundary condition, and the rotation rate. We analyze the effect on the flows, differential rotation, and the magnetic field evolution. Even though the solar corona most likely has limited influence on the dynamics of subsurface flows and the evolution of the magnetic fields in the Sun, these studies are important for investigating different influences and effects on convective dynamo simulations. Every simulation, in which we better understand the mechanism causing flow and magnetic field evolution, will bring us a step closer toward understanding the dynamics of the interior of the Sun and other stars.

2. Model and setup

Our setup is similar to the one-layer model of Käpylä et al. (2012, 2013) and the two-layer model of Warnecke et al. (2013a), both of which have recently also been used in Warnecke et al. (2014). We use a wedge in spherical polar coordinates (r,θ,φ), in which the layer below the surface (r₀ ≤ r ≤ R) represents the convection zone. Here, R is the solar radius and r₀ corresponds to the bottom of the convection zone at r = 0.7 R. The layer above the surface (R ≤ r ≤ R_C) represents a simplified coronal envelope, which extends to outer radius R_C. The domain spans 15° ≤ θ ≤ 165° in colatitude and 0° ≤ φ ≤ 90° in longitude. We solve the equations of compressible magnetohydrodynamics, (1)(2)(3)(4)where the magnetic field is defined via the vector potential B = ∇ × A, ensuring the solenoidality of B at all times, J = ∇ × B/μ₀ is the current density with μ₀ being the vacuum permeability, η is the magnetic diffusivity, ν is the kinematic viscosity, u is the plasma velocity, ρ is the mass density, s is the specific entropy, and D/Dt = ∂/∂t + u·∇ is the Lagrangian derivative. The traceless rate-of-strain tensor is given by (5)where semicolons denote covariant differentiation; see Mitra et al. (2009) for details. The gravitational acceleration is given by (6)where G is Newton’s gravitational constant and M is the mass of the Sun. In addition Ω₀ = Ω₀(cosθ,−sinθ,0) is the rotation vector, where Ω₀ is the rotation rate of the comoving frame. Using the ideal gas law, the pressure is given by p = (γ−1)ρe, where γ = c_P/c_V = 5/3 is the ratio of specific heats at constant pressure and constant volume, respectively, and e = c_VT is the internal energy density, which is related to the temperature T. The two diffusive heat fluxes are defined as (7)where is the radiative heat flux with the radiative heat conductivity K and is the subgrid scale (SGS) heat flux that carries the unresolved turbulent heat flux of convection with the SGS heat diffusivity χ_SGS; see Käpylä et al. (2013) and Warnecke et al. (2013a) for details. Finally, the function Γ_cool relaxes the temperature toward a predefined profile T_ref(r)(8)where Γ₀ is a cooling luminosity. f(r) is a profile function tending to unity in r>R and going smoothly to zero in r ≤ R, see Warnecke et al. (2013a) for details. Figure 1 shows the corresponding temperature and density stratifications for the runs in Set A, see Table 1.

Fig. 1 Radial profiles of azimuthally and latitudinally averaged temperature ⟨ T ⟩ θφa) and density ⟨ ρ ⟩ θφb) normalized by their values at the bottom of the domain, T0 or ρ0, respectively for Set A. The inlays show the radial profiles near the surface.

We use isentropic, hydrostatic initial conditions, as in previous models (Käpylä et al. 2013; Warnecke et al. 2013a). This initial setup is not in thermal equilibrium, but the flux at the lower boundary exceeds the flux leaving at the outer boundary, resulting the onset of convective instability. Furthermore, we initialize the magnetic field as a white noise seed field in the convection zone. We apply periodic boundary conditions in the azimuthal (φ) direction. For the velocity field we apply stress-free boundary conditions at the radial and latitudinal boundaries. The magnetic field follows a perfect conductor condition at the lower radial and at the latitudinal boundaries. We force the field to be radial at the top boundary. Furthermore, the temperature gradient at the bottom boundary is fixed to have a constant heat flux into the domain, and the latitudinal boundaries are impermeable for heat fluxes. On the upper radial boundary we either apply a black body condition, (9)or a constant temperature (10)In the former case the heat is transported out of the domain via an enhanced SGS flux; see Käpylä et al. (2013) and in the latter case via a cooling flux; see Warnecke et al. (2013a).

We apply roughly 10–50 times larger values of viscosity and magnetic diffusivity in the coronal envelope compared to the convection zone to avoid high velocity amplitudes and strong shear flows resulting in small-scale magnetic field enhancements. For the transition region we use a hyperbolic tangent function with a radial width of 0.01 R at r = 1.06 R, which is the largest value feasible. In the case of Run A5, the viscosity in the coronal envelope is just three times higher than in the convection zone, but we apply an additional shock viscosity and shock diffusivity; see Haugen et al. (2004) and Gent et al. (2013) for details regarding their implementation. The heat conductivity K in the coronal envelope is chosen such that the heat diffusivity χ = K/ρc_P is constant.

We characterize the runs by the values of the input parameters: Prandtl number Pr = ν/χ^m, sub-grid scale Prandtl number , magnetic Prandtl number Pr_M = ν/η, Taylor number Ta = (2Ω₀(R−r₀)²/ν)², and Rayleigh number, (11)which is obtained from a hydrostatic one-dimensional model for the same initial setup with χ^m = χ(r = 0.85 R) and . Furthermore we define the fluid and magnetic Reynolds numbers, Re = u_rms/νk_f and Re_M = u_rms/ηk_f, respectively, the Coriolis number Co = 2Ω₀/u_rmsk_f and the Péclet number , where k_f = 2π/ (R−r₀) ≈ 21 /R is used as a reference wavenumber and where u_rms is the typical turbulent velocity in the convection zone defined as (12)which corrects for the removal of the differential rotation-dominated u_φ. The slow mean meridional flows, however, are not removed. Azimuthal averages combined with time averages in the saturated stage are referred to as mean and indicated with an overbar, e.g., , while other averages are indicated as ⟨ . ⟩ with the spatial directions as indices. The index 0 refers to the value at the bottom of the domain, that is ρ₀ and T₀. We also use the meridional distribution of turbulent velocities (13)where the fluctuating velocity is defined via . Thus, meridional mean flows are here removed. The mesh is equidistant in all directions. The grid resolutions are given in Table 1.

We express our results in physical units following Käpylä et al. (2013, 2014) and Warnecke et al. (2014) by choosing a normalized rotation rate , where Ω_⊙ = 2.7 × 10^-6s^-1 is the solar rotation rate. The simulations were performed with the Pencil Code¹, which uses a high-order finite difference method for solving the compressible equations of magnetohydrodynamics.

Fig. 2 Radial profiles of azimuthally and latitudinally averaged turbulent rms velocity a) and ⟨ s ⟩ θφb) in m/s or normalized by cP respectively for Set A. Mean radial mass flux c) through the surface (r = R) for Set A. The horizontal line indicates the zero value and the three vertical thin lines indicate the equator (θ = 90°) and the intersection with the inner tangent cylinder (θ−90° ≈ ± 45° latitude).

3. Results

In this work we compare and analyze 13 runs divided into two sets based on the rotation rate. In Set A the normalized rotation rate is whereas the runs in Set B have a rotation rate of . For both sets we investigate the effects of a cooling layer and the blackbody boundary condition as well as the size of the coronal envelope. A summary of the runs can be found in Table 1 and the stratifications of temperature, density and entropy of Runs A are shown in Figs. 1 and 2. Run A1 is nearly the same as Run B4m of Käpylä et al. (2012), Run C1 of Käpylä et al. (2013), Run I of Warnecke et al. (2014), and the run of Warnecke et al. (2016b), where a blackbody boundary condition is used. However, we choose a slightly higher stratification and a slightly lower value of Pr_SGS, namely 2 instead of 2.5; therefore Run A1 is the same as Run D3 of Käpylä et al. (2016c). In Runs A1c and A1c2, the blackbody boundary condition is replaced by a shallow (R ≤ r ≤ R_C = 1.01 R) cooling layer, where in the former case the temperature minimum is below the surface (r = R) and in the latter above the surface. These two runs have been recently used in Warnecke et al. (2014) as Runs III and IV, respectively. The only difference between Runs A1pc and A1 is the use of a perfect conductor condition instead of a radial field condition for the magnetic field at the top boundary.

The other runs of Set A have a coronal envelope with different outer radii R_C. Runs A2, A3, and A4 have the same cooling function as Run A1c, where the temperature reaches a minimum below the surface. The temperature increases to a constant coronal value, which is more than twice the value at the bottom of the convection zone; see Fig. 1a. This results in a positive entropy gradient above r = 0.97 R, where the convection ceases and drops by a factor of two; see Figs. 2a, b. In Run A3t, the same cooling function is applied as in Run A1c2, leading to a temperature minimum above the surface at r ≈ 1.01 R. This causes the entropy gradient to become positive at the surface (r = R) and an increase of all the way to the surface; see Figs. 2a, b. Already here, we can state that the use of cooling profiles in Runs A1c and A3t reproduce most properties of the density, temperature and entropy stratification as the blackbody boundary conditions in Run A1. In Run A5, we lower Pr_SGS and Pr_M to 0.5 and therefore increase the fluid Reynolds number. As the heat flux at the boundary is the same as in the other runs, increases only slightly, see Fig. 2a, leading to a reduced Coriolis number; otherwise Run A5 is similar to Run A3t. However, the use of a lower value of viscosity in combination with a shock viscosity allows higher velocities in the coronal envelope as shown in Fig. 2a. The runs of Set B are essentially the same as the corresponding runs in Set A with a lower rotation rate. The radial temperature, density, entropy and velocity profiles behave similarly as in Set A and are therefore not shown here. In the following we investigate the influence of the coronal layer on mass flux and temperature distribution (Sect. 3.1) as well as on differential rotation and meridional circulation (Sect. 3.2). Furthermore we discuss Reynolds stresses and the Λ effect (Sect. 3.3) and their contribution to the thermal wind balance (Sect. 3.4). Then we investigate the influence of the magnetic top boundary on the field structure near the surface (Sect. 3.5) and the magnetic field evolution (Sect. 3.6).

3.1. Mass flux and temperature distribution at the boundary

We begin by inspecting the influence of the top boundary on the mass flux and the temperature distribution; see Fig. 2c. Except for Run A1, the mass flux at r = R is nonvanishing, showing a strong latitudinal dependency. Near the equator it is positive (outflow) and at latitudes around ±30° negative (inflow) suggesting a circulation in the coronal envelope. At mid-latitudes the mass flux becomes positive again, but is fluctuating around zero toward higher latitudes. The flow structure is strongly influenced by the rotation, as seen from the alignment with the inner tangent cylinder; see Fig. 2c. Although the mass flux in the runs with an extension above the surface has non-zero values, they are small in comparison to at the surface. Furthermore, there is no qualitative difference between runs with a cooling layer R_C = 1.01 R (Runs A1c and A1c2) and a coronal envelope R_C ≥ 1.2 R (Runs A2, A3, A3t, and A4). In fact, Run A1c2 has the highest mass flux through the surface; see Fig. 2c. Furthermore, Run A5 is similar to the other runs showing more fluctuations as a function of latitude but with a similar magnitude of variations as in the other runs. This shows that the influence of a coronal envelope via a radial mass flux is small and can be neglected. Furthermore, we find no indication that the viscosity profile has a major influence on the mass flux.

As a second step we investigate the latitudinal temperature variation, , at two radii; see Fig. 3. In general the surface perturbations are strongest near the poles, decrease toward a minimum at mid-latitudes, and increase again below ± 20° latitude. However, in the middle of the convection zone the temperature minimum is at the equator. This is a clear indication of a strong rotational influence on the temperature distribution. Run A1 shows the largest relative temperature perturbations; over 0.3 near the poles and up to −0.1 at ± 20° latitude. Here, the radiative boundary condition lets the temperature at the surface evolve more freely, which leads to this strong variation.

Fig. 3 Latitudinal dependence of at radius r = 0.98 Ra) and at radius r = 0.85 Rb) for Set A. Linestyles as in Fig. 1. The thin black lines indicate the zero value, the equator (θ = π/ 2) and the location of the inner tangent cylinder.

Runs A1c and A1c2 have a similar distribution, but with around three times smaller values. There the cooling layers cool the temperature to a certain latitudinally independent value with a relaxation time equal to the turnover time. This leads to a reduction of the temperature perturbations near the surface. The temperature difference ΔT is significantly reduced for all runs with a coronal envelope. However, the temperature is still higher near the poles and lower near the equator than the latitudinal average. In these runs the coronal envelope with its mass and heat capacity serves as a buffer in smoothing the temperature at the surface. The influence of the cooling layer and the coronal envelope seems to penetrate also deeper down and influences the temperature variations in the middle of convection zone. The difference in the temperature profiles caused by the coronal envelope can also influence differential rotation; see Sect. 3.2.

In the higher Re and thus more turbulent Run A5, the temperature variation is similar to that of the other runs with a coronal envelope. The runs with slower rotation (Set B, not shown here) show a similar behavior, but the latitudinal temperature perturbations are weaker than in the more rapidly rotating runs, due to the reduction of the rotational influence.

On the solar surface, systematic latitudinal temperature variations have not been observed. However, variations in the range of a few kelvin, so less than 0.1% of the surface temperature, are below the measurable range. This value is exceeded in all of our simulations. The temperature difference between the poles and equator can influence the differential rotation in the Sun (Rüdiger 1989) and in simulations (e.g., Miesch et al. 2006; Warnecke et al. 2013a).

3.2. Differential rotation and meridional circulation

Fig. 4 Angular velocity Ω(r,θ)/Ω0 for Runs A1 and A3 (top row), Runs A3t and A5 (middle row) and Runs B1 and B3t (bottom row). Here, is the local rotation rate. The dashed lines indicate the surface (r = R).

Fig. 5 Differential rotation Ω(r,θ)/Ω0 as a function of radius for all runs of Set A at θ = 75° (15° latitude) a), and at the equator θ = π/ 2b), and Set B at mid-latitudes θ = 75°c), and at the equator θ = π/ 2d). The inlay in b) shows the angular velocity near the surface for Runs A2, A3, A3t, A4, and A5.

As shown in Warnecke et al. (2013a), simulations with a coronal envelope are more capable of reproducing a spoke-like differential rotation profile than runs without a coronal envelope using similar parameters (Käpylä et al. 2013). In this work, the parameters of the runs are nearly identical, so we can isolate the influence of the coronal envelope. In Fig. 4 we plot the angular velocity, , in the meridional plane and in Fig. 5 as a function of radius at two latitudes. In all of the runs the equator rotates faster than the poles, similarly to the other simulations with similar Coriolis numbers (Brown et al. 2010; Käpylä et al. 2013; Warnecke et al. 2013a; Augustson et al. 2015). Furthermore, all runs possess a local minimum in the angular velocity Ω at mid-latitudes, which has been shown to facilitate equatorward migration (Warnecke et al. 2014; Käpylä et al. 2016c). In runs without coronal envelopes the rotation profile has a similar structure; the contours of constant angular velocity show a strong alignment with the rotation axis, following the Taylor-Proudman theorem. Run A1 possesses the strongest differential rotation. In particular at the equator the surface rotates faster than in the other runs. Runs A1c and A1c2 show a similar radial dependency, but differential rotation is weaker than in Run A1; see Figs. 5a and b.

In the runs with a coronal envelope the Taylor-Proudman balance is broken and the rotation profile has a more spoke-like shape; see Fig. 4. The differential rotation and the local minimum at mid-latitudes in these runs is much weaker. The minimum also occurs at a greater depth than in the runs without a corona. The dependence on the size of the coronal extent on differential rotation is weak – similarly to what is seen for the latitudinal temperature distribution; see Fig. 3. The cooling profile has a minor influence on the rotation profiles such that in Run A3t the contours of constant rotation are slightly more radial than in Run A3 In Run A5, the differential rotation is even weaker with more radial contours of rotation.

The overall rotation profiles of Set B are similar to those of Set A and also show a local minimum of Ω at mid-latitudes although it is weaker and at a greater depth; see Fig. 4. Also here, the presence of a coronal envelope leads to more radial contours of rotation. However, the amount of differential rotation is reduced only for Run B3t, but not for Run B3.

Fig. 6 Meridional circulation in terms of the mass flux for Runs A1 and A3t. The dashed lines indicate the surface (r = R) and the red solid line the inner tangent cylinder.

If we compare the differential rotation profiles with those of Warnecke et al. (2013a), we see that the contours of constant rotation are more cylindrical here. In particular, moderate rotation runs had more spoke-like differential rotation than rapidly rotating ones which is not seen in the present work. Therefore, we relate the more spoke-like rotation profiles in Warnecke et al. (2013a) to a lower density stratification as the other parameters are similar.

For most of the runs (A2, A3, A3t, A4, B1, B3, and B3t) the maximum of rotation at the equator is actually below the surface, indicating a near-surface shear layer with negative radial shear. The logarithmic gradient of rotation, dlnΩ/dlnr, is around −0.2 near the surface for Run A4, and −0.15 for Runs A2, A3, and A3t. This is much weaker than the value for the Sun, which is dlnΩ/dlnr ≈ −1 for all latitudes (Barekat et al. 2014, 2016). In all runs of Set B, except for B1c, the near-surface shear region is more extended than the ones in Set A, and the gradient is stronger; dlnΩ/dlnr reaches values of −0.8 for Run B3 and −0.5 for Run B3t. It is expected that for runs with lower rotation rate, a near-surface shear layer is stronger due to the weaker influence of the Coriolis force near the surface; see Sects. 3.3 and 3.4. In agreement with Λ effect theory (Rüdiger 1980, 1989), the double-logarithmic gradient should only be close to −1 very near the surface where the local Coriolis number is small (Kitchatinov & Rüdiger 2005; Kitchatinov 2013; Rüdiger et al. 2014; Kitchatinov 2016). While this is true for the Sun, it is not the case in our simulations owing to limited stratification.

To investigate the influence of the cooling profile on the temperature variation and differential rotation, we have performed two additional runs, in which we either increased or decreased the cooling luminosity compared to Run A1c2 (not shown). A decrease of the cooling luminosity by a factor of two leads to a shift of the temperature minimum at higher radii. The mean temperature increases slightly, resulting in a higher density at the surface and a decrease of the density stratification in the convection zone. An increase of the cooling luminosity has the opposite effect, leading to a temperature minimum at a greater depth, and a lower temperature and density in the convection zone. Weaker cooling causes stronger differential rotation, especially at higher latitudes, while stronger cooling does not show a significant effect. The gradient dlnΩ/dlnr at the equator becomes more negative with a weaker cooling.

Differential rotation is also generated in the coronal envelopes. Below r = 1.01 R, the differential rotation follows the rotational behavior of the convection zone. Above r = 1.01 R the plasma rotates nearly uniformly near the equator with a rotational speed close the Ω₀. The mid-latitudes rotate faster than the equator and at high latitudes the coronal envelopes decrease to slower rotation. All runs exhibit strong cylindrical and radial shear. This is consistent with Warnecke et al. (2013a), where runs with lower stratification showed a similar behavior. The change of behavior from below to above r = 1.01 R is caused by the change in temperature and density stratifications and not by the viscosity profiles. This can be seen by comparing Runs A3t and A5, where the same temperature and density profiles are used, but a much lower viscosity is applied in Run A5 than in Run A3t.

In Fig. 6 we plot the meridional circulation in terms of the mass flux in the meridional plane, where is the meridional flow. The meridional circulation has a multi-cellular structure in all runs. Near the equator at the surface the flow is poleward, but it can become equatorward at high latitudes; see Fig. 6. The strongest contribution to the mass flux carried by the meridional circulation occurs within the bulk of the convection zone. There the flow is aligned with the rotation axis and streaming toward the equator along the inner tangent cylinder and toward higher latitudes further away from the rotation axis. These mass flows seem to stream toward the local minima of Ω at mid-latitudes. From there, most of the runs develop a flow toward the equator following the θ direction. The stronger meridional flows in Run A1 are due to the higher density, see Fig. 1b, while the actual flow is quite similar in all runs of Set A, see Fig. 8h. The runs of Set B generate stronger meridional circulation, similarly to what was found in Warnecke et al. (2013a). At these rotation rates, slower rotation leads to an increase of meridional circulation as found in mean-field models (Köhler 1970; Rüdiger 1989) and numerical simulations (e.g., Brown et al. 2008; Augustson et al. 2012). In general the meridional flow pattern does not change due to the influence of the coronal envelope.

3.3. Reynolds stresses and Λ effect

Fig. 7 From left to right: off-diagonal components of the Reynolds stress Qrθ, Qrφ and Qθφ normalized by νtΩ0 (first three columns), the turbulent viscosity in terms of molecular viscosity νt/ν (fourth column) and the three components of the Λ effect ΛM, ΛV, and ΛH (fifth to seventh column) normalized by νt for Runs A1 (top row), A3t (middle) and A5 (bottom).

Fig. 8 Off-diagonal components of the Reynolds stress Qrθa), Qrφb) and Qθφc) normalized by νtΩ0, the turbulent viscosity in terms of molecular viscosity νt/νd) and the three components of the Λ effect ΛMe), ΛVf) and ΛHg) normalized by νt as well as the meridional flow h) for Set A in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (r = R).

Differential rotation and meridional circulation in the Sun and other stars is generated by the interaction of turbulent convection and rotation (Rüdiger 1989). Reynolds stresses become anisotropic due to an angle between the gravity and the axis of rotation. The non-diffusive contribution of the Reynolds stress tensor Q_ij can be expressed via the Λ effect which produces equatorial acceleration if the angular momentum transport is directed equatorward (see e.g., Rüdiger 1989; Kitchatinov & Rüdiger 1995, 2005). It has been recently shown that there is strong evidence for the Λ effect operating in the Sun, causing the observed rapidly rotating equator (Rüdiger et al. 2014) and the latitude-independent surface shear (Kitchatinov 2016).

We calculate the three off-diagonal components of the Reynolds stress tensor , , and . Both Q_rφ and Q_θφ contribute to the angular momentum balance and their non-diffusive parts are associated with the horizontal and vertical Λ effects, respectively. Even though Q_rθ does not directly contribute to angular momentum transport, it has been argued to be important in generating a near-surface shear layer in global convection simulations (Hotta et al. 2015). The Reynolds stresses can be written as (e.g., Rüdiger 1989; Rüdiger et al. 2013) (14)(15)where Λ_V and Λ_H are the vertical and horizontal components of the Λ effect and ν_t is the turbulent viscosity (assumed isotropic). Following Käpylä et al. (2014) we approximate ν_t as (16)where α_MLT = 5/3 has been assumed for the mixing length parameter, and H_p = − [∂lnp(r,θ) /∂r] ^-1 is the pressure scale height.

A meridional Λ effect also exists (Pulkkinen et al. 1993; Rieutord et al. 1994; Käpylä et al. 2004; Käpylä & Brandenburg 2008) which is related to Q_rθ via (17)In the earlier definitions, no Ω factor was included in the first term because, on theoretical grounds, one expects the rotational effects on the meridional part of the Reynolds stress Q_rθ to be proportional to sinθcosθ (Rüdiger 1989; Rieutord et al. 1994). However, to obtain the same units for the coefficient Λ_M as for the other components of Λ, we include here an Ω factor. In Appendix A, we give a simplistic derivation for these coefficients, which confirms that Q_rθ is not only proportional to a higher power of Ω than the other two components of Λ, but that it also picks up a contribution proportional to . However, under nearly isotropic conditions, the contribution of Λ_M should have the same sign as ; see Appendix A. In our simulations this is indeed the case; see Fig. 7.

We can now solve Eqs. (15)–(17) for Λ_M, Λ_V and Λ_HWe plot the Reynolds stresses (Q_rθ, Q_rφ, and Q_θφ), the turbulent viscosity ν_t, and the three components of the Λ effect (Λ_M, Λ_V and Λ_H) in the meridional plane for Runs A1, A3t, and A5 in Fig. 7 and as a latitudinal cut (15° latitude) for all runs of Set A in Fig. 8. For Run A1, the Reynolds stresses show the usual behavior of rotating convection (Käpylä et al. 2011, 2016c). Q_rθ and Q_θφ are antisymmetric over the equator. In the northern hemisphere, Q_rθ is negative at the surface and positive deeper down, but only outside the inner tangent cylinder. The latitudinal variation of Q_rθ agrees with that found both by Pulkkinen et al. (1993) and Rieutord et al. (1994). Q_θφ is positive in the northern hemisphere near the surface and weakly negative deeper down at low latitudes. As expected, Q_rφ is symmetric over the equator (Rüdiger 1980). This stress component is mostly positive (negative) at low (high) latitudes with a small negative region at the equator. The meridional structure of all Reynolds stress components agrees qualitatively with Run A6 of Käpylä et al. (2011) and Hotta et al. (2015) with significantly lower and higher density stratifications, respectively. However, the peak values are half of those in Run A6 of Käpylä et al. (2011), which is likely due to the slower rotation in their study, as the stresses are known to be quenched for faster rotation (e.g., Rüdiger 1989; Rüdiger et al. 2013).

The turbulent viscosity ν_t has a maximum at the equator and towards the bottom of the convection zone. The value of ν_t/ν is close to Reynolds number Re = 35. Λ_M is non-zero only at low latitudes, being positive near the surface in the topmost 5–10 Mm, negative a bit deeper and mostly positive with a lower amplitude even deeper down. Λ_V shows strong alignment with the rotation axis, being positive (negative) outside (inside) a cylindrical radius of 0.85 R, with a negative region near the surface at the equator. Λ_H has a similar structure, but shows a concentration near the equator above 0.85 R, where Λ_H is three times stronger than Λ_V. A similar structure, but with three times lower values, has been found by Käpylä et al. (2014) and Karak et al. (2015) using the same technique, but for runs that rotate slower.

The presence of a coronal envelope significantly alters the Reynolds stresses and therefore ν_t and the Λ effect; see Figs. 7 and 8. The Reynolds stresses lose their alignment with the rotation axis. Q_rθ changes sign and is positive (negative) at lower latitudes in the northern (southern) hemisphere, which is similar to Hotta et al. (2015). This behavior has been found for all runs with an extended coronal envelope (A2, A3, A3t, A4, and A5), whereas in Runs A1, A1c, and A1c2, Q_rθ is negative close to the surface in the northern hemisphere; see Fig. 8a. This seems to confirm the presence of a correlation between a positive (negative) value of Q_rθ in the northern (southern) hemisphere and the generation of near-surface negative shear. The inclusion of a coronal layer changes Q_rφ such that the minimum at the equator disappears and the overall magnitude of the stress is reduced by roughly a factor of two; see Fig. 7. At 15° latitude, all runs show a similar behavior with the main variation being the location of the maximum near the surface depending on the depth of the cooling layer; see Fig. 8b. Q_θφ changes similarly as Q_rφ. The maxima near the surface are again shifted deeper in the runs with a deeper acting cooling layer; see Fig. 8c. The profiles of Q_rφ and Q_θφ of Run A3t are similar to Hotta et al. (2015).

The turbulent viscosity ν_t varies less as a function latitude in runs with a coronal envelope; see Fig. 7. Figure 8d shows that ν_t is reduced at the bottom of the convection zone and near the surface in the runs with a coronal envelope which is consistent with the lower turbulent velocities in those cases; see Fig. 2a. In Run A5, the profile is similar to the other runs with coronal envelope, but ν_t/ν is higher due to a lower value of ν.

All of the Λ-coefficients are reduced in runs with a coronal envelope; see Figs. 7 and 8e–g. This is consistent with a weaker differential rotation in these runs; see Fig. 5. Similarly as Q_rθ, also Λ_M changes sign in the equatorial regions in the runs with a coronal envelope. This is consistent with the mean meridional flows changing sign in these runs as can be seen in Fig. 8h. Λ_V does not change significantly apart from the reduced amplitude. The vertical Λ effect is thought to be responsible for generating radial shear. This is consistent with a stronger Λ_V producing a stronger radial differential rotation in Run A1 and with a correspondingly weaker Λ_V and weaker radial differential rotation in Runs A2, A3, A3t, A4, and A5; see Fig. 5. The horizontal Λ effect shows a more broad maximum as a function of latitude in the near surface layers in Run A3t in comparison to Run A1, whereas the amplitude is reduced by a factor of four. Figure 8 reveals that the values of Λ_H near the surface are lower in the runs with a coronal envelope or a cooling function penetrating the surface, consistent with weaker latitudinal shear. The changes due to the coronal envelope are even more pronounced in the anisotropy parameters A_M, A_V, A_H, which are related to the Λ effect; see Appendix B and Fig. B.1. A_M and A_V change from being mostly positive to mostly negative (also within the coronal layer), whereas the changes in A_H are not as strong. The results for the more turbulent Run A5 are very similar with those of Run A3t; see Fig. 7.

In Set B, where rotation is slower, the rotational influence on convection is reduced by roughly a factor of two. However, the only major difference to the runs in Set A occurs for Q_rθ, where the sign at low latitudes near the surface is positive already without a coronal envelope, although the difference in the meridional Λ effect is rather minor. The main effect of the coronal envelope on the Reynolds stresses as well as the components of the Λ effect is their reduced amplitude similarly as in Set A.

According to theory (e.g., Kitchatinov & Rüdiger 1995, 2005; Rüdiger et al. 2014; Kitchatinov 2016), the near-surface shear layer in the Sun is caused by the vanishing horizontal Λ effect and the sole contribution of a negative vertical Λ effect. In our simulations, we do indeed find a weaker Λ_H in some of the runs, where dlnΩ/dlnr is negative, see Figs. 7 and 8g, but for most of the runs the relation is inconclusive. In Sect. 3.4, we will investigate this in more detail.

In this section we have not discussed the influence of the magnetic field on differential rotation generators. This has been discussed in detail in Karak et al. (2015) for similar simulations, but without coronal envelope. They found that the large-scale and turbulent Maxwell stresses are around two orders of magnitude smaller than the Reynolds stresses and therefore do not influence the angular momentum transport significantly. However, the recent study by Käpylä et al. (2016c) has revealed that the turbulent Maxwell stresses in similar setups without coronal envelopes are of the same order of magnitude as the Reynolds stresses at high magnetic Reynolds numbers. The current runs are in the intermediate range in Re_M where the total stress is already affected but the qualitative character of the hydrodynamic results remains unchanged.

3.4. Thermal wind balance

The Taylor-Proudman balance can be broken by a non-zero baroclinic term in the mean azimuthal vorticity equation, also know as the thermal wind balance or meridional circulation evolution equation (e.g., Brandenburg et al. 1992; Kitchatinov & Rüdiger 1995; Warnecke et al. 2013a), (21)where ∂/∂z = cosθ∂/∂r−r^-1sinθ∂/∂θ is the derivative along the rotation axis and is the mean vorticity. We neglect here the contribution of the Maxwell stress arising from the correlations of the fluctuating magnetic field and those from the mean magnetic field. The contributions of meridional flows turn out to be small and can be neglected. Warnecke et al. (2013a) could show that the baroclinic term, that is the second term in Eq. (21), is responsible for spoke-like differential rotation profile. Following the study of Warnecke et al. (2013a), we plot the two dominant terms on the right hand side of Eq. (21) together with the residual in Fig. 9 for Runs A1, A3t, A5, B1, and B3t. The meridional distributions of and rsinθ∂Ω²/∂z of Runs A1 and Runs A3t are very similar to Figs. 5 and 9 of Warnecke et al. (2013a), respectively. For most of the convection zone the two terms match well and ΔL is close to zero. In the upper 0.1 R, ΔL is non-zero. In the case of Run A1, ΔL is positive in the region above r = 0.94 R. However, right at the surface, it is negative, probably because of boundary effects. By contrast, Run B1 develops regions, where ΔL is positive (0.90 R ≤ r ≤ 0.94 R) and where ΔL is even more strongly negative (0.95 R ≤ r ≤ 0.99 R). In Runs A3t, A5, and B3t, the overall magnitude of terms in the thermal wind balance is lower, resulting in a significantly smaller residual ΔL. However, ΔL seems to be close to zero in all three runs and becomes non-zero only near the surface; see Fig. 9. Thus, there is no clear change in behavior depending on rotational influence or Prandtl numbers in the coronal runs. As the meridional circulation is stationary, apart from cycle dependent variations, the left-hand side of Eq. (21) is on average zero. Therefore, in the region close to the surface the contribution of the Reynolds stress, that is the third term on the rhs of Eq. (21), plays a more important role and will be related to the residual ΔL. We can rewrite this term (22)The latitudinal component of the divergence of is given by (23)and the radial component of the divergence is given by (24)where and ∂/∂φ = 0 because of the azimuthal mean. This means that the contributions related to , and are zero. Furthermore, we find that terms related to are too small to have a strong effect.

Fig. 9 The two dominant terms of Eq. (21) and their difference for Runs A1, A3t, A5, B1, and B3t in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (r = R).

Fig. 10 The three dominant contributions of the Reynolds stress (25), (26) and (27) as well as their sum for Runs A1, A3t, A5, B1, and B3t in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (r = R).

In Eq. (23) the first three terms and in Eq. (24) only the first term have a strong contribution to the thermal wind balance. We summarize them in the following three expressions: (25)where the first term on the lhs is the dominant one, (26)and (27)In Fig. 10, we plot all three terms and their sum for the same runs as in Fig. 9. In Run A1, their sum is positive, resulting in a negative contribution of the Reynolds stresses to the thermal wind balance; see Eq. (21). In Run B1, these terms are positive around r = 0.94 R and negative closer to the surface. is in all runs negative giving a positive contribution to the thermal wind balance. The runs with a coronal envelope seem to have a weaker and negative contribution from the Reynolds stresses. This is probably related to the low density stratification near the surface. If we compare the plots of Fig. 10 with Fig. 9, there seems to be some agreement of the residual ΔL with the contributions of the Reynolds stresses , , and . However, there seems to be a shift in radius, which might be related to still missing contributions, which we could not identify in the present work. This might be related to the fact that in our simulations, in particular near the surface, is not zero as assumed in many mean-field models (e.g., Kitchatinov & Rüdiger 1995, 2005). Furthermore, we have here neglected the contribution of the Maxwell stresses, which might also explain this inconsistency. However, we can confirm the result from Hotta et al. (2015) that Q_rθ is important near the surface and gives an important contribution to the thermal wind balance near the top boundary. However, this contribution is not directly related to the stress Q_rθ itself but to its radial gradient. Furthermore, also Q_θθ and Q_rr have a large contribution to the thermal wind balance in terms of and . Similar to Q_rθ in Fig. 8, changes sign in the runs where we have identified a near-surface shear layer; see Sect. 3.2. This effect is even stronger near the equator, where it could play a role in generating the near-surface shear layer.

3.5. Radial field boundary condition

Fig. 11 Inclination of the magnetic field near the surface plotted as a 2D histogram over radius r/R for Runs A1, A1c2, A3, and A3t. 0° means fully radial and 90° fully horizontal. The red line indicates the average for each radii.

In the case of Runs A1 and B1, we employ a radial field condition at the surface, r = R, whereas in the other cases the radial condition is enforced either at r = 1.01 R (Runs A1c, A1c2 and B1) or high above the convection zone (Runs A2, A3, A3t, A4, A5, B3, B3t). This affects the magnetic field distribution and therefore the dynamo in the convection zone. In the majority of the upper convection zone of Run A1, the angle between the field and the radial direction is distributed nearly uniformly, as indicated by a mean of the absolute angle being close to 60°. The field is fully radial just in the two uppermost grid points; see Fig. 11. In Runs A1c and A1c2 the field distribution is the same, except that the magnetic boundary condition is placed higher. For the runs with coronal envelope, the field is less radial at the surface than in Run A1. Furthermore, the field is more horizontal in the coronal envelope than in the convection zone. The comparison of Runs A3 and A3t suggests that the cooling layer, which reaches deeper in Run A3, allows the field to become less radial at and above the surface. This leads us to conclude that the radial boundary condition does not significantly affect the structure and inclination below r = 0.98 R.

Warnecke et al. (2016b) find turbulent radial downward pumping from a similar run, which can cause radial alignment of the magnetic field, but only close to the surface (r> 0.98 R) and at low latitudes, which is consistent with our findings. In local convection simulations, it has been found that turbulent pumping is a major effect causing horizontal magnetic field to move downward (e.g., Nordlund et al. 1992; Tobias et al. 1998; Ossendrijver et al. 2002). It is also believed that this effect can cause the field near the surface to become dominantly radial (Cameron et al. 2012) and is therefore the main reason for the success of surface flux transport models (see Mackay & Yeates 2012, for a review).

3.6. Cyclic dynamo solutions with dynamo wave propagation

Fig. 12 Mean toroidal magnetic field evolution as time-latitude (butterfly) diagram, plotted at a radius r = 0.98 R (top row) and time-radius diagram plotted at a 25° latitude (middle row) in kG during a 20 yr interval in the saturated stage for Runs A1, A1pc, A3, A3t, and B1. For Run B1, we plot the butterfly diagram at radius r = 0.84 R and its time interval is 40 yr. The dashed horizontal lines mark the equator (θ = π/ 2) and the radii r = R, r = 0.98 R and r = 0.85 R, respectively. Propagation of the mean magnetic field for the same runs (bottom row). Color coded is plotted during the saturated stage together with white arrows showing the direction of migration of an αΩ dynamo wave (Parker 1955; Yoshimura 1975); see Warnecke et al. (2014). We suppress the arrows above r = 1.05. The black solid lines indicate isocontours of at 2.0 kG. The dashed white lines indicate the surface (r = R).

In all runs of both sets, convective motions, overall rotation, and their interaction contribute to generating a large-scale magnetic field. Most of the runs produce cyclic magnetic fields in the saturated phase. In Fig. 12, we plot the time evolution of the mean toroidal magnetic field for a selection of runs. To investigate the cause of the propagation direction of the magnetic field, we apply the same technique as in Warnecke et al. (2014) and Käpylä et al. (2016a) and we calculate the propagation direction using the so-called Parker-Yoshimura rule (Parker 1955; Yoshimura 1975).

As discussed in Warnecke et al. (2014) in detail, Runs A1 and A1c2 show a solar-like equatorward migration of the mean field. This is probably caused by an equatorward propagating dynamo wave generated by the region of negative shear in the middle of convection zone; see Fig. 4 and bottom row of Fig. 12. We have confirmed this interpretation using all turbulent transport coefficients obtained with the test-field method (Warnecke et al. 2016b). To reiterate, Runs A1 and A1c2 have two dominant dynamo modes; equatorward and poleward migrating branches at mid-latitudes and high latitudes, respectively; see top row of Fig. 12. The dynamo wave seems to be formed in the middle of the convection zone and propagates toward the bottom and top of the convection zone; see bottom row of Fig. 12. The cycle period is around five years (Warnecke et al. 2014). Near the surface there is a weaker mode, which has a cycle period of 1.5 yr that propagates poleward at mid-latitudes and vanishes at higher latitudes. At the bottom of the convection zone exists another dynamo mode, see Fig. 12, which might cause a long-term variability. Such competing dynamo modes have been analyzed in detail in Käpylä et al. (2016a). Furthermore, the magnetic field in Run A1c shows no equatorward migration. As described by Warnecke et al. (2014), this is most likely caused by the much weaker negative radial gradient of Ω in the middle of convection zone, also visible in Fig. 5a and the suppression of turbulent convective motion near the surface, as discussed in Sect. 3.1. In this paper, we focus on the differences between the runs with and without a coronal envelope.

Before we discuss the runs with a coronal envelope, we investigate how the magnetic evolution changes, if we change the radial field condition to a perfect conductor condition at surface (Run A1pc). This is sufficient to cause the dynamo modes to change. Close to the surface the toroidal magnetic field becomes the strongest, reaching values of more than 10 kG which are more than two times higher than the maximum values in Runs A1 and A1c2; see top and middle row of Fig. 12. There the predicted migration direction is poleward, which agrees with the actual migration near the surface; see bottom row of Fig. 12. Furthermore, we find indication of a poleward migrating subdominant dynamo mode similar as in Runs A1, A1c, and A1c2. However, we find an equatorward propagating mode in the middle of the convection zone, which is most likely related to the smaller and weaker concentration of found in the middle of the convection zone which is predicted to propagate equatorward. It seems that the hydro-thermal setup tends to produce strong toroidal magnetic fields near the surface, but the radial field condition at the top boundary of Run A1, A1c, and A1c2 prevents this. The toroidal magnetic field is concentrated at mid-latitudes and produces nearly no field above 35° latitude. This is very surprising, given that the other runs without a coronal envelope produce strong magnetic fields at high latitudes. This absence of polar field in Run A1pc suggests a relation between the polar field and the radial field at the surface. The relation might be interesting in view of the Babcock-Leighton dynamo framework (e.g., Babcock 1961; Leighton 1964; Dikpati & Charbonneau 1999).

Adding a coronal envelope on top of the convection zone changes the magnetic field evolution significantly; see Fig. 12. For the runs of Set A, the mean toroidal magnetic field migrates mostly poleward in a region between the equator and ± 40° latitude. At high latitudes the magnetic field is weaker, but shows a tendency of equatorward migration, in particular in Run A3t this equatorward propagation is clearly visible. In the middle of convection zone, the equatorward propagation also dominates at lower latitudes. As a main difference compared to Run A1, the mean toroidal field occurs close to the surface instead of in bulk of the convection zone. At this location, the radial shear is positive causing a dynamo wave to propagate poleward, see bottom row of Fig. 12. This agrees with the magnetic field propagation. Furthermore, we find some indication of predicted equatorward migration in the middle of the convection zone, however this is not as clear as in Run A1. We do not find any clear difference between Runs A2, A3, and A4 (only Run A3 is shown). This implies that the size of the coronal envelope is not important for the magnetic field evolution. However, in Run A3t with a modified cooling profile the field is weaker in the band near the equator and more concentrated at or even above the surface in comparison with the other runs with coronal envelopes. Furthermore, the equatorward migration near the poles is more pronounced.

Also the migration periods of the large-scale fields change due to the coronal envelope. The cycle period of the poleward migrating field in Runs A2, A3, and A4 are around two years, which is shorter than in Runs A1 and A1c2, where it is around five years (Warnecke et al. 2014). However, the equatorward branch near the poles seems to appear only every second poleward cycle in Runs A2, A3, and A4, which give a similar period as in Runs A1 and A1c2. In Run A3t, the magnetic field near the equator does not show a regular behavior with a clear cycle. The equatorward branch near the poles has a period of around two years, The magnetic field evolution of Run A5 shows the same features as A3t: poleward migration or quasi-stationary behavior at low latitudes and equatorward migration at high latitudes.

In Set B with slower rotation, the magnetic field evolution shows a similar dependence on the cooling profile and the coronal envelope as in Set A. In Run B1 the magnetic field shows equatorward migration, in particular in the middle of the convection zone, as shown in Fig. 12. The radius-time diagram is similar to that of Run A1, but due to the slower rotation, the cycle period is extended to around ten years. This is consistent with the predicted dynamo wave scenario; see bottom row of Fig. 12. Analogously to Run A1c, the cooling layer in Run B1c causes the magnetic field to lose its equatorward migration mode. The other runs have mostly a stationary mode compared to their rapidly rotating counterparts, in particular at low latitudes. As in Set A, all migration directions of the magnetic field in runs of Set B agree with the predicted propagation.

In general it seems as if the toroidal magnetic field tends to be strong near the surface in the runs with an extended coronal envelope (Runs A1pc, A2, A3, A3t, A4, B3, and B3t), but in the runs with the vertical field boundary condition (Runs A1, A1c, A1c2, B1, and B1c) the boundary condition does not allow this. Overall, there is good agreement between the predicted and actual direction of migration of the mean toroidal field.

The magnetic field evolution is similar to that found in Warnecke et al. (2013a) whose runs have more than two times higher SGS Prandtl numbers and lower stratification (ρ/ρ_surf = 14), but the other parameters (Re, Pr_M, Co) are comparable with the runs of this work. In contrast to this work in the rapid rotating runs of Warnecke et al. (2013a) the field near the equator seems to continue to propagate equatorward. In their slower rotating case the field becomes quasi-stationary in the saturated stage also at high latitudes. At this point it is not clear if these differences are related to the larger SGS Prandtl numbers or to the weaker stratification.

4. Conclusion

In this work we have studied the influence of the upper boundary on convectively driven dynamos in spherical wedges. For this purpose we have added a convectively stable coronal envelope of different sizes on top of the convectively unstable dynamo region. This coronal envelope effectively corresponds to having a free boundary as opposed to a stress-free radial field boundary condition used in many earlier one-layer dynamo simulations. We confirm the result of Warnecke et al. (2013a) that the coronal envelope has an influence on the dynamo region leading to a change in differential rotation and magnetic field evolution. If the radius R_C of the coronal envelope is just 1% of the solar radius R, its influence is small and can entirely be related to small changes in the radial density and temperature profiles. If the size of the corona is larger (R_C ≥ 1.2 R), the influence is stronger. However, runs with coronal sizes extending higher than R_C = 1.2 R are nearly identical.

Regarding the hydrothermal properties, the influence of the corona can be summarized as follows: (i) the radial mass flux across the free surface does not show any major difference between runs with a small and an extended coronal envelope. The radial mass flux in our simulations is too small to have an influence on the flow properties inside the convection zone; (ii) the latitudinal temperature variations due to rotation are significantly weakened due to the presence of the coronal envelope; (iii) this seems to cause the differential rotation profile to become more spoke-like and weaker in runs with a coronal envelope; (iv) this effect can also be seen in the change of the off-diagonal Reynolds stress components due to the coronal envelope which can be explained by a change in the Λ effect and the anisotropy.

Furthermore, in the cases with and a coronal envelope, as well as in all runs with , we find the generation of a weak near-surface shear layer at low latitudes. We have related this generation to a change of sign in the meridional Reynolds stress tensor component Q_rθ near the surface in these runs. This component contributes to the meridional Λ effect and turns out to be non-zero in all simulations. Additionally, we have shown that the radial gradients of Q_rθ and Q_rr as well as the latitudinal gradient of Q_θθ are important near the surface to balance against differential rotation and the baroclinic term. A change of sign in these terms can be associated with the generation of a near-surface shear layer.

The coronal envelope serves as a free top boundary for the magnetic field. We find that the dynamo properties are generally strongly influenced by the choice of boundary conditions. This can be particularly important when the magnetic field is strong near the boundary, as was seen recently in connection with the latitudinal boundary condition for simple α² mean-field dynamos (Cole et al. 2016) in the wedge geometry that we consider here too. In all of our simulations a toroidal field tends to form preferentially near the surface of the convection zone but it is pushed down by a radial field boundary condition. However, with the radial field boundary condition, the field is only radial over a few grid points below the boundary; otherwise it is distributed for all simulations isotropically within the convection zone. As in Warnecke et al. (2014), we compare the migration of the mean toroidal magnetic field with the predicted propagation direction of the αΩ dynamo wave following the Parker-Yoshimura rule. It turns out that this rule can explain all the different migration directions found in our simulations. This is a remarkable result given the variety of the simulations.

We must emphasize that the combination of a dynamo region with a coronal envelope is still far from realistic. This is mostly because of the rather low density contrast and the strong viscous coupling between the two layers. We must therefore be aware of the possibility of artifacts, for example the occurrence of a differentially rotating coronal envelope at higher latitudes may be an example. This is not normally expected to be the case in the Sun (Timothy et al. 1975), although the observations have only been limited to regions where one sees a rigidly rotating coronal magnetic field. Indeed, recent work on coronal holes by Lionello et al. (2005) claims that the rigid rotation is only apparent.

Our work may have an impact on our understanding of the properties of the convection zone and the solar dynamo, because changing the properties of the boundary and studying its influence teaches us about the physical properties and dynamical effects within the solar convection zone. However, we also have to keep in mind that our model, similarly to other models of spherical convection for the Sun (e.g., Brun et al. 2004; Hotta et al. 2015), uses an SGS model for the unresolved convective heat flux. In addition, these models represent the strong radiative cooling at the surface just by a relaxation term, which ignores the effects of strong driving of convection by what is known as entropy rain (Spruit 1997; Brandenburg 2016). This can have dramatic effects on the typical length scales of stellar convection (Cossette & Rast 2016), which in turn would affect the differential rotation and dynamo properties.

In our simulations, we have confirmed the influence of the latitudinal temperature distribution on the differential rotation as described by mean-field models (e.g., Rüdiger 1989) and backed up by simulations (e.g., Miesch et al. 2006; Warnecke et al. 2013a). Only a small temperature difference between pole and equator is needed to cause the differential rotation to become more spoke-like. However, if the rotation causes too large latitudinal temperature variations, the differential rotation becomes dominated by the Taylor-Proudman theorem and has cylindrical contours. On the other hand, the strength of the differential rotation benefits from large latitudinal temperature variations. Furthermore, our work confirms the results of Hotta et al. (2015) in that Q_rθ may be important for the generation of the near-surface shear layer. However, we went a step further and identified the terms in the thermal wind balance which are important near the surface. These results appear to be in conflict with those of Kitchatinov & Rüdiger (1995, 2005), which explain the near-surface shear layer solely by the vanishing horizontal Λ effect near the surface. We confirm in some simulations that Λ_H goes to zero near the surface, but there we find no clear relation between Λ_H and a negative radial gradient of rotation near the surface. The associated Reynolds stress component also does not play a role in the thermal wind balance. Future helioseismic measurements may yield information about the presence and distribution of Q_rθ in the Sun.

The fact that the Parker-Yoshimura rule can explain the migration direction of mean toroidal magnetic field found in our simulations has an impact on the interpretation of equatorward migration of the magnetic field in the Sun. Applying this rule to the Sun will lead to the generation of toroidal field in the near-surface shear layer, where negative shear can cause equatorward migration. This would also imply that sunspots are formed near the surface by a local flux concentration mechanism as proposed by Brandenburg et al. (2011), Stein & Nordlund (2012), Warnecke et al. (2013b, 2016a), and Käpylä et al. (2016b). Furthermore, we find that in our simulations turbulent pumping is not strong enough to cause a preferred radial orientation of magnetic field near the surface as in the Sun. This is possibly due to insufficient stratification in our simulations. However, a detailed analysis of the effect of turbulent pumping crucially depends on the test-field method (e.g., Schrinner et al. 2007; Warnecke et al. 2016b).

Acknowledgments

The simulations have been carried out on supercomputers at GWDG, on the Max Planck supercomputer at RZG in Garching and in the facilities hosted by the CSC–IT Center for Science in Espoo, Finland, which are financed by the Finnish ministry of education. J.W. acknowledges funding by the Max-Planck/Princeton Center for Plasma Physics and funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement No. 623609. This work was partially supported by the Swedish Research Council grants No. 621-2011-5076 and 2012-5797 (A.B.), and the Academy of Finland Centre of Excellence ReSoLVE 272157 (M.J.K., P.J.K. and J.W.) and grants 136189, 140970, 272786 (P.J.K).

References

Appendix A: Meridional Λ effect

To obtain an expression for the Λ effect, we use the minimal tau approximation (Blackman & Field 2002, 2003); see also the derivation in Käpylä & Brandenburg (2008). We denote partial time derivatives by a dot and compute (A.1)in a non-rotating frame of reference (), using where the three dots denote nonlinear and gradient terms of u_r, u_θ, and u_φ/rsinθ that will be neglected. Below we also consider the case where gradient terms of u_θ will be included. In other words, the base state corresponds to rigid rotation with a meridional flow u_θ ∝ rsinθ, and gradients around this state are neglected.

Inserting Eqs. (A.2)–(A.4) into Eq. (A.1), we obtain an expression of the form (A.5)where L_ijk is a rank 3 tensor (related to the coefficients of the Λ effect) and R_ij stands for all remaining terms, in particular those that stem from the triple correlations. In the minimal tau approximation we replace those by a relaxation term with the turbulent correlation time τ, that is R_ij = −Q_ij/τ. Inserting this into Eq. (A.5) and assuming a steady state, that is , using that the background velocity correlation is of the form , we have where has been replaced by Ωrsinθ and gradient terms of u_θ/r are assumed to vanish, so that Eq. (A.3) yields . We note that, while Q_rφ and Q_θφ are proportional to Ω, the meridional component Q_rθ is proportional to (≡). If we were to allow u_θ/r to have non-vanishing gradients, we would have .

In that case, under isotropic conditions (), Λ_V = Λ_H = 0, but Λ_M is non-vanishing and would be positive. Therefore, the non-diffusive contribution to Q_rθ would be negative for a poleward flow. This is in agreement with the profiles of Λ_M and Q_rθ in the surface regions of our simulations; see Fig. 7.

Appendix B: The anisotropy parameter

The Λ effect is related to anisotropy parameters (Rüdiger 1980), defined as

Fig. B.1 Anisotropy parameters AM, AV, AH for Runs A1 (top row) and A3t (bottom).

All Tables

All Figures

	Fig. 1 Radial profiles of azimuthally and latitudinally averaged temperature ⟨ T ⟩ θφa) and density ⟨ ρ ⟩ θφb) normalized by their values at the bottom of the domain, T0 or ρ0, respectively for Set A. The inlays show the radial profiles near the surface.
In the text

Fig. 2 Radial profiles of azimuthally and latitudinally averaged turbulent rms velocity a) and ⟨ s ⟩ θφb) in m/s or normalized by cP respectively for Set A. Mean radial mass flux c) through the surface (r = R) for Set A. The horizontal line indicates the zero value and the three vertical thin lines indicate the equator (θ = 90°) and the intersection with the inner tangent cylinder (θ−90° ≈ ± 45° latitude).

In the text

	Fig. 3 Latitudinal dependence of at radius r = 0.98 Ra) and at radius r = 0.85 Rb) for Set A. Linestyles as in Fig. 1. The thin black lines indicate the zero value, the equator (θ = π/ 2) and the location of the inner tangent cylinder.
In the text

	Fig. 4 Angular velocity Ω(r,θ)/Ω0 for Runs A1 and A3 (top row), Runs A3t and A5 (middle row) and Runs B1 and B3t (bottom row). Here, is the local rotation rate. The dashed lines indicate the surface (r = R).
In the text

	Fig. 5 Differential rotation Ω(r,θ)/Ω0 as a function of radius for all runs of Set A at θ = 75° (15° latitude) a), and at the equator θ = π/ 2b), and Set B at mid-latitudes θ = 75°c), and at the equator θ = π/ 2d). The inlay in b) shows the angular velocity near the surface for Runs A2, A3, A3t, A4, and A5.
In the text

	Fig. 6 Meridional circulation in terms of the mass flux for Runs A1 and A3t. The dashed lines indicate the surface (r = R) and the red solid line the inner tangent cylinder.
In the text

	Fig. 7 From left to right: off-diagonal components of the Reynolds stress Qrθ, Qrφ and Qθφ normalized by νtΩ0 (first three columns), the turbulent viscosity in terms of molecular viscosity νt/ν (fourth column) and the three components of the Λ effect ΛM, ΛV, and ΛH (fifth to seventh column) normalized by νt for Runs A1 (top row), A3t (middle) and A5 (bottom).
In the text

	Fig. 8 Off-diagonal components of the Reynolds stress Qrθa), Qrφb) and Qθφc) normalized by νtΩ0, the turbulent viscosity in terms of molecular viscosity νt/νd) and the three components of the Λ effect ΛMe), ΛVf) and ΛHg) normalized by νt as well as the meridional flow h) for Set A in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (r = R).
In the text

	Fig. 9 The two dominant terms of Eq. (21) and their difference for Runs A1, A3t, A5, B1, and B3t in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (r = R).
In the text

	Fig. 10 The three dominant contributions of the Reynolds stress (25), (26) and (27) as well as their sum for Runs A1, A3t, A5, B1, and B3t in the northern hemisphere at 15° latitude. The thin black lines indicate the zero value and the surface (r = R).
In the text

	Fig. 11 Inclination of the magnetic field near the surface plotted as a 2D histogram over radius r/R for Runs A1, A1c2, A3, and A3t. 0° means fully radial and 90° fully horizontal. The red line indicates the average for each radii.
In the text

Fig. 12 Mean toroidal magnetic field evolution as time-latitude (butterfly) diagram, plotted at a radius r = 0.98 R (top row) and time-radius diagram plotted at a 25° latitude (middle row) in kG during a 20 yr interval in the saturated stage for Runs A1, A1pc, A3, A3t, and B1. For Run B1, we plot the butterfly diagram at radius r = 0.84 R and its time interval is 40 yr. The dashed horizontal lines mark the equator (θ = π/ 2) and the radii r = R, r = 0.98 R and r = 0.85 R, respectively. Propagation of the mean magnetic field for the same runs (bottom row). Color coded is plotted during the saturated stage together with white arrows showing the direction of migration of an αΩ dynamo wave (Parker 1955; Yoshimura 1975); see Warnecke et al. (2014). We suppress the arrows above r = 1.05. The black solid lines indicate isocontours of at 2.0 kG. The dashed white lines indicate the surface (r = R).

In the text

	Fig. B.1 Anisotropy parameters AM, AV, AH for Runs A1 (top row) and A3t (bottom).
In the text