Issue 
A&A
Volume 670, February 2023



Article Number  A101  
Number of page(s)  18  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202244990  
Published online  10 February 2023 
Threedimensional nonkinematic simulation of the postemergence evolution of bipolar magnetic regions and the BabcockLeighton dynamo of the Sun^{⋆}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: bekki@mps.mpg.de
Received:
16
September
2022
Accepted:
7
December
2022
Context. The BabcockLeighton fluxtransport model is a widely accepted dynamo model of the Sun that can explain many observational aspects of solar magnetic activity. This dynamo model has been extensively studied in a twodimensional (2D) meanfield framework in both kinematic and nonkinematic regimes. Recent threedimensional (3D) models have been restricted to the kinematic regime. In these models, the surface poloidal flux is produced by the emergence of bipolar magnetic regions (BMRs) that are tilted according to Joy’s law.
Aims. We investigate the prescription for emergence of a BMR in 3D nonkinematic simulations. In particular, we examine the effect of the radial extent of the BMR. We also report our initial results based on a cyclic BabcockLeighton dynamo simulation.
Methods. We extended a conventional 2D meanfield model of the BabcockLeighton fluxtransport dynamo into 3D nonkinematic regime, in which a full set of magnetohydrodynamic (MHD) equations are solved in a spherical shell using a YinYang grid. The largescale mean flows, such as differential rotation and meridional circulation, are not driven by rotationally constrained convection, but rather by the parameterized Λeffect in this model. For the induction equation, we used a BabcockLeighton αeffect source term by which the surface BMRs are produced in response to the dynamogenerated toroidal field inside the convection zone.
Results. We find that in the 3D nonkinematic regime, the tilt angle of a newlyemerged BMR is very sensitive to the prescription for the subsurface structure of the BMR (particularly, its radial extent). AntiJoy tilt angles are found unless the BMR is deeply embedded in the convection zone. We also find that the leading spot tends to become stronger (higher field strengths) than the following spot. The antiJoy’s law trend and the morphological asymmetry of the BMRs can be explained by the Coriolis force acting on the Lorentzforcedriven flows. Furthermore, we demonstrate that the solarlike magnetic cycles can be successfully obtained if Joy’s law is explicitly given in the BabcockLeighton αeffect. In these cyclic dynamo simulations, a strong Lorentz force feedback leads to cycle modulations in the differential rotation (torsional oscillation) and meridional circulation. The simulations, however, do not include radiative effects (e.g., enhanced cooling by faculae) that are required to properly model the torsional oscillations. The nonaxisymmetric components of the flows are found to exist as inertial modes such as the equatorial Rossby modes.
Key words: Sun: interior / Sun: activity / Sun: magnetic fields / Sun: helioseismology
Animation associated to Fig. 11 is available at https://www.aanda.org.
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model.
Open access funding provided by Max Planck Society.
1. Introduction
The Sun exhibits an 11yr cyclic magnetic activity that is sustained by the dynamo processes in the convection zone (e.g., Charbonneau 2020). The BabcockLeighton fluxtransport model is one of the most promising solar dynamo models available at present that can explain many observational features (e.g., Dikpati & Charbonneau 1999). In this model, the migration of sunspot groups toward the equator is attributed to the meridional flow near the base of the convection zone (Wang et al. 1991; Choudhuri et al. 1995). This model is supported by recent helioseismic observations in which the meridional flow is found to be poleward at the surface and equatorward at the base (Gizon et al. 2020). Another characteristic feature of this dynamo model is that the main process generating poloidal fields from toroidal fields is the socalled BabcockLeighton αeffect, whereby the surface poloidal fields are generated by the poleward advection and equatorial cancellation of the bipolar sunspots that are tilted with respect to the eastandwest directions (Babcock 1961; Leighton 1964).
The tendency that the leading spot is located closer to the equator than the following one has been dubbed Joy’s law (Hale et al. 1919). The physical origin of Joy’s law is still under debate: thin flux tube simulations have demonstrated that Joy’s law can be explained by the Coriolis force acting on the buoyantly rising flux tubes through the convection zone (D’Silva & Choudhuri 1993; Fan et al. 1994; Weber et al. 2011). On the other hand, recent observations have shown that the active regions emerge with eastwest alignment (with zero tilt), on average, and that the tilts characterized by Joy’s law are generated by the northsouth separation motions after emergence (Schunker et al. 2020).
Numerical investigations of the BabcockLeighton fluxtransport dynamo model have mainly been carried out in a twodimensional (2D) kinematic meanfield framework (e.g., Dikpati & Charbonneau 1999; Nandy & Choudhuri 2001; Chatterjee et al. 2004; Hazra et al. 2014; Karak & Cameron 2016). In these models, the BabcockLeighton αeffect is modeled as the axisymmetric poloidal source term which is localized near the surface. Although there are some nonkinematic studies where the dynamogenerated fields are allowed to give feedback on the mean flows (Rempel 2006; Ichimura & Yokoyama 2017; Inceoglu et al. 2017), the longitudinal component of the Lorentz force has been ignored because of the axisymmetry of the system.
Several recent studies have aimed to realize the BabcockLeighton process in a threedimensional (3D) fullspherical domain. Yeates & MuñozJaramillo (2013) first presented a kinematic model in which the upward velocity perturbation associated with the magnetic buoyancy is explicitly prescribed to produce the tilted bipolar magnetic regions (BMRs) at the surface. This method has also been used in Kumar et al. (2019) and Whitbread et al. (2019). Furthermore, Miesch & Dikpati (2014) developed a different model of the BabcockLeighton dynamo, in which the BMRs are artificially placed at the surface in response to the toroidal field at the base under the constraint of Joy’s law. In fact, this method is regarded as a 3D realization of the socalled “doublering” algorithm used in 2D meanfield models (Durney 1997; Nandy & Choudhuri 2001; MuñozJaramillo et al. 2010). The same model has also been used to study the longterm cycle variability (Karak & Miesch 2017). However, all of these models are restricted to kinematic regime. Therefore, it still remains unclear how the Lorentzforces of the BMRs affect their postemergence evolution and the resulting dynamo solution in the nonkinematic regime.
The models that have the least assumptions are provided by magnetohydrodynamic (MHD) convective dynamo simulations in a spherical shell (e.g., Brun et al. 2004; Ghizaru et al. 2010; Brown et al. 2010; Fan & Fang 2014; Hotta et al. 2016; Strugarek et al. 2017). However, they have difficulty in reproducing the largescale mean flows that we observe when the solar parameters are used (this problem is known as the convective conundrum, e.g., Nelson et al. 2018). Moreover, they still cannot capture the full dynamics of the fluxemergence and the resulting formation of BMRs at the surface comprehensively (Nelson et al. 2011; Fan & Fang 2014; Chen et al. 2017). Therefore, it is still helpful to use meanfield^{1} models in which the largescale meanflows are largely controlled with parameterizations of smallscale convective angular momentum transport (the Λeffect; see Kitchatinov & Ruediger 1995) and the flux emergence is modeled via a parametrization.
In this paper, we present a new numerical framework to study the BabcockLeighton dynamo processes of the Sun in a 3D nonkinematic regime, which takes advantage of both the meanfield approach for the solar largescale mean flows and the 3D realization of the BabcockLeighton process. Therefore, our model extends both the 2D nonkinematic meanfield models (e.g., Rempel 2006) and the 3D kinematic models (e.g., Miesch & Teweldebirhan 2016). Although our model is still less complete than 3D MHD convective dynamo models, we can solve the MHD dynamo equations under the constraints of the observed differential rotation and meridional circulation (which are difficult to obtain in 3D MHD convective dynamo models). We believe that our model can potentially provide many future applications, such as data assimilation and cycle prediction.
We note that a similar 3D meanfield model was recently presented by Pipin (2022), in which the meanfield induction equation is solved together with the meanfield hydrodynamic equations. Although the model considers the nonaxisymmetric magnetic field, such as newlyemerged BMRs, all the hydrodynamic variables were assumed to be axisymmetric and only the longitudinallyaveraged Lorentz force was taken into account. In our model, by contrast, we consider the nonaxisymmetric Lorentz force feedback on the nonaxisymmetric flows, which we find crucial with regard to the postemergence evolution of the BMRs.
Recently, various kinds of inertial modes have been discovered and identified on the Sun (e.g., Löptien et al. 2018; Gizon et al. 2021). Since these inertial modes have different mode properties from those of acoustic (p) modes, they are expected to be useful as an alternative tool to probe the interior of the Sun (Gizon et al. 2021; Bekki et al. 2022a). However, it remains largely uncertain how these modes are affected by the dynamogenerated magnetic fields. We believe that our 3D nonkinematic dynamo model can also be used to study the effects of magnetic fields on various inertial modes in the Sun’s convection zone in the nonlinear regime.
The organization of this paper is as follows. The numerical model is explained in detail in Sect. 2. In Sect. 3, we show how the postemergence evolution of the BMRs differs from the previous models. Our initial results of the cyclic dynamo are then presented in Sect. 4. We conclude by summarizing our results and discussing the future prospects in Sect. 5.
2. Model
2.1. Governing equations
We start by numerically solving a set of MHD equations in a spherical coordinate (r, θ, ϕ):
where g, ρ_{0}, p_{0}, and H_{p} denote the gravitational acceleration, density, pressure, and pressure scale height of the background state that is in an adiabaticallystratified hydrostatic equilibrium. We use the same radial profiles for the background stratification as the model presented in Rempel (2005) and Bekki & Yokoyama (2017). The quantities with subscript 1, ρ_{1}, and p_{1}, are the perturbations with respect to the background that are assumed to be sufficiently small, namely, p_{1}/p_{0}≈ρ_{1}/ρ_{0}≪1, so that the equation of state is linearized
where γ = 5/3 is the specific heat ratio and s_{1} is entropy perturbation from the adiabatic background. The rotation rate of the radiative core Ω_{0}/2π = 431.3 nHz is used for a system rotation rate.
The tensor Π represents the turbulent Reynolds stress associated with smallscale (subgridscale) convective motions that are not explicitly resolved in our model. This contains, in principle, the effects of turbulent diffusion and turbulent momentum transport (Λeffect, see Kitchatinov & Ruediger 1995). Therefore, the Reynolds stress is expressed as:
where S_{ik} and δ_{ik} denote the velocity deformation tensor and the Kroneckerdelta unit tensor. The detailed expression of S_{ik} in a spherical coordinate can be found in Fan & Fang (2014).
In our model, turbulent viscous, thermal, and magnetic diffusivities are all assumed to be isotropic. We use the same radial profiles for the viscous (ν_{vis}), thermal (κ), and magnetic (η) diffusivities as of Rempel (2006). The magnetic diffusivity η is 10^{12} cm^{2} s^{−1} at the top boundary and is on the order of 10^{11} cm^{2} s^{−1} in the bulk of the convection zone (see Fig. 2 in Rempel 2006). Although this diffusivity value is smaller than the estimate by the local mixinglength model (MuñozJaramillo et al. 2011), it is still about one order magnitude larger than the typical diffusivity value used in the kinematic fluxtransport models in the advectiondominated regime (e.g., Dikpati & Charbonneau 1999; Guerrero 2007).
In order to break the TaylorProudman’s constraint on the differential rotation via the thermal wind balance, a negative (positive) latitudinal entropy gradient in the northern (southern) hemisphere is required. Although there are several proposed mechanisms to generate this latitudinal entropy gradient (e.g., Kitchatinov & Ruediger 1995; Masada 2011; Hotta 2018), in this paper, we adopt the idea proposed by Rempel (2005) that the latitudinal entropy variation is generated by the radial meridional flows when the base of the convection zone is weakly subadiabatic. To this end, we give the superadiabaticity δ = ∇ − ∇_{ad}, with ∇ = dlnT/dlnp, as follows:
where 𝒯 denotes a transition function defined by:
We set the superadiabaticity at the poles and at the equator as δ_{pl} = −1.5 × 10^{−5} and δ_{eq} = −2 × 10^{−5}, respectively, at the base of the convection zone. The depths where the stratification changes from subadiabatic to adiabatic are given as r_{pl} = 0.725 R_{⊙} and r_{eq} = 0.735 R_{⊙}. The weakly subadiabatic layer near the base is thought to be an outcome of a nonlocal energy transport of strongly magnetized convection (Skaley & Stix 1991; Brandenburg 2016) and has been reported in some recent numerical simulations (Käpylä et al. 2017; Hotta 2017; Bekki et al. 2017). The subadiabaticity is slightly enhanced in the equatorial area owing to the latitudinal variation of the Coriolis force acting on lowentropy downdrafts (Karak et al. 2018).
The dimensionless tensor, Λ_{ik}, specifies the amplitude and direction of the turbulent momentum transport. In this model, we only consider the turbulent angular momentum transport. Therefore, we parameterize Λ_{rϕ}(=Λ_{ϕr}) and Λ_{θϕ}(=Λ_{ϕθ}) similarly to the model presented in Rempel (2005),
The overall amplitude of the Λeffect is given by Λ_{0} = 0.85. The inclination is set to λ = +(−)15° in the northern (southern) hemisphere. Thus, the associated angular momentum flux becomes mainly equatorward and is weakly directed away from the rotational axis. The spatial distribution of the Λeffect is specified as:
where d_{l} = 0.025 R_{⊙}. With this parameterization, the profiles of differential rotation and meridional circulation become similar to observations (Howe 2009; Gizon et al. 2020)
The quantities ζ_{r} and ζ_{θ} denote the random fluctuations due to the unresolved turbulent convection. In our model, the random fields ζ_{r} and ζ_{θ} are separately constructed by simply superposing multiple Gaussians as:
where the locations of Gaussian peaks (r_{i}, θ_{i}, ϕ_{i}) are randomly chosen and their amplitudes, c_{i}, are also randomly determined within the range of −2 < c_{i} < 2. The spatial scale of each gaussian is set as (δr, δθ, δϕ) = (0.03 R_{⊙}, 5° ,5° ). In our reference calculation, we set the number of Gaussians as N = 30. We generate the random field, ζ, at every time step and therefore it is uncorrelated in time. We note that the nonaxisymmetric flows can be partially driven by these random fluctuations of the Λeffect.
In the induction equation (Eq. (3)), we add an electromotiveforce ℰ to model the BabcockLeighton αeffect, by which the poloidal field is generated near the surface from the toroidal field near the base of the convection zone. We note that this term is only switched on when the cyclic dynamo is simulated in Sect. 4. A detail formulation of ℰ will be given in Sect. 4.1.
2.2. Numerical scheme
We numerically solved Eqs. (1)–(4) using a fourthorder centereddifference method for space and a fourstep RungeKutta scheme for time integration (Vögler et al. 2005). To avoid the severe CFL constraint for time step, we used the reduced speed of sound technique (Rempel 2005), so that the background sound speed is artificially reduced by a factor of ξ = 200, which still ensures that flows remain sufficiently subsonic (Hotta et al. 2014). Moreover, we used the hyperbolic divergence cleaning method (9wave method) for minimizing the numerical error resulting from the divergence of magnetic field (Dedner et al. 2002).
The numerical domain is a fullspherical shell extending from r_{min} = 0.65 R_{⊙} up to r_{max} = 0.985 R_{⊙}. The base of the convection zone is located at r_{bc} = 0.71 R_{⊙}. In order to avoid the singularities in a spherical coordinate at the poles, we used the YinYang grid (Kageyama & Sato 2004). For more details about the implementation of the YinYang grid, refer Bekki et al. (2022b). The grid resolution is 72(N_{r}) × 72(N_{θ}) × 216(N_{ϕ}) × 2 (Yin and Yang grids). The code is parallelized using message passing interface (MPI). At both radial boundaries, impenetrable and stressfree boundary condition is used for velocity. The magnetic field is assumed to be radial at the top and horizontal at the bottom.
We first carried out a hydrodynamic simulation until the largescale mean flows become quasistationary. Figure 1 shows profiles of the differential rotation, ⟨Ω⟩=Ω_{0} + ⟨v_{ϕ}⟩/(r sin θ), and meridional circulation, v_{m} = ⟨v_{r}⟩e_{r} + ⟨v_{θ}⟩e_{θ}, obtained from our hydrodynamic simulation. Here, ⟨ ⟩ denotes the longitudinal average. We then add magnetic fields to carry out MHD calculations.
Fig. 1. Results of the mean flows from the hydrodynamic calculation. (a) Differential rotation ⟨Ω⟩=Ω_{0} + ⟨v_{ϕ}⟩/(r sin θ). The black dotted curve shows the location of the base of the convection zone at r = 0.71 R_{⊙}, below which the stratification is subadiabatic. (b) Meridional circulation ⟨v_{θ}⟩. The black solid and dashed curves show contours of the streamfunction Ψ defined by ρ_{0}v_{m} = ∇ × (Ψe_{ϕ}) where v_{m} = ⟨v_{r}⟩e_{r} + ⟨v_{θ}⟩e_{θ}. The meridional circulation is counterclockwise (clockwise) in the northern (southern) hemisphere, i.e., the flow is poleward (equatorward) at the surface (base). 
3. Postemergence evolution of BMRs
In this section, we carry out a set of numerical simulations to study how the postemergence evolution of a BMR is dependent on how it is injected into the simulation. To this end, we solved the MHD Eqs. (1)–(4) starting from different initial magnetic field configurations for a newly emerged single BMR. For simplicity, we only considered the shortterm evolution of a BMR and do not discuss the longterm buildup of the polar fields and the resulting dynamo cycles. Hence, we set ℰ = 0.
3.1. Initial condition of the magnetic fields
For simplicity, we simulate the evolution of a BMR with zero initial tilt, that is, the leading and following polarity spots are perfectly eastwest aligned. The initial magnetic field is given as
where the vector potential, A_{ic}, is given by:
with
Here, we set the colatitude θ^{*} = 70° and ϕ^{*} = 0° and, thus, a BMR is located at the latitude of 20° in the northern hemisphere at the central meridian. The parameters Δr_{bmr}, Δθ_{bmr}, and Δϕ_{bmr} specify the radial, latitudinal, and longitudinal extent of the BMR. We fix Δθ_{bmr} = 5° but vary both Δr_{bmr} and Δϕ_{bmr} as free model parameters to change the subsurface shape of the BMR, as summarized in Table 1. In case 1, Δr_{bmr} is relatively small and Δϕ_{bmr} is large, indicating that the subsurface structure of the BMR is very shallow in radius but stretched in longitude. In case 3, on the other hand, the subsurface field morphology is changed to a verticallyelongated half ellipse. Case 2 is an intermediate case between case 1 and 3. In all cases, the amplitude of the initial field is determined such that the maximum radial field at the top boundary is 4 kG.
Model parameters for the postemergence BMRs simulations.
3.2. Temporal evolution
We first take a look at case 1 where the results are most drastically changed from the previous studies. Figure 2 shows the evolution of the radial field, flows at the surface, and the subsurface field of the BMR over the first several days after the emergence from case 1. As soon as the BMR is inserted, there emerge strong upflows at the surface because the mass is expelled from the flux tube due to the pressure imbalance. At the same time, the strong Lorentz force of the BMR drives strong longitudinal converging flows towards the polarity inversion line and latitudinal diverging flows along the polarity inversion line. These are clearly seen in Fig. 2b (middle panel). Consequently, the two spots (that are initially separated in longitude) are quickly pulled together and stretched in latitude, as shown in Fig. 2c (top panel). For comparison, we show the same snapshot from the corresponding kinematic calculation in Fig. 2d. The result reveals that the temporal evolution of the BMR is substantially changed in a nonkinematic regime, where the Lorentz force and the Coriolis force are taken into account. This type of evolution is not observed on the Sun.
Fig. 2. Temporal evolution of a BMR from case 1 at selected temporal points (a) t = 0 days, (b) t = 2.1 days, and (c) t = 6.4 days. The kinematic simulation with the same initial condition is shown in panel (d) at t = 6.4 days. Top panels show the radial field B_{r} at the surface r = 0.985 R_{⊙}. Thick black solid curves show the contours at B_{r} = 0.3 kG. The blue and red cross marks represent the locations of the fluxweighted center for the leading and following spots, and the grey straight lines are drawn to connect these two cross marks. Middle panels show the radial flows, v_{r}, near the surface, r = 0.98 R_{⊙}, by colored contour and the horizontal velocities (v_{θ}, v_{ϕ}) at the surface by vector arrows. Bottom panels show cross sections of the magnetic field strength B at the fixed latitude of 20°, denoted by black dashed lines in the top and middle panels. 
3.3. Tilt angle
In order to measure the tilt angle of the BMR, we compute the fluxweighted center locations of the leading and following polarity regions (θ_{L}, ϕ_{L}) and (θ_{F}, ϕ_{F}), respectively. The tilt angle, γ, is then defined as
Since we consider the BMR located in the northern hemisphere, Joy’s law is satisfied when γ > 0 by definition. It is seen from Fig. 2c (top) that the leading polarity spot is located at slightly higher latitude than the following polarity spot on average, indicating that the associated tilt angle is negative (γ ≈ −10° < 0 at t = 6.2 days). This is against Joy’s law. This is because the Coriolis force acts on the longitudinal converging flows (towards the polarity inversion line), as schematically illustrated in Fig. 3a. It should be emphasized that this generation of antiJoy’s law tilt is essentially 3D nonkinematic effect, and thus, it cannot be captured neither by the 2D nonkinematic models (which ignore the longitudinal dependence), nor the 3D kinematic models (which ignore the Lorentz force feedback).
Fig. 3. Schematic illustrations explaining the generation of antiJoy’s law tilt and the morphological asymmetry of the field strengths the BMRs. (a) Crosssection at the surface seen from the top. The red and blue arrows show the directions of the Lorentz force and the Coriolis force, respectively. (b) Crosssection at the fixed latitude seen from the equator to the north pole. (c) Threedimensional view of the evolution of a BMRs. 
Needless to say, the generation of negative tilts is contrary to the solar observations. One possible way of reconciling these simulations and observations is to assume that Joy’s law tilts are generated during the rise of the toroidal flux tubes (e.g., D’Silva & Choudhuri 1993) and thus are already embedded in the emerging BMRs, which overcomes the tendency to generate the antiJoy’s law tilts. In Sect. 4, we demonstrate that this is possible. Another possibility is that the BMRs emerge with nearly zero tilt but acquire the positive tilts after the emergence by asyetuncaptured physical process. For example, MartinBelda & Cameron (2016) proposed that the net positive tilt can be generated from the initial zerotilt state due to the coupled effects of differential rotation and active region inflows. This effect is not considered in this study since our code does not include the effect of radiative cooling in the BMRs, which geostrophically drives the active region inflows (e.g., Spruit 2003).
3.4. Morphological asymmetry
We also find a significant asymmetry between the leading and following spots: the leading polarity region tends to retain its compact shape and its strong field strength, whereas the following spot tends to gradually expand and the field becomes substantially weaker as time passes. To qualitatively assess this field strength asymmetry, we measured the maximum field strengths in the leading and following polarity regions B_{L} and B_{F}. In our simulation case 1, the leading spot has about twice stronger field than the following spot, B_{L}/B_{F}≈2 at t = 6.2 days. This morphological asymmetry of the BMRs has been a wellknown feature observed on the Sun (e.g., Bray & Loughhead 1979; Fisher et al. 2000) and often explained by the differential stretching of the rising Ωloop due to the Coriolis force (e.g., Fan et al. 1993). We refer to Fan (2021) for a more comprehensive review on the observational and theoretical studies on the morphological asymmetry.
Here, we provide a different but related explanation for this asymmetry. As illustrated in Fig. 3b, the two spots are attracted with each other by the Lorentz force and the Coriolis force acting on these longitudinal converging flows drives an downflow (upflow) inside the leading (following) polarity regions of the flux tube. This can be confirmed by the contour map of the radial motion in Fig. 2c (middle). Owing to the mass conservation, these up and downflows are accompanied by the horizontal converging and diverging motions, respectively, at the surface. Therefore, the leading spot becomes compressed and the field gets stronger, whereas the following spot becomes broader and the field gets weaker. This is schematically illustrated in Fig. 3c.
3.5. Dependence on the subsurface structure of a BMR
The results from the case 2 and case 3 are shown in Figs. 4 and 5, respectively. We can clearly see that the temporal evolution of a BMR is sensitively dependent upon the initial field structure of the BMR. From case 1 (where the BMR is localized in a shallow surface and the two spots are largely separated in longitude) to case 3 (where the BMR extends deeper in the convection zone and the longitudinal separation of the spots is small), the dominant component of the magnetic tension force is changed from horizontal (longitudinally converging at the surface) to vertical (upward at the bottom apex of the flux tube). Since the longitudinal component of the Lorentz force at the surface is dominated by the magnetic tension force, the surface driving of the longitudinal converging flows decreases from case 1 to case 3. This is explained in more detail in Appendix A. Consequently, the generation of a negative tilt is significantly suppressed from case 1 to case 3, as shown in Fig. 6a.
On the other hand, the field strength asymmetry between the leading and following spots can be found in all cases as shown in Fig. 6b. This is because, in case 3, the Lorentzforcedriven rising motion of the bottom apex of the flux tube drives the counterrotating flow inside the flux tube, which can enhance the converging (diverging) motion in the leading (following) polarity regions at the surface. Therefore, regardless of the subsurface shape of the BMR, the observed morphological asymmetry can be reproduced. In fact, the asymmetry is most significant in case 2, where the Coriolis force can act both on the longitudinally converging motion at the surface and on the radially upward motion of the deep flux tube. This is schematically illustrated in Figs. 3b and c.
Fig. 6. Temporal evolution of (a) the tilt angle γ and (b) the ratio of the maximum field strengths between the leading and following magnetic regions B_{L}/B_{F}. The Red, green, and blue curves correspond to cases 1, 2, and 3. 
Although the simulations reported in this section are based on a very simplified model of the halftorusshaped BMR with zero initial tilt, we find that its temporal evolution is extremely sensitive to its subsurface structure in the 3D nonkinematic regime. In particular, this study warns that the shallow BMRs model (which is conventionally used in many 2D nonkinematic models and 3D kinematic models) will lead to drastically different dynamo results when all these realistic 3D nonkinematic effects are included.
4. Cyclic dynamo with BabcockLeighton αeffect
In the previous section, we show that the newly emerged BMRs with shallow subsurface root have a general tendency to produce the antiJoy’s law tilts in the 3D nonkinematic regime. In this section, we demonstrate that despite this trend, the cyclic dynamo solution can be obtained if the tilt of Joy’s law is explicitly imposed for the emerging BMRs.
4.1. BabcockLeighton αeffect source term
Now, we switch on the electromotiveforce ℰ in Eq. (3) that represents the source of the BabcockLeighton αeffect, by which the surface poloidal field is produced as a result of the northsouth tilt of the BMRs (Babcock 1961; Leighton 1964). In our model, the emergence of BMRs at the surface is assumed to occur in response to the dynamogenerated toroidal field deep inside the convection zone: The tilted BMRs are instantaneously generated when the toroidal field near the base exceeds a threshold field strength. Our approach differs from the method presented in Yeates & MuñozJaramillo (2013), Kumar et al. (2019), and Pipin (2022), where the upward velocity associated with the magnetic buoyancy of the toroidal flux is prescribed. Rather, our method is similar to that used in Miesch & Dikpati (2014) and Miesch & Teweldebirhan (2016) where the BMRs are explicitly spotted at the surface.
We take the following steps to construct ℰ. First, the toroidal field near the base of the convection zone is computed at every time step:
where the average is taken over a narrow radial range near the base of the convection zone from r_{a} = 0.71 R_{⊙} to r_{b} = 0.735 R_{⊙}. Then, we determine the location of the flux emergence in a spherical surface (θ^{*}, ϕ^{*}). In order to suppress the emergence at high latitudes as suggested by observations, we apply a latitudinal mask to such that
where θ_{em} = 17.5° and Δθ_{tran} = 8.5°. We impose a necessary condition for the BMRs emergence to occur, that exceeds a threshold field strength B_{crit} = 750 G. When the above condition is satisfied on multiple points, the location of emergence (θ^{*}, ϕ^{*}) is randomly chosen between points satisfying the criterion.
Eventually, ℰ is expressed as follows being proportional to ,
where represents the spatial distribution of BMRs,
Here, Δr_{bmr}, Δθ_{bmr}, and Δϕ_{bmr} denote the radial, latitudinal, and longitudinal size of the BMRs. In order to demonstrate that the cyclic dynamo is possible even with the presence of strong nonkinematic effects (discussed in Sect. 3), we set Δr_{bmr} = 0.04 R_{⊙}. Therefore, the BMRs are confined in the shallow surface layer. We find that this is actually necessary to avoid the poleward dynamo wave propagation (see Appendix B.2). In the reference calculation, we set Δθ_{bmr} = Δϕ_{bmr} = 6°, which is consistent with observations suggesting the typical size of BMRs of r_{bmr} ≈ 5 − 100 Mm (e.g., Solanki 2003) that implies Δθ_{bmr} = Δϕ_{bmr} ≈ 2 r_{bmr}/R_{⊙} ≈ 0.4 − 8°. The overall amplitude of the BabcockLeighton αeffect is set to a_{0} = 100 km s^{−1}. This value, in combination with the typical toroidal field strength near the base kG (Dikpati & Charbonneau 1999), leads to the total magnetic flux of BMRs of 10^{22} − 10^{23} Mx, which is consistent with observations (Schrijver & Harvey 1994).
The northsouth tilt of the BMRs (ψ^{*}) obeys Joy’s law such that
where denotes the random fluctuation of the tilt angle (Hale et al. 1919; Howard 1991; Stenflo & Kosovichev 2012; Wang et al. 2015). For simplicity, we assume that the probability distribution of is roughly given by the following Gaussian distribution,
with σ_{f} = 15°. Unlike the kinematic model of Karak & Miesch (2017), a quenching term is not necessary in our model because the saturation of the dynamo occurs selfconsistently owing to the Lorentzforce feedback (Rempel 2006; Ichimura & Yokoyama 2017). Figure 7 shows examples of ℰ and the resulting tilted BMRs produced by our BabcockLeighton αeffect source where we assume sufficiently strong positive (negative) toroidal field near the base of the convection zone.
Fig. 7. Example structure of BMRs per each hemisphere produced from our BabcockLeighton αeffect model. Radial field at the solar surface is shown where red (blue) points represent positive (negative) B_{r}. The Solid black arrows denote the direction of the electromotiveforce, ℰ, defined by Eq. (22) with the appropriate Joy’s law. Positive (negative) toroidal field in the northern (southern) hemisphere is implicitly assumed near the base of the convection zone. 
In order to prevent overlapping emergence events on the same location in a very short time span, we introduce the following time delay algorithm as presented in Miesch & Dikpati (2014), Miesch & Teweldebirhan (2016), Karak & Miesch (2017): We used a cumulative lognormal distribution function of the emergence events defined as:
where Δ_{t} = t − t_{s} is the time lag since the last emergence event at t_{s}. A flux emergence event is allowed only when the cumulative 𝒞_{em} exceeds a number ∈[0, 1) randomly chosen at every time step. The standard deviation σ_{t} and the mean μ_{t} are specified by τ_{p} and τ_{s} as follows.
Here, we set τ_{p, 0} = 0.8 days, τ_{s, 0} = 1.9 days, Δτ_{p} = 0.75 days, and Δτ_{s} = 3.0 days. The quantity represents the horizontallyaveraged and B_{τ} = 1.5 kG denotes the threshold value of the toroidal field strength for characterizing phase in the activity cycle. Figure 8 shows the two examples of the cumulative 𝒞_{em}(Δ_{t}) each corresponding to the solar activity minima and maxima. Therefore, the flux emergence becomes more frequent during the activity maxima () and less frequent during the activity minima ().
Fig. 8. Cumulative lognormal distribution of the emergence events 𝒞_{em}(Δ_{t}) used in our model during the activity maxima (red) and activity minima (blue). 
We must note that our BabcockLeighton αeffect model is strongly spatiallylocalized and temporally intermittent. This is clearly different from the conventional 2D models with spatiallydistributed and temporallycontinuous source term (Choudhuri et al. 1995; Dikpati & Charbonneau 1999; Rempel 2006). Taking into account the tilt angle inclination of ℰ, the localization in longitudes, and the emergence frequency of the BMRs, we can estimate the corresponding α_{0} value within the 2D meanfield framework as
where we use the typical values ψ = 17.5°, Δ_{t} = 5 days, and Δt_{CFL} = 17 min. This α_{0} value is consistent with the previous 2D meanfield models.
4.2. Initial condition of magnetic fields
We add an axisymmetric dipolar field into the fullydeveloped hydrodynamic calculation shown in Fig. 1. The simulation is evolved until initial transients disappear and the dynamo cycles with quasisteady amplitudes are obtained. In the following subsections, we analyze the last three cycles of our simulation.
4.3. Dynamo cycles and the Lorentz force feedback
Figures 9a and b show the timelatitude plots of the longitudinallyaveraged radial field ⟨B_{r}⟩ at the surface and the toroidal field ⟨B_{ϕ}⟩ near the base of the convection zone, represented in terms of the wellknown magnetic butterfly diagram. We can clearly see the cyclic polarity reversals that occur roughly at every nine years, which is slightly shorter than the solar cycle period yet comparable. In each cycle, there is an equatorward migration of sunspot groups (BMRs) and the buildup of the polar field by poleward advection of the magnetic fluxes associated with the trailing sunspots. These are owing to the singlecell meridional circulation achieved in our model, which has an amplitude of about 15 m s^{−1} at the surface and 2 m s^{−1} near the base of the convection zone. The black solid lines in Fig. 9b denote the range of the emergence latitudes of BMRs at each time (the socalled active region belt). The phase of the equatorward advection of the toroidal field at the base corresponds to that of the emergence of the BMRs at the surface.
Fig. 9. Temporal evolution of the longitudinallyaveraged magnetic fields and horizontal velocities. (a) Azimuthal mean of the radial field ⟨B_{r}⟩ at the surface r = 0.985 R_{⊙}. (b) Azimuthal mean of the longitudinal field ⟨B_{ϕ}⟩ near the base of the convection zone r = 0.715 R_{⊙}. Black solid lines are the contours of the emerged BMRs at each time. (c) Torsional oscillation pattern δ⟨Ω⟩=⟨Ω⟩−⟨Ω⟩_{t} at the surface where ⟨⟩_{t} denotes the azimuthal and temporal average. (d) Azimuthal mean of the latitudinal velocity ⟨v_{θ}⟩ at the surface. Red (blue) in the northern hemisphere represents the equatorward (poleward) flow. (e) The same as (d) but near the base of the convection zone. Black dashed lines denote the contours of the toroidal field at the base (8.5 kG). (f) Entropy perturbation δ⟨s_{1}⟩=⟨s_{1}⟩−⟨s_{1}⟩_{t} at the surface. 
In our nonkinematic model, the dynamogenerated fields have strong impacts on flows via the Lorentz force feedback. Figure 9c shows the timelatitude plot of the fluctuation of the differential rotation δ⟨Ω⟩=⟨Ω⟩−⟨Ω⟩_{t}, where ⟨ ⟩_{t} denotes the longitudinal and temporal average. This is commonly known as torsional oscillations (Howard & Labonte 1980). We clearly find both poleward and equatorward propagating oscillation patterns with the typical amplitude of about 5 nHz at the surface. Figures 9d and e show the timelatitude plots of the latitudinal velocity ⟨v_{θ}⟩ at the top and bottom of the convection zone, respectively. Although the poleward flow at the surface and the equatorward flow near the base are strongly suppressed and disturbed during the activity maxima, the feedback is not large enough to switch off the advective transport of the magnetic fluxes (Ichimura & Yokoyama 2017).
Figure 9f shows the timelatitude plot of the entropy perturbation δ⟨s_{1}⟩=⟨s_{1}⟩−⟨s_{1}⟩_{t} at the surface with typical variation amplitude of about 250 erg g^{−1} K^{−1}, which corresponds to the temperature fluctuation of about 1.4 K. The positive entropy fluctuation can be seen along with the active region belt, implying that the surface is heated whenever the BMRs emerge due to the strong magnetic diffusion in our model. We note, however, that this is not likely in the real Sun: The surface is expected to be cooled by enhanced radiation in the BMRs, leading to lower temperature due to the radiative loss in the active region belt (Spruit 2003). Theories suggest that this radiative loss at the surface can produce the lowlatitude branches of the torsional oscillation by inducing the geostrophical flows around the BMRs (thermal forcing; Rempel 2006; Gizon & Rempel 2008). This effect is not included in our model. It should be emphasized that in our simulation reported here, the artificial heating in the active region belt may be responsible for the lowlatitude torsional oscillation branches due to the thermal forcing with the opposite sign.
Figure 10 shows the volumeintegrated kinetic and magnetic energies of various components. Their definitions are as follows:
Fig. 10. Temporal evolution of the volumeintegrated kinetic and magnetic energies. The red, green, and black lines denote the kinetic energies of the differential rotation KE_{DR}, the meridional circulation KE_{MC}, and the nonaxisymmetric flows KE_{m ≠ 0}. The blue, purple, and orange lines denote the magnetic energies of the mean toroidal field KE_{tor}, the mean poloidal field ME_{pol}, and the nonaxisymmetric fields ME_{m ≠ 0}. 
where the integrals are taken over the whole volume of the convection zone. The two largest energy reservoirs in our simulation are the differential rotation kinetic energy, KE_{DR}, and the toroidal field magnetic energy, ME_{tor}. When the toroidal field is amplified by the Ωeffect, KE_{DR} is converted to ME_{tor}. The toroidal field eventually becomes superequipartition (ME_{tor} > KE_{DR}) with respect to the differential rotation on average. In our reference simulation, this nonlinear Lorentz force feedback on differential rotation leads to a dynamo saturation. The nonaxisymmetric magnetic energy, ME_{m ≠ 0}, is greater than the mean poloidal field energy ME_{pol} because the BMRs are strongly nonaxisymmetric. The nonaxisymmetric fields drive strong nonaxisymmetric flows, whose kinetic energy, KE_{m ≠ 0}, is also greater than ME_{pol}. This implies that the nonaxisymmetric components of the magnetic fields and the flows are important for the convection zone dynamics.
Figure 11 shows snapshots of the magnetic fields and mean flows over the course of a magnetic cycle. The two leftmost panels show the mollweide projections of the radial magnetic field, B_{r}, at the surface and the toroidal field, B_{ϕ}, at the bottom convection zone. As prescribed in our BabcockLeighton source term, BMRs emerge at low latitudes obeying the Hale’s and Joy’s laws. Therefore, the radial magnetic field at the surface is substantially nonaxisymmetric. On the other hand, toroidal field near the base of the convection zone is found to be almost axisymmetric. Meridional plots on the third, fourth, and fifth columns of Fig. 11 show the azimuthally mean profiles of the poloidal and toroidal magnetic fields, differential rotation, and meridional circulation, respectively. When longitudinally averaged, the dynamo solution shows a qualitatively similar time evolution pattern as the previous 2D meanfield models (e.g., Rempel 2006; Ichimura & Yokoyama 2017), although our model has a stronger torsional oscillation and a stronger meridional circulation modulation. These presumably depend on the radial structure we have assumed for flux emergence.
Fig. 11. Time evolution of magnetic field and velocity. Shown are the snapshots at t = 12.9 yr (from (a) to (e)), t = 14.9 yr (from (f) to (j)), t = 17.9 yr (from (k) to (o)), t = 19.9 yr (from (p) to (t)) in Fig. 9. The mollweide projections in the first and second columns show the radial field B_{r} at the surface r = 0.985 R_{⊙} and longitudinal field B_{ϕ} near the base of the convection zone r = 0.715 R_{⊙}, respectively. The meridional plot in the third column represents the mean toroidal field (color scales) and poloidal field (contours). The meridional plots in the fourth and fifth columns represent the mean differential rotation and streamfunction of the meridional circulation, respectively. An animation of this figure is available online. 
4.4. Nonaxisymmetric flows
In our model, nonaxisymmetric flows are driven largely by the nonaxisymmetric Lorentz forces and partially by the random fluctuations in the Λeffect. Figure 12 shows a snapshot of the radial field at the surface (top rows) and the radial velocity near the surface (bottom rows). Black arrows represent the horizontal flow motions at the surface. Strong horizontal flows exist only in the vicinity of the BMRs: When a BMR emerge at the surface which happens instantaneously in our model, horizontal converging flows are driven towards the polarity inversion line with the typical amplitudes of about 100 m s^{−1}. This is owing to a strong magnetic tension force of the BMRs that pulls the two spots together. This strong converging flow further drives both horizontal outflows and radial downflows along the polarity inversion line, as shown in Figs. 12c and d. Due to these strong horizontal flows at the surface, a newlyemerged BMR that initially consists of two roundshaped sunspots is quickly deformed into an elongated shape along the polarity inversion line, as seen in Fig. 12a. This temporal evolution is similar to the simulation reported in Sect. 3.2 with horizontallyelongated BMRs (case 1). We find that this elongated feature of the BMRs, which is not observed on the Sun (van DrielGesztelyi & Green 2015), can be greatly suppressed if the BabcockLeighton αeffect is much weaker (Appendix B.1) or if BMRs have much deeper radial extent (Appendix B.2). It should also be noted that we do not include the radiative cooling associated with the active regions in our simulations. Thus, our model currently lacks the physics required to properly produce the observed inflows associated with active regions (Gizon et al. 2001; Spruit 2003).
Fig. 12. Snapshots of the radial field B_{r} at the surface (top panels) and the radial velocity v_{r} near the surface (bottom panels) at t = 10.3 yr in Fig. 9. The black arrows represent the horizontal flow (v_{θ}, v_{ϕ}) at the surface. Panels (c) and (d) are the zoomin of the panels (a) and (b), focusing on the single BMR, denoted by red thick solid lines. 
Other interesting nonaxisymmetric flow features are lowfrequency inertial modes of oscillation, in particular, the equatorial Rossby modes that have recently been detected on the Sun (e.g., Löptien et al. 2018). Figure 13 shows the equatorial power spectrum of latitudinal velocity v_{θ} near the surface from our nonkinematic dynamo simulation. We note that the spectrum is computed in a frame rotating at Ω_{ref}/2π = 431.3 nHz. We can clearly see the existence of the equatorial Rossby modes as represented by a clear power ridge along the expected dispersion relations (red points) for 3 ≤ m ≤ 12. In our simulation, these Rossby modes are excited both by the nonaxisymmetric random fluctuations in the Λeffect and by the nonaxisymmetric Lorentzforce; this is different from the rotating convection simulation of Bekki et al. (2022b), where they are stochastically excited by turbulent convective motions alone. It is implied that our code can be used to study the magnetic cycle dependence of the Rossby modes (or inertial modes in general) in the future.
Fig. 13. Equatorial power spectrum of latitudinal velocity v_{θ} near the surface r = 0.95 R_{⊙}. The spectra are computed in a frame rotating at Ω_{0}/2π = 431.3 nHz. The blue solid line represents the differential rotation rate at the surface ω = m(Ω_{sf} − Ω_{0}) where Ω_{sf} = ⟨Ω(0.95 R_{⊙}, π/2)⟩. The red points denote the theoretical dispersion relation of the sectoral Rossby modes, ω = −2Ω_{sf}/(m + 1)+m(Ω_{sf} − Ω_{0}). 
5. Summary and discussions
In this paper, we have developed a new BabcockLeighton fluxtransport dynamo model of the Sun. In our model, we do not solve the smallscale convection and focus on the largescale flows and magnetic structures in a full spherical shell. The solarlike largescale mean flows are driven by proper parameterization of the Λeffect. The model operates in a 3D nonkinematic regime, and therefore, is more realistic than the 2D nonkinematic models (e.g., Rempel 2006; Ichimura & Yokoyama 2017) and 3D kinematic models (e.g., Yeates & MuñozJaramillo 2013; Miesch & Dikpati 2014).
To better illustrate the major differences from the conventional 2D nonkinematic models and the 3D kinematic models, we first carried out a set of simulations for a single BMR with different initial subsurface structure. We find that when the BMR has a shallow subsurface structure and a large longitudinal separation, the postemergence evolution of the BMR becomes significantly changed from those from the conventional models: Even if the initial BMR is perfectly eastwest aligned (zero tilt angle), it begins to acquire a negative tilt angle (which is opposite to that of Joy’s law). The strength of the negative tilt angle decreases as the model bipole is embedded deeper in the solar convection zone. Furthermore, we find a strong asymmetry in the field strengths between the leading and following polarity regions. The leading polarity field becomes stronger whereas that of the following spot becomes weaker, which is similar to the observations. These results can be explained by the Coriolis force acting on the flows driven by the Lorentz force of the BMR (see Fig. 3).
We also carry out the cyclic solar dynamo simulation using the source term of the BabcockLeighton αeffect which is implemented in a 3D manner where Joy’s law tilts are explicitly given. We have successfully demonstrated that many observational features are reproduced in our model such as the activity cycles with decadal periods, the equatorward migration of the sunspot groups (BMRs), and the poleward transport of the surface radial fields. The nonlinear saturation of the dynamo occurs due to a strong Lorentz force feedback: The magnetic energy of the toroidal field amplified by the Ωeffect is found to exceed the kinetic energy of the differential rotation. This strong Lorentz force feedback can be seen in the cyclic modulations of the differential rotation (torsional oscillations) and the meridional circulation. We note, however, that our study does not exclude other nonlinear dynamo saturation mechanisms such as variability in the BabcockLeighton process (Weber et al. 2011; Karak & Miesch 2017) and magnetic quenching of the turbulent transport processes (Kitchatinov et al. 1994; Cattaneo & Hughes 1996; Yousef et al. 2003).
We note that our model is highly sensitive to various model parameters and that there are still several disagreements with the solar observations such as a slightly shorter cycle period of nine years, stronger radial field strengths at the surface of typical amplitudes of about 200 − 300 G, and slightly larger torsional oscillations. Obviously, the model parameters associated with the subsurface structure of the newlyemerged BMRs will be highly influential, as expected from the discussion in Sect. 3. The magnetic diffusivity, η, also has a substantial impact on the dynamo cycle properties by regulating the diffusive transport of magnetic fluxes in the Sun. The other important parameters would include: Λ_{0}, which determines the amplitudes of the differential rotation and meridional circulation, and a_{0}, which determines the field strengths of BMRs at the surface. A detailed parameter study is required in the future.
Due to the 3D nonkinematic nature of our model, we find the substantial nonaxisymmetric flows that are driven by the Lorentzforce of the BMRs. These flows have a spatial extent of about 10°, similar to what is seen in observations (Gizon et al. 2001; Löptien et al. 2017). Such flows are known to affect not only the evolution of the associated active regions but also the global magnetic field configuration through interaction with other nearby active regions (MartinBelda & Cameron 2016). This nonlinear interaction is sensitive to the details of the surface flows which are not yet consistent with the solar observations.
An important physical ingredient still missing in the present model is the enhanced radiative cooling associated with the BMRs, which will affect both the shortterm postemergence evolution of the BMRs and the longterm cyclic dynamo behaviors. This radiative loss will substantially affect the surface horizontal motions by geostrophycally inducing inflows around the active regions (e.g., Gizon et al. 2001; Gizon & Rempel 2008). These active region inflows are expected to affect the tilt angle (MartinBelda & Cameron 2016), regulate the poleward transport of the poloidal fluxes, and limit the buildup of the polar fields (Jiang et al. 2010), thus affecting the cycle amplitudes in the BabcockLeighton solar dynamo (Cameron & Schüssler 2012). Furthermore, it is often argued that the lowlatitude branches of the torsional oscillation are attributed to the thermallyinduced flows due to the enhanced surface cooling of the BMRs (thermal forcing; Spruit 2003; Rempel 2006, 2007). In the future model, we plan to include this effect to study how the postemergence of the BMRs and the nonlinear saturation of the dynamo are changed.
Lastly, we note that our code can be used to examine the impact of magnetic fields on various kinds of inertial modes that we found to exist in our simulation (see Fig. 13). Recent observations suggest that the amplitudes and frequencies of some of the solar equatorial Rossby modes exhibit a cycle dependence (Liang et al. 2019). If we properly understand the effects of deepseated magnetic fields on the mode frequencies and eigenfunctions of the equatorial Rossby modes, observations could potentially be used to infer the location and strength of the magnetic fields hidden in the Sun.
In the solar dynamo community, the “meanfield” models are conventionally regarded as 2D axisymmetric models where the mean is taken over longitudes. In this paper, however, we use the term “meanfield” in a more general sense: the mean should be regarded as an ensemble average or a spatial average over small portions in the convection zone that satisfy the Reynolds’ averaging rules. See Pipin (2022, Sect. 2.1 therein) for more details.
Acknowledgments
We thank an anonymous referee for constructive comments. We also thank B. Karak for helpful comments on the initial manuscript. Y.B. was enrolled in the International MaxPlanck Research School for Solar System Science at the University of Göttingen (IMPRS). Y.B. also acknowledges a support from a longterm scholarship program for degreeseeking graduate students abroad from the Japan Student Services Organization (JASSO). We acknowledge a support from ERC Synergy Grant WHOLE SUN 810218. All the numerical computations were performed at GWDG and the MaxPlanck supercomputer RZG in Garching.
References
 Babcock, H. W. 1961, ApJ, 133, 572 [Google Scholar]
 Bekki, Y., & Yokoyama, T. 2017, ApJ, 835, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Bekki, Y., Cameron, R. H., & Gizon, L. 2022a, A&A, 662, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bekki, Y., Cameron, R. H., & Gizon, L. 2022b, A&A, 666, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brandenburg, A. 2016, ApJ, 832, 6 [Google Scholar]
 Bray, R. J., & Loughhead, R. E. 1979, Sunspots (New York: Dover Publications) [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
 Cameron, R. H., & Schüssler, M. 2012, A&A, 548, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cattaneo, F., & Hughes, D. W. 1996, Phys. Rev. E, 54, R4532 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P. 2020, Liv. Rev. Sol. Phys., 17, 4 [Google Scholar]
 Chatterjee, P., Nandy, D., & Choudhuri, A. R. 2004, A&A, 427, 1019 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chen, F., Rempel, M., & Fan, Y. 2017, ApJ, 846, 149 [Google Scholar]
 Choudhuri, A. R., Schussler, M., & Dikpati, M. 1995, A&A, 303, L29 [NASA ADS] [Google Scholar]
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
 Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508 [NASA ADS] [CrossRef] [Google Scholar]
 D’Silva, S., & Choudhuri, A. R. 1993, A&A, 272, 621 [NASA ADS] [Google Scholar]
 Durney, B. R. 1997, ApJ, 486, 1065 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2021, Liv. Rev. Sol. Phys., 18, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., Fisher, G. H., & Deluca, E. E. 1993, ApJ, 405, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., Fisher, G. H., & McClymont, A. N. 1994, ApJ, 436, 907 [Google Scholar]
 Fisher, G. H., Fan, Y., Longcope, D. W., Linton, M. G., & Pevtsov, A. A. 2000, Sol. Phys., 192, 119 [CrossRef] [Google Scholar]
 Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Rempel, M. 2008, Sol. Phys., 251, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Duvall, T. L., Jr., & Larsen, R. M. 2001, in Recent Insights into the Physics of the Sun and Heliosphere: Highlights from SOHO and Other Space Missions, eds. P. Brekke, B. Fleck, & J. B. Gurman, 203, 189 [Google Scholar]
 Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020, Science, 368, 1469 [Google Scholar]
 Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guerrero, G., & de Gouveia Dal Pino, E. M., 2007, A&A, 464, 341 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hale, G. E., Ellerman, F., Nicholson, S. B., & Joy, A. H. 1919, ApJ, 49, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Hazra, G., Karak, B. B., & Choudhuri, A. R. 2014, ApJ, 782, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
 Hotta, H. 2018, ApJ, 860, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2014, ApJ, 786, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2016, Science, 351, 1427 [NASA ADS] [CrossRef] [Google Scholar]
 Howard, R. F. 1991, Sol. Phys., 136, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Howard, R., & Labonte, B. J. 1980, ApJ, 239, L33 [Google Scholar]
 Howe, R. 2009, Liv. Rev. Sol. Phys., 6, 1 [Google Scholar]
 Ichimura, C., & Yokoyama, T. 2017, ApJ, 839, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Inceoglu, F., Arlt, R., & Rempel, M. 2017, ApJ, 848, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, J., Işik, E., Cameron, R. H., Schmitt, D., & Schüssler, M. 2010, ApJ, 717, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Kageyama, A., & Sato, T. 2004, Geochem. Geophys. Geosyst., 5, Q09005 [Google Scholar]
 Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23 [CrossRef] [Google Scholar]
 Karak, B. B., & Cameron, R. 2016, ApJ, 832, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., & Miesch, M. 2017, ApJ, 847, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., Miesch, M., & Bekki, Y. 2018, Phys. Fluids, 30, 046602 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Ruediger, G. 1995, A&A, 299, 446 [NASA ADS] [Google Scholar]
 Kitchatinov, L. L., Ruediger, G., & Kueker, M. 1994, A&A, 292, 125 [NASA ADS] [Google Scholar]
 Kumar, R., Jouve, L., & Nandy, D. 2019, A&A, 623, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leighton, R. B. 1964, ApJ, 140, 1547 [Google Scholar]
 Liang, Z.C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löptien, B., Birch, A. C., Duvall, T. L., et al. 2017, A&A, 606, A28 [Google Scholar]
 Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
 MartinBelda, D., & Cameron, R. H. 2016, A&A, 586, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masada, Y. 2011, MNRAS, 411, L26 [CrossRef] [Google Scholar]
 Miesch, M. S., & Dikpati, M. 2014, ApJ, 785, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Teweldebirhan, K. 2016, Adv. Space Res., 58, 1571 [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt, H. K. 1978, Magnetic Field Generation in Electrically Conducting Fluids (Cambridge: Cambridge University Press) [Google Scholar]
 MuñozJaramillo, A., Nandy, D., Martens, P. C. H., & Yeates, A. R. 2010, ApJ, 720, L20 [CrossRef] [Google Scholar]
 MuñozJaramillo, A., Nandy, D., & Martens, P. C. H. 2011, ApJ, 727, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Nandy, D., & Choudhuri, A. R. 2001, ApJ, 551, 576 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 739, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, N. J., Featherstone, N. A., Miesch, M. S., & Toomre, J. 2018, ApJ, 859, 117 [Google Scholar]
 Parker, E. N. 1955, ApJ, 122, 293 [Google Scholar]
 Pipin, V. V. 2022, MNRAS, 514, 1522 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2006, ApJ, 647, 662 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2007, ApJ, 655, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., & Harvey, K. L. 1994, Sol. Phys., 150, 1 [Google Scholar]
 Schunker, H., Baumgartner, C., Birch, A. C., et al. 2020, A&A, 640, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Skaley, D., & Stix, M. 1991, A&A, 241, 227 [NASA ADS] [Google Scholar]
 Solanki, S. K. 2003, A&A Rev, 11, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C. 2003, Sol. Phys., 213, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Stenflo, J. O., & Kosovichev, A. G. 2012, ApJ, 745, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Strugarek, A., Beaudoin, P., Charbonneau, P., Brun, A. S., & do Nascimento, J. D., 2017, Science, 357, 185 [NASA ADS] [CrossRef] [Google Scholar]
 van DrielGesztelyi, L., & Green, L. M. 2015, Liv. Rev. Sol. Phys., 12, 1 [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [Google Scholar]
 Wang, Y. M., Sheeley, N. R., Jr., & Nash, A. G. 1991, ApJ, 383, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y. M., Colaninno, R. C., Baranyi, T., & Li, J. 2015, ApJ, 798, 50 [Google Scholar]
 Weber, M. A., Fan, Y., & Miesch, M. S. 2011, ApJ, 741, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Whitbread, T., Yeates, A. R., & MuñozJaramillo, A. 2019, A&A, 627, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yeates, A. R., & MuñozJaramillo, A. 2013, MNRAS, 436, 3366 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshimura, H. 1975, ApJ, 201, 740 [Google Scholar]
 Yousef, T. A., Brandenburg, A., & Rüdiger, G. 2003, A&A, 411, 321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Lorentz force of a model BMR
In Sect. 3, we show that the postemergence evolution of a model BMR is very sensitive to its subsurface shape and this is because of the change in the Lorentz force of the BMR. In this appendix, we show the radial and longitudinal components of the Lorentz force associated with the model BMR discussed in Sect. 3. The initial magnetic field is specified by Eqs. (16)–(18). We decompose the Lorentz force per unit mass into the magnetic pressure force (F^{MP}) and the magnetic tension force (F^{MT}) as:
Figures A.1a–d show the longitudinal component of the magnetic pressure and magnetic tension forces for case 1 and case 3, respectively. In case 1 where the BMR is localized in a shallow surface layer and the two spots are distantly separated in longitude, the longitudinal Lorentz force at the surface is dominated by the magnetic tension force , which acts to pull the two spots together. On the other hand, in case 3, where the BMR is anchored deeper in the convection zone and the two spots are located close by each other, the longitudinal component of the magnetic tension force significantly decreases because the field lines are no longer curved at the surface. Figure A.1e shows the ratio between and at the surface, implying that the longitudinal Lorentz force is dominated by the magnetic tension force when Δr_{bmr} is small, that is, Δr_{bmr} ≲ 0.05 R_{⊙}. Within this parameter regime, the BMR is expected to gain a negative (antiJoy’s law) tilt owing to the strong tensionforcedriven longitudinal converging flows and the subsequent Coriolis force acting on them.
Fig. A.1. Longitudinal component of the magnetic pressure and magnetic tension forces, and . (a,b) and for the initial BMR for the case 1, respectively. The crosssections at θ = θ^{*} are shown. (c,d) The same as panels (a) and (b) for the case 3. (e) Ratio between the maximum amplitudes of and at the surface (r = 0.985 R_{⊙}) with different combinations of the model parameters, Δr_{bmr} and Δϕ_{bmr}. The parameters used in cases 1–3 are denoted by red, green, and blue circles. The white dotted line represents the contour line of unity. 
The radial component of the magnetic pressure and tension forces are shown in Fig. A.2. It can be seen that the radial Lorentz force is largely dominated by the magnetic pressure force in case 1, by which the plasma inside the flux tube is pushed outward (both radially downward and upward). The associated pressure disturbances are propagated as (magneto) acoustic waves, which can be seen in Figs. 2a–4a. This initial relaxation occurs on timescales shorter than the dynamical timescale of the BMR evolution. On the other hand, in case 3, the radial Lorentz force is dominated by the magnetic tension force which is directing upward. This upward motion, coupled with the Coriolis effect, is expected to induce the retrograde plasma flows inside the flux tube and produce the field asymmetry between the leading and the following spots in cases 2 and 3, as illustrated in Fig. 3b.
Fig. A.2. Radial component of the magnetic pressure and magnetic tension forces and shown in the same way as in Fig. A.1. (e) Ratio between the maximum amplitudes of and at the polarityinversion line (ϕ = 0). 
Appendix B: Cyclic dynamo simulations with different model parameters
In Sect. 4, we describe the results from our reference simulation. In this appendix, we report two additional cyclic dynamo simulations with different model parameters.
B.1. a_{0} = 50 km s^{−1}
We carried out the simulation with weaker BabcockLeighton αeffect source, in which a_{0} is decreased to 50 km s^{−1} while all the other parameters remain unchanged from the reference calculation reported in Sect. 4. Figure B.1a shows the magnetic butterfly diagram at the surface. The amplitude of the longitudinallyaveraged radial field ⟨B_{r}⟩ can be reduced by a factor of 3 (see Fig. 9a). It is seen that the simulation can still reproduce many observed properties of the solar dynamo such as the decadal polarity reversals and the equatorward migration of the activity belt. We find that the Lorentz force feedback is substantially weaker in this simulation. Figure B.1b shows the volumeintegrated kinetic and magnetic energies as functions of time. It is shown the kinetic energy of the differential rotation, KE_{DR}, is always greater than the toroidal magnetic energy, ME_{tor}, by about one order of magnitude, suggesting that the Ωeffect operates as if in the kinematic regime without being quenched by Lorentz force feedback. This is in clear contrast to the reference case shown in Fig. 10. In fact, the typical torsional oscillation amplitude is found to be about 1 − 2 nHz, which is smaller than those of our reference model and the solar observations.
Fig. B.1. Temporal evolution of the cyclic dynamo simulation with weak BabcockLeighton αeffect (a_{0} = 50 km s^{−1}). (a) Magnetic butterfly diagram at the surface r = 0.985 R_{⊙}. (b) Temporal evolution of the volumeintegrated kinetic and magnetic energies. The same definition of colors as in Fig. 10 is used. 
Figure B.2 shows the temporal snapshot of the radial magnetic field, B_{r}, and the nonaxisymmetric flows at the surface. The latitudinal elongation of the BMRs (which is commonly seen in our reference calculation but not observed on the Sun) can be significantly alleviated.
Fig. B.2. Snapshot of the radial field B_{r} at the surface (top panels) and the radial velocity v_{r} near the surface (bottom panels) for the dynamo simulation with weak BabcockLeighton αeffect source term (a_{0} = 50 km s^{−1}). The same notation as Fig. 12 is used. 
We find that in the simulation with weak BabcockLeighton αeffect source, some inertial modes exhibit a very strong dependence on the magnetic activity cycles. This, as well as the other dynamo properties of this simulation, will be reported and discussed in great detail in a future publication (Bekki et al., in prep.).
B.2. Δr_{bmr} = 0.12 R_{⊙}
We further carried out the additional simulation where the newlyemerged BMRs have a much deeper radial extent by setting Δr_{bmr} = 0.12 R_{⊙}. All the other model parameters are unchanged from the reference simulation. Figures B.3a and b show the timelatitude plots of the longitudinally averaged radial magnetic field ⟨B_{r}⟩ at the surface and the toroidal field ⟨B_{ϕ}⟩ at the base of the convection zone, respectively. It is shown that when the BMRs are anchored deep in the convection zone, the emergence latitude of the BMRs tends to propagate poleward in time, in striking contrast to the solar observations. The cycle period is about six or seven years, which is shorter than that of the reference simulation and in the observation. The torsional oscillation δ⟨Ω⟩ at the surface also reflects this poleward migration of activity belt, forming clear poleward highlatitude branches, as shown in Fig. B.3c. We confirm that the longitudinallyaveraged meridional flow, ⟨v_{θ}⟩, is always equatorward near the base of the convection zone and does not flip the sign via the Lorentz force during the activity cycles (Fig. B.3d). Therefore, the poleward migration of the activity belt cannot be attributed to the advection by the poleward meridional flow.
Fig. B.3. Temporal evolution of the longitudinallyaveraged magnetic fields and horizontal velocities for the dynamo simulation with deeper BabcockLeighton αeffect (Δr_{bmr} = 0.12 R_{⊙}). (a) Azimuthallyaveraged radial field ⟨B_{r}⟩ at the surface r = 0.985 R_{⊙}. (b) Azimuthallyaveraged toroidal field ⟨B_{ϕ}⟩ near the base of the convection zone r = 0.715 R_{⊙}. Black solid lines are the contours of the emerged BMRs at each time. (c) Torsional oscillation pattern δ⟨Ω⟩ at the surface. (d) Azimuthallyaveraged latitudinal velocity ⟨v_{θ}⟩ near the base of the convection zone. 
We attribute the origin of the poleward migration of activity belt to the αΩ dynamo waves. When Δr_{bmr} is large, the mean poloidal field generation by the BabcockLeighton αeffect (which occurs in response to the toroidal field at the base of the convection zone) is no longer confined in a thin surface layer. If the locations of the αeffect and the Ωeffect are not spatially separated, the dynamo waves cannot be avoided. In the simulation reported here, the BabcockLeighton αeffect has an effective α_{ϕϕ} which is positive (negative) in the northern (southern) hemisphere. According to the ParkerYoshimura sign rule (Parker 1955; Yoshimura 1975), the propagation direction of the αΩ dynamo wave s_{αΩ} is given by:
Therefore, the dynamo waves propagate poleward in low to middle latitudes where ∂⟨Ω⟩/dr is positive in the convection zone. We must note that the αΩ dynamo waves seen in our simulation are distinct from the conventional ones where the αeffect represents the smallscale helical turbulence (e.g., Parker 1955; Moffatt 1978).
Figure B.4 shows the temporal snapshot of the radial magnetic field B_{r} and the nonaxisymmetric flows at the surface from the simulation with deep BabcockLeighton αeffect (Δr_{bmr} = 0.12 R_{⊙}). We can see that the latitudinal elongation of the BMRs which is characteristic in the model with shallow BMRs (see Fig. 12) can be significantly relaxed.
Fig. B.4. Snapshot of the radial field B_{r} at the surface (top panels) and the radial velocity v_{r} near the surface (bottom panels) for the dynamo simulation with deeper BabcockLeighton αeffect (Δr_{bmr} = 0.12 R_{⊙}). The same notation as Fig. 12 is used. 
Movie
Movie 1 associated with Fig. 11 (animation_cyclic_dynamo) Access here
All Tables
All Figures
Fig. 1. Results of the mean flows from the hydrodynamic calculation. (a) Differential rotation ⟨Ω⟩=Ω_{0} + ⟨v_{ϕ}⟩/(r sin θ). The black dotted curve shows the location of the base of the convection zone at r = 0.71 R_{⊙}, below which the stratification is subadiabatic. (b) Meridional circulation ⟨v_{θ}⟩. The black solid and dashed curves show contours of the streamfunction Ψ defined by ρ_{0}v_{m} = ∇ × (Ψe_{ϕ}) where v_{m} = ⟨v_{r}⟩e_{r} + ⟨v_{θ}⟩e_{θ}. The meridional circulation is counterclockwise (clockwise) in the northern (southern) hemisphere, i.e., the flow is poleward (equatorward) at the surface (base). 

In the text 
Fig. 2. Temporal evolution of a BMR from case 1 at selected temporal points (a) t = 0 days, (b) t = 2.1 days, and (c) t = 6.4 days. The kinematic simulation with the same initial condition is shown in panel (d) at t = 6.4 days. Top panels show the radial field B_{r} at the surface r = 0.985 R_{⊙}. Thick black solid curves show the contours at B_{r} = 0.3 kG. The blue and red cross marks represent the locations of the fluxweighted center for the leading and following spots, and the grey straight lines are drawn to connect these two cross marks. Middle panels show the radial flows, v_{r}, near the surface, r = 0.98 R_{⊙}, by colored contour and the horizontal velocities (v_{θ}, v_{ϕ}) at the surface by vector arrows. Bottom panels show cross sections of the magnetic field strength B at the fixed latitude of 20°, denoted by black dashed lines in the top and middle panels. 

In the text 
Fig. 3. Schematic illustrations explaining the generation of antiJoy’s law tilt and the morphological asymmetry of the field strengths the BMRs. (a) Crosssection at the surface seen from the top. The red and blue arrows show the directions of the Lorentz force and the Coriolis force, respectively. (b) Crosssection at the fixed latitude seen from the equator to the north pole. (c) Threedimensional view of the evolution of a BMRs. 

In the text 
Fig. 4. Temporal evolution of a BMR for the simulation case 2. Details are the same as in Fig. 2. 

In the text 
Fig. 5. Temporal evolution of a BMR for the simulation case 3. Details are the same as in Fig. 2. 

In the text 
Fig. 6. Temporal evolution of (a) the tilt angle γ and (b) the ratio of the maximum field strengths between the leading and following magnetic regions B_{L}/B_{F}. The Red, green, and blue curves correspond to cases 1, 2, and 3. 

In the text 
Fig. 7. Example structure of BMRs per each hemisphere produced from our BabcockLeighton αeffect model. Radial field at the solar surface is shown where red (blue) points represent positive (negative) B_{r}. The Solid black arrows denote the direction of the electromotiveforce, ℰ, defined by Eq. (22) with the appropriate Joy’s law. Positive (negative) toroidal field in the northern (southern) hemisphere is implicitly assumed near the base of the convection zone. 

In the text 
Fig. 8. Cumulative lognormal distribution of the emergence events 𝒞_{em}(Δ_{t}) used in our model during the activity maxima (red) and activity minima (blue). 

In the text 
Fig. 9. Temporal evolution of the longitudinallyaveraged magnetic fields and horizontal velocities. (a) Azimuthal mean of the radial field ⟨B_{r}⟩ at the surface r = 0.985 R_{⊙}. (b) Azimuthal mean of the longitudinal field ⟨B_{ϕ}⟩ near the base of the convection zone r = 0.715 R_{⊙}. Black solid lines are the contours of the emerged BMRs at each time. (c) Torsional oscillation pattern δ⟨Ω⟩=⟨Ω⟩−⟨Ω⟩_{t} at the surface where ⟨⟩_{t} denotes the azimuthal and temporal average. (d) Azimuthal mean of the latitudinal velocity ⟨v_{θ}⟩ at the surface. Red (blue) in the northern hemisphere represents the equatorward (poleward) flow. (e) The same as (d) but near the base of the convection zone. Black dashed lines denote the contours of the toroidal field at the base (8.5 kG). (f) Entropy perturbation δ⟨s_{1}⟩=⟨s_{1}⟩−⟨s_{1}⟩_{t} at the surface. 

In the text 
Fig. 10. Temporal evolution of the volumeintegrated kinetic and magnetic energies. The red, green, and black lines denote the kinetic energies of the differential rotation KE_{DR}, the meridional circulation KE_{MC}, and the nonaxisymmetric flows KE_{m ≠ 0}. The blue, purple, and orange lines denote the magnetic energies of the mean toroidal field KE_{tor}, the mean poloidal field ME_{pol}, and the nonaxisymmetric fields ME_{m ≠ 0}. 

In the text 
Fig. 11. Time evolution of magnetic field and velocity. Shown are the snapshots at t = 12.9 yr (from (a) to (e)), t = 14.9 yr (from (f) to (j)), t = 17.9 yr (from (k) to (o)), t = 19.9 yr (from (p) to (t)) in Fig. 9. The mollweide projections in the first and second columns show the radial field B_{r} at the surface r = 0.985 R_{⊙} and longitudinal field B_{ϕ} near the base of the convection zone r = 0.715 R_{⊙}, respectively. The meridional plot in the third column represents the mean toroidal field (color scales) and poloidal field (contours). The meridional plots in the fourth and fifth columns represent the mean differential rotation and streamfunction of the meridional circulation, respectively. An animation of this figure is available online. 

In the text 
Fig. 12. Snapshots of the radial field B_{r} at the surface (top panels) and the radial velocity v_{r} near the surface (bottom panels) at t = 10.3 yr in Fig. 9. The black arrows represent the horizontal flow (v_{θ}, v_{ϕ}) at the surface. Panels (c) and (d) are the zoomin of the panels (a) and (b), focusing on the single BMR, denoted by red thick solid lines. 

In the text 
Fig. 13. Equatorial power spectrum of latitudinal velocity v_{θ} near the surface r = 0.95 R_{⊙}. The spectra are computed in a frame rotating at Ω_{0}/2π = 431.3 nHz. The blue solid line represents the differential rotation rate at the surface ω = m(Ω_{sf} − Ω_{0}) where Ω_{sf} = ⟨Ω(0.95 R_{⊙}, π/2)⟩. The red points denote the theoretical dispersion relation of the sectoral Rossby modes, ω = −2Ω_{sf}/(m + 1)+m(Ω_{sf} − Ω_{0}). 

In the text 
Fig. A.1. Longitudinal component of the magnetic pressure and magnetic tension forces, and . (a,b) and for the initial BMR for the case 1, respectively. The crosssections at θ = θ^{*} are shown. (c,d) The same as panels (a) and (b) for the case 3. (e) Ratio between the maximum amplitudes of and at the surface (r = 0.985 R_{⊙}) with different combinations of the model parameters, Δr_{bmr} and Δϕ_{bmr}. The parameters used in cases 1–3 are denoted by red, green, and blue circles. The white dotted line represents the contour line of unity. 

In the text 
Fig. A.2. Radial component of the magnetic pressure and magnetic tension forces and shown in the same way as in Fig. A.1. (e) Ratio between the maximum amplitudes of and at the polarityinversion line (ϕ = 0). 

In the text 
Fig. B.1. Temporal evolution of the cyclic dynamo simulation with weak BabcockLeighton αeffect (a_{0} = 50 km s^{−1}). (a) Magnetic butterfly diagram at the surface r = 0.985 R_{⊙}. (b) Temporal evolution of the volumeintegrated kinetic and magnetic energies. The same definition of colors as in Fig. 10 is used. 

In the text 
Fig. B.2. Snapshot of the radial field B_{r} at the surface (top panels) and the radial velocity v_{r} near the surface (bottom panels) for the dynamo simulation with weak BabcockLeighton αeffect source term (a_{0} = 50 km s^{−1}). The same notation as Fig. 12 is used. 

In the text 
Fig. B.3. Temporal evolution of the longitudinallyaveraged magnetic fields and horizontal velocities for the dynamo simulation with deeper BabcockLeighton αeffect (Δr_{bmr} = 0.12 R_{⊙}). (a) Azimuthallyaveraged radial field ⟨B_{r}⟩ at the surface r = 0.985 R_{⊙}. (b) Azimuthallyaveraged toroidal field ⟨B_{ϕ}⟩ near the base of the convection zone r = 0.715 R_{⊙}. Black solid lines are the contours of the emerged BMRs at each time. (c) Torsional oscillation pattern δ⟨Ω⟩ at the surface. (d) Azimuthallyaveraged latitudinal velocity ⟨v_{θ}⟩ near the base of the convection zone. 

In the text 
Fig. B.4. Snapshot of the radial field B_{r} at the surface (top panels) and the radial velocity v_{r} near the surface (bottom panels) for the dynamo simulation with deeper BabcockLeighton αeffect (Δr_{bmr} = 0.12 R_{⊙}). The same notation as Fig. 12 is used. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.