© The Authors 2023

1. Introduction

The Sun exhibits an 11-yr cyclic magnetic activity that is sustained by the dynamo processes in the convection zone (e.g., Charbonneau 2020). The Babcock-Leighton flux-transport 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 so-called Babcock-Leighton α-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 east-and-west 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 east-west alignment (with zero tilt), on average, and that the tilts characterized by Joy’s law are generated by the north-south separation motions after emergence (Schunker et al. 2020).

Numerical investigations of the Babcock-Leighton flux-transport dynamo model have mainly been carried out in a two-dimensional (2D) kinematic mean-field framework (e.g., Dikpati & Charbonneau 1999; Nandy & Choudhuri 2001; Chatterjee et al. 2004; Hazra et al. 2014; Karak & Cameron 2016). In these models, the Babcock-Leighton α-effect is modeled as the axisymmetric poloidal source term which is localized near the surface. Although there are some non-kinematic studies where the dynamo-generated 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 Babcock-Leighton process in a three-dimensional (3D) full-spherical domain. Yeates & Muñoz-Jaramillo (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 Babcock-Leighton 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 so-called “double-ring” algorithm used in 2D mean-field models (Durney 1997; Nandy & Choudhuri 2001; Muñoz-Jaramillo et al. 2010). The same model has also been used to study the long-term cycle variability (Karak & Miesch 2017). However, all of these models are restricted to kinematic regime. Therefore, it still remains unclear how the Lorentz-forces of the BMRs affect their post-emergence evolution and the resulting dynamo solution in the non-kinematic 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 large-scale 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 flux-emergence 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 mean-field¹ models in which the large-scale mean-flows are largely controlled with parameterizations of small-scale 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 Babcock-Leighton dynamo processes of the Sun in a 3D non-kinematic regime, which takes advantage of both the mean-field approach for the solar large-scale mean flows and the 3D realization of the Babcock-Leighton process. Therefore, our model extends both the 2D non-kinematic mean-field 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 mean-field model was recently presented by Pipin (2022), in which the mean-field induction equation is solved together with the mean-field hydrodynamic equations. Although the model considers the non-axisymmetric magnetic field, such as newly-emerged BMRs, all the hydrodynamic variables were assumed to be axisymmetric and only the longitudinally-averaged Lorentz force was taken into account. In our model, by contrast, we consider the non-axisymmetric Lorentz force feedback on the non-axisymmetric flows, which we find crucial with regard to the post-emergence 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 dynamo-generated magnetic fields. We believe that our 3D non-kinematic 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 post-emergence 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, θ, ϕ):

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_1}{\partial t}& =-\mathrm{\nabla }·({\rho }_0\mathit{v}),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{v}}{\partial t}& =-\mathit{v}·\mathrm{\nabla }\mathit{v}-\frac{\mathrm{\nabla }{p}_1}{{\rho }_0}-\frac{{\rho }_1}{{\rho }_0}g{\mathit{e}}_{\mathrm{r}}+2\mathit{v}\times {\mathbf{\Omega }}_{\mathbf{0}}\hfill \\ \hfill & +\frac{1}{4\pi {\rho }_0}(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}+\frac{1}{{\rho }_0}\mathrm{\nabla }·\mathbf{\Pi },\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =\mathrm{\nabla }\times (\mathit{v}\times \mathit{B}+\mathcal{E}-\eta \mathrm{\nabla }\times \mathit{B}),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial s_1}{\partial t}& =\mathit{v}·\mathrm{\nabla }{s}_1+c_{\mathrm{p}}\delta \frac{v_r}{H_{\mathrm{p}}}+\frac{1}{{\rho }_0{T}_0}\mathrm{\nabla }·({\rho }_0{T}_0\kappa \mathrm{\nabla }{s}_1)\hfill \\ \hfill & +\frac{1}{{\rho }_0{T}_0}[(\mathbf{\Pi }·\mathrm{\nabla })·\mathit{v}+\frac{\eta }{4\pi }{|\mathrm{\nabla }\times \mathit{B}|}^2],\hfill \end{array}$$(4)

where g, ρ₀, p₀, and H_p denote the gravitational acceleration, density, pressure, and pressure scale height of the background state that is in an adiabatically-stratified 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, ρ₁, and p₁, are the perturbations with respect to the background that are assumed to be sufficiently small, namely, |p₁/p₀|≈|ρ₁/ρ₀|≪1, so that the equation of state is linearized

$$\begin{array}{c}\hfill p_1=p_0(\gamma \frac{{\rho }_1}{{\rho }_0}+\frac{s_1}{c_{\mathrm{v}}}),\end{array}$$(5)

where γ = 5/3 is the specific heat ratio and s₁ is entropy perturbation from the adiabatic background. The rotation rate of the radiative core Ω₀/2π = 431.3 nHz is used for a system rotation rate.

The tensor Π represents the turbulent Reynolds stress associated with small-scale (subgrid-scale) 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:

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{\mathit{ik}}={\rho }_0{\nu }_{\mathrm{vis}}(S_{\mathit{ik}}-\frac{2}{3}{\delta }_{\mathit{ik}}\mathrm{\nabla }·\mathit{v}+{\mathrm{\Lambda }}_{\mathit{ik}}{\mathrm{\Omega }}_0),\end{array}$$(6)

where S_ik and δ_ik denote the velocity deformation tensor and the Kronecker-delta 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¹² cm² s⁻¹ at the top boundary and is on the order of 10¹¹ cm² s⁻¹ 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 mixing-length model (Muñoz-Jaramillo et al. 2011), it is still about one order magnitude larger than the typical diffusivity value used in the kinematic flux-transport models in the advection-dominated regime (e.g., Dikpati & Charbonneau 1999; Guerrero 2007).

In order to break the Taylor-Proudman’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:

$$\begin{array}{cc}& \delta (r,\theta )={\mathcal{T}}_{-}(r;r_{\mathrm{sub}},d_{\mathrm{sub}}){\delta }_{\mathrm{sub}}(\theta ),\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\delta }_{\mathrm{sub}}(\theta )={\delta }_{\mathrm{pl}}+({\delta }_{\mathrm{eq}}-{\delta }_{\mathrm{pl}}){sin}^2\theta ,\hfill \end{array}$$(8)

$$\begin{array}{cc}& r_{\mathrm{sub}}(\theta )=r_{\mathrm{pl}}+(r_{\mathrm{eq}}-r_{\mathrm{pl}}){sin}^2\theta ,\hfill \end{array}$$(9)

where 𝒯 denotes a transition function defined by:

$$\begin{array}{c}\hfill {\mathcal{T}}_{\pm }(x;x_0,d)=\frac{1}{2}[1\pm tanh\left(\frac{x-x_0}{d}\right)].\end{array}$$(10)

We set the superadiabaticity at the poles and at the equator as δ_pl = −1.5 × 10⁻⁵ and δ_eq = −2 × 10⁻⁵, 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 non-local 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 low-entropy 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),

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{r\varphi }=+{\mathrm{\Lambda }}_0{\stackrel{\sim }{f}}_{\mathrm{l}}(r,\theta )cos(\theta +\lambda )[1+{\zeta }_r(r,\theta ,\varphi )],\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\theta \varphi }=-{\mathrm{\Lambda }}_0{\stackrel{\sim }{f}}_{\mathrm{l}}(r,\theta )sin(\theta +\lambda )[1+{\zeta }_{\theta }(r,\theta ,\varphi )].\hfill \end{array}$$(12)

The overall amplitude of the Λ-effect is given by Λ₀ = 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:

$$\begin{array}{cc}& {\stackrel{\sim }{f}}_{\mathrm{l}}(r,\theta )=\frac{f_{\mathrm{l}}(r,\theta )}{\max |f_{\mathrm{l}}(r,\theta )|},\hfill \end{array}$$(13)

$$\begin{array}{cc}& f_{\mathrm{l}}(r,\theta )={sin}^2\theta cos\theta tanh\left(\frac{r_{\max }-r}{d_{\mathrm{l}}}\right),\hfill \end{array}$$(14)

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:

$$\begin{array}{c}\hfill \zeta (r,\theta ,\varphi )={\displaystyle \sum _{i=1}^N{c}_{i}exp[-{\left(\frac{r-r_i}{\delta r}\right)}^2-{\left(\frac{\theta -{\theta }_i}{\delta \theta }\right)}^2-{\left(\frac{\varphi -{\varphi }_i}{\delta \varphi }\right)}^2],}\end{array}$$(15)

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 non-axisymmetric flows can be partially driven by these random fluctuations of the Λ-effect.

In the induction equation (Eq. (3)), we add an electro-motive-force ℰ to model the Babcock-Leighton α-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 fourth-order centered-difference method for space and a four-step Runge-Kutta 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 (9-wave method) for minimizing the numerical error resulting from the divergence of magnetic field (Dedner et al. 2002).

The numerical domain is a full-spherical 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 Yin-Yang grid (Kageyama & Sato 2004). For more details about the implementation of the Yin-Yang 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 stress-free 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 large-scale mean flows become quasi-stationary. Figure 1 shows profiles of the differential rotation, ⟨Ω⟩=Ω₀ + ⟨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 ρ0vm = ∇ × (Ψeϕ) where vm = ⟨vr⟩er + ⟨vθ⟩eθ. The meridional circulation is counter-clockwise (clockwise) in the northern (southern) hemisphere, i.e., the flow is poleward (equatorward) at the surface (base).

3. Post-emergence evolution of BMRs

In this section, we carry out a set of numerical simulations to study how the post-emergence 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 short-term evolution of a BMR and do not discuss the long-term 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 east-west aligned. The initial magnetic field is given as

$$\begin{array}{c}\hfill \mathit{B}\propto \mathrm{\nabla }\times (A_{\mathrm{ic}}{\mathit{e}}_{\theta }),\end{array}$$(16)

where the vector potential, A_ic, is given by:

$$\begin{array}{c}\hfill A_{\mathrm{ic}}(r,\theta ,\varphi )=\{1-tanh\left(\frac{l_{\mathrm{bmr}}(r,\varphi )-1}{0.5}\right)\}exp[-{\left(\frac{\theta -{\theta }^{\ast }}{\mathrm{\Delta }{\theta }_{\mathrm{bmr}}}\right)}^2],\end{array}$$(17)

with

$$\begin{array}{c}\hfill l_{\mathrm{bmr}}(r,\varphi )=\sqrt{{\left(\frac{r-r_{\max }}{\mathrm{\Delta }{r}_{\mathrm{bmr}}}\right)}^2+{\left(\frac{\varphi -{\varphi }^{\ast }}{\mathrm{\Delta }{\varphi }_{\mathrm{bmr}}}\right)}^2}.\end{array}$$(18)

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 vertically-elongated 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.

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 non-kinematic 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 Br at the surface r = 0.985 R⊙. Thick black solid curves show the contours at |Br| = 0.3 kG. The blue and red cross marks represent the locations of the flux-weighted 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, vr, 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 flux-weighted center locations of the leading and following polarity regions (θ_L, ϕ_L) and (θ_F, ϕ_F), respectively. The tilt angle, γ, is then defined as

$$\begin{array}{cc}\hfill \gamma & ={tan}^{-1}\left(\frac{sin{\theta }_{\mathrm{L}}-sin{\theta }_{\mathrm{F}}}{cos{\theta }_{\mathrm{L}}sin{\theta }_{\mathrm{L}}-cos{\theta }_{\mathrm{F}}sin{\theta }_{\mathrm{F}}}\right)\hfill \\ \hfill & \approx {tan}^{-1}\left[\frac{{\theta }_{\mathrm{L}}-{\theta }_{\mathrm{F}}}{({\varphi }_{\mathrm{L}}-{\varphi }_{\mathrm{F}})cos{\theta }^{\ast }}\right]·\hfill \end{array}$$(19)

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 anti-Joy’s law tilt is essentially 3D non-kinematic effect, and thus, it cannot be captured neither by the 2D non-kinematic 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 anti-Joy’s law tilt and the morphological asymmetry of the field strengths the BMRs. (a) Cross-section 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) Cross-section at the fixed latitude seen from the equator to the north pole. (c) Three-dimensional 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 anti-Joy’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 as-yet-uncaptured physical process. For example, Martin-Belda & Cameron (2016) proposed that the net positive tilt can be generated from the initial zero-tilt 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 well-known 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.

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

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

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 Lorentz-force-driven rising motion of the bottom apex of the flux tube drives the counter-rotating 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 |BL/BF|. 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 half-torus-shaped BMR with zero initial tilt, we find that its temporal evolution is extremely sensitive to its subsurface structure in the 3D non-kinematic regime. In particular, this study warns that the shallow BMRs model (which is conventionally used in many 2D non-kinematic models and 3D kinematic models) will lead to drastically different dynamo results when all these realistic 3D non-kinematic effects are included.

4. Cyclic dynamo with Babcock-Leighton α-effect

In the previous section, we show that the newly emerged BMRs with shallow subsurface root have a general tendency to produce the anti-Joy’s law tilts in the 3D non-kinematic 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. Babcock-Leighton α-effect source term

Now, we switch on the electro-motive-force ℰ in Eq. (3) that represents the source of the Babcock-Leighton α-effect, by which the surface poloidal field is produced as a result of the north-south 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 dynamo-generated 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ñoz-Jaramillo (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:

$$\begin{array}{c}\hfill {\overline{B}}_{\mathrm{tor}}(\theta ,\varphi )=\frac{1}{r_{\mathrm{b}}-r_{\mathrm{a}}}{\displaystyle {\int }_{r_{\mathrm{a}}}^{r_{\mathrm{b}}}{B}_{\varphi }(r,\theta ,\varphi )\mathrm{d}r,}\end{array}$$(20)

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 ${\overline{B}}_{\mathrm{tor}}(\theta ,\varphi )$ such that

$$\begin{array}{cc}\hfill {\overline{B}}_{\mathrm{tor}}^{\ast }(\theta ,\varphi )& ={\mathcal{T}}_{+}(\theta ;\pi /2-{\theta }_{\mathrm{em}},\mathrm{\Delta }{\theta }_{\mathrm{tran}})\hfill \\ \hfill & \times {\mathcal{T}}_{-}(\theta ;\pi /2+{\theta }_{\mathrm{em}},\mathrm{\Delta }{\theta }_{\mathrm{tran}}){\overline{B}}_{\mathrm{tor}}(\theta ,\varphi ),\hfill \end{array}$$(21)

where θ_em = 17.5° and Δθ_tran = 8.5°. We impose a necessary condition for the BMRs emergence to occur, that $|{\overline{B}}_{\mathrm{tor}}^{\ast }(\theta ,\varphi )|$ 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 ${\overline{B}}_{\mathrm{tor}}^{\ast }({\theta }^{\ast },{\varphi }^{\ast })$,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\mathcal{E}}_r\\ {\mathcal{E}}_{\theta }\\ {\mathcal{E}}_{\varphi }\end{array}\right)=a_0{\stackrel{\sim }{f}}_{\alpha }^{\ast }(r,\theta ,\varphi )\left(\begin{array}{c}0\\ -cos{\psi }^{\ast }\\ sin{\psi }^{\ast }\end{array}\right){\overline{B}}_{\mathrm{tor}}^{\ast }({\theta }^{\ast },{\varphi }^{\ast }),\end{array}$$(22)

where ${\stackrel{\sim }{f}}_{\alpha }^{\ast }$ represents the spatial distribution of BMRs,

$$\begin{array}{c}\hfill {\stackrel{\sim }{f}}_{\alpha }^{\ast }(r,\theta ,\varphi )=exp[-{\left(\frac{r-r_{\max }}{\mathrm{\Delta }{r}_{\mathrm{bmr}}}\right)}^2-{\left(\frac{\theta -{\theta }^{\ast }}{\mathrm{\Delta }{\theta }_{\mathrm{bmr}}}\right)}^2-{\left(\frac{\varphi -{\varphi }^{\ast }}{\mathrm{\Delta }{\varphi }_{\mathrm{bmr}}}\right)}^2].\end{array}$$(23)

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 non-kinematic 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 Babcock-Leighton α-effect is set to a₀ = 100 km s⁻¹. This value, in combination with the typical toroidal field strength near the base ${\overline{B}}_{\mathrm{tor}}^{\ast }\approx 5-20$ kG (Dikpati & Charbonneau 1999), leads to the total magnetic flux of BMRs of 10²² − 10²³ Mx, which is consistent with observations (Schrijver & Harvey 1994).

The north-south tilt of the BMRs (ψ^*) obeys Joy’s law such that

$$\begin{array}{c}\hfill {\psi }^{\ast }={35}^{°}cos{\theta }^{\ast }+{\psi }_{\mathrm{f}}^{\prime },\end{array}$$(24)

where ${\psi }_{\text{f}}^{\prime }$ 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 ${\psi }_{\text{f}}^{\prime }$ is roughly given by the following Gaussian distribution,

$$\begin{array}{c}\hfill P_{\mathrm{f}}({\psi }_{\mathrm{f}}^{\prime })=\frac{1}{{\sigma }_{\mathrm{f}}\sqrt{2\pi }}exp[-{\psi }_{\mathrm{f}}^{\prime 2}/(2{\sigma }_{\mathrm{f}}^2)],\end{array}$$(25)

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 self-consistently owing to the Lorentz-force feedback (Rempel 2006; Ichimura & Yokoyama 2017). Figure 7 shows examples of ℰ and the resulting tilted BMRs produced by our Babcock-Leighton α-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 Babcock-Leighton α-effect model. Radial field at the solar surface is shown where red (blue) points represent positive (negative) Br. The Solid black arrows denote the direction of the electro-motive-force, ℰ, 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 log-normal distribution function of the emergence events defined as:

$$\begin{array}{cc}& {\mathcal{C}}_{\mathrm{em}}({\mathrm{\Delta }}_{\mathrm{t}})={\displaystyle {\int }_{t=t_{\mathrm{s}}}^{t_{\mathrm{s}}+{\mathrm{\Delta }}_{\mathrm{t}}}{P}_{\mathrm{em}}({\mathrm{\Delta }}_{\mathrm{t}})\mathrm{d}t,}\hfill \end{array}$$(26)

$$\begin{array}{cc}& P_{\mathrm{em}}({\mathrm{\Delta }}_{\mathrm{t}})=\frac{1}{{\sigma }_{\mathrm{t}}{\mathrm{\Delta }}_{\mathrm{t}}\sqrt{2\pi }}exp[-\frac{{(ln{\mathrm{\Delta }}_{\mathrm{t}}-{\mu }_{\mathrm{t}})}^2}{2{\sigma }_{\mathrm{t}}^2}],\hfill \end{array}$$(27)

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.

$$\begin{array}{cc}& {\sigma }_{\mathrm{t}}^2=\frac{2}{3}ln\left(\frac{{\tau }_{\mathrm{s}}}{{\tau }_{\mathrm{p}}}\right),{\mu }_{\mathrm{t}}={\sigma }_{\mathrm{t}}^2+ln{\tau }_{\mathrm{p}},\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill & {\tau }_{\mathrm{p}}={\tau }_{\mathrm{p},0}+\mathrm{\Delta }{\tau }_{\mathrm{p}}{e}^{-{({\overline{B}}_{\mathrm{em}}^{\ast }/B_{\tau })}^2},\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill & {\tau }_{\mathrm{s}}={\tau }_{\mathrm{s},0}+\mathrm{\Delta }{\tau }_{\mathrm{s}}{e}^{-{({\overline{B}}_{\mathrm{em}}^{\ast }/B_{\tau })}^2}.\hfill \end{array}$$(30)

Here, we set τ_p, 0 = 0.8 days, τ_s, 0 = 1.9 days, Δτ_p = 0.75 days, and Δτ_s = 3.0 days. The quantity ${\overline{B}}_{\mathrm{em}}^{\ast }$ represents the horizontally-averaged ${\overline{B}}_{\mathrm{tor}}^{\ast }$ 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 (${\overline{B}}_{\mathrm{em}}^{\ast }>B_{\tau }$) and less frequent during the activity minima (${\overline{B}}_{\mathrm{em}}^{\ast }<B_{\tau }$).

Fig. 8. Cumulative log-normal 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 Babcock-Leighton α-effect model is strongly spatially-localized and temporally intermittent. This is clearly different from the conventional 2D models with spatially-distributed and temporally-continuous 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 α₀ value within the 2D mean-field framework as

$$\begin{array}{cc}\hfill {\alpha }_0& \approx sin\psi \left(\frac{\mathrm{\Delta }{\varphi }_{\mathrm{bmr}}}{2\pi }\right)\left(\frac{\mathrm{\Delta }{t}_{\mathrm{CFL}}}{{\mathrm{\Delta }}_{\mathrm{t}}}\right)a_0\hfill \\ \hfill & \approx 1.2\mathrm{m}{\mathrm{s}}^{-1},\hfill \end{array}$$(31)

where we use the typical values ψ = 17.5°, Δ_t = 5 days, and Δt_CFL = 17 min. This α₀ value is consistent with the previous 2D mean-field models.

4.2. Initial condition of magnetic fields

We add an axisymmetric dipolar field into the fully-developed hydrodynamic calculation shown in Fig. 1. The simulation is evolved until initial transients disappear and the dynamo cycles with quasi-steady 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 time-latitude plots of the longitudinally-averaged radial field ⟨B_r⟩ at the surface and the toroidal field ⟨B_ϕ⟩ near the base of the convection zone, represented in terms of the well-known 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 build-up of the polar field by poleward advection of the magnetic fluxes associated with the trailing sunspots. These are owing to the single-cell meridional circulation achieved in our model, which has an amplitude of about 15 m s⁻¹ at the surface and 2 m s⁻¹ 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 so-called 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 longitudinally-averaged magnetic fields and horizontal velocities. (a) Azimuthal mean of the radial field ⟨Br⟩ 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 δ⟨s1⟩=⟨s1⟩−⟨s1⟩t at the surface.

In our non-kinematic model, the dynamo-generated fields have strong impacts on flows via the Lorentz force feedback. Figure 9c shows the time-latitude 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 time-latitude 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 time-latitude plot of the entropy perturbation δ⟨s₁⟩=⟨s₁⟩−⟨s₁⟩_t at the surface with typical variation amplitude of about 250 erg g⁻¹ K⁻¹, 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 low-latitude 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 low-latitude torsional oscillation branches due to the thermal forcing with the opposite sign.

Figure 10 shows the volume-integrated kinetic and magnetic energies of various components. Their definitions are as follows:

$$\begin{array}{cc}& {\mathrm{KE}}_{\mathrm{DR}}={\displaystyle {\int }_V\frac{{\rho }_0}{2}{⟨v_{\varphi }⟩}^2\mathrm{d}V,}\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\mathrm{KE}}_{\mathrm{MC}}={\displaystyle {\int }_V\frac{{\rho }_0}{2}({⟨v_r⟩}^2+{⟨v_{\theta }⟩}^2)\mathrm{d}V,}\hfill \end{array}$$(33)

$$\begin{array}{cc}& {\mathrm{KE}}_{m\ne 0}={\displaystyle {\int }_V\frac{{\rho }_0}{2}{(\mathit{v}-⟨\mathit{v}⟩)}^2\mathrm{d}V,}\hfill \end{array}$$(34)

$$\begin{array}{cc}& {\mathrm{ME}}_{\mathrm{tor}}={\displaystyle {\int }_V\frac{1}{8\pi }{⟨B_{\varphi }⟩}^2\mathrm{d}V,}\hfill \end{array}$$(35)

$$\begin{array}{cc}& {\mathrm{ME}}_{\mathrm{pol}}={\displaystyle {\int }_V\frac{1}{8\pi }({⟨B_r⟩}^2+{⟨B_{\theta }⟩}^2)\mathrm{d}V,}\hfill \end{array}$$(36)

$$\begin{array}{cc}& {\mathrm{ME}}_{m\ne 0}={\displaystyle {\int }_V\frac{1}{8\pi }{(\mathit{B}-⟨\mathit{B}⟩)}^2\mathrm{d}V,}\hfill \end{array}$$(37)

Fig. 10. Temporal evolution of the volume-integrated kinetic and magnetic energies. The red, green, and black lines denote the kinetic energies of the differential rotation KEDR, the meridional circulation KEMC, and the non-axisymmetric flows KEm ≠ 0. The blue, purple, and orange lines denote the magnetic energies of the mean toroidal field KEtor, the mean poloidal field MEpol, and the non-axisymmetric fields MEm ≠ 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 non-axisymmetric magnetic energy, ME_m ≠ 0, is greater than the mean poloidal field energy ME_pol because the BMRs are strongly non-axisymmetric. The non-axisymmetric fields drive strong non-axisymmetric flows, whose kinetic energy, KE_m ≠ 0, is also greater than ME_pol. This implies that the non-axisymmetric 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 Babcock-Leighton 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 non-axisymmetric. 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 mean-field 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 Br 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. Non-axisymmetric flows

In our model, non-axisymmetric flows are driven largely by the non-axisymmetric 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⁻¹. 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 newly-emerged BMR that initially consists of two round-shaped 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 horizontally-elongated BMRs (case 1). We find that this elongated feature of the BMRs, which is not observed on the Sun (van Driel-Gesztelyi & Green 2015), can be greatly suppressed if the Babcock-Leighton α-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 Br at the surface (top panels) and the radial velocity vr 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 zoom-in of the panels (a) and (b), focusing on the single BMR, denoted by red thick solid lines.

Other interesting non-axisymmetric flow features are low-frequency 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 non-kinematic 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 non-axisymmetric random fluctuations in the Λ-effect and by the non-axisymmetric Lorentz-force; 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 Babcock-Leighton flux-transport dynamo model of the Sun. In our model, we do not solve the small-scale convection and focus on the large-scale flows and magnetic structures in a full spherical shell. The solar-like large-scale mean flows are driven by proper parameterization of the Λ-effect. The model operates in a 3D non-kinematic regime, and therefore, is more realistic than the 2D non-kinematic models (e.g., Rempel 2006; Ichimura & Yokoyama 2017) and 3D kinematic models (e.g., Yeates & Muñoz-Jaramillo 2013; Miesch & Dikpati 2014).

To better illustrate the major differences from the conventional 2D non-kinematic 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 post-emergence evolution of the BMR becomes significantly changed from those from the conventional models: Even if the initial BMR is perfectly east-west 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 Babcock-Leighton α-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 Babcock-Leighton 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 newly-emerged 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: Λ₀, which determines the amplitudes of the differential rotation and meridional circulation, and a₀, which determines the field strengths of BMRs at the surface. A detailed parameter study is required in the future.

Due to the 3D non-kinematic nature of our model, we find the substantial non-axisymmetric flows that are driven by the Lorentz-force 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 (Martin-Belda & 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 short-term post-emergence evolution of the BMRs and the long-term 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 (Martin-Belda & 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 Babcock-Leighton solar dynamo (Cameron & Schüssler 2012). Furthermore, it is often argued that the low-latitude branches of the torsional oscillation are attributed to the thermally-induced 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 post-emergence 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 deep-seated 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.

Movie

Movie 1 associated with Fig. 11 (animation_cyclic_dynamo) Access here

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 Max-Planck Research School for Solar System Science at the University of Göttingen (IMPRS). Y.B. also acknowledges a support from a long-term scholarship program for degree-seeking 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 Max-Planck supercomputer RZG in Garching.

References

Appendix A: Lorentz force of a model BMR

In Sect. 3, we show that the post-emergence 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:

$$\begin{array}{c}\hfill \frac{1}{4\pi }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}=\underset{{\mathit{F}}^{\mathrm{MP}}}{\underbrace{-\mathrm{\nabla }\left(\frac{{\mathit{B}}^2}{8\pi }\right)}}+\underset{{\mathit{F}}^{\mathrm{MT}}}{\underbrace{\frac{1}{4\pi }(\mathit{B}·\mathrm{\nabla })\mathit{B}}}.\end{array}$$(A.1)

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 $F_{\varphi }^{\mathrm{MT}}$, 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 $F_{\varphi }^{\mathrm{MT}}$ and $F_{\varphi }^{\mathrm{MP}}$ 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 (anti-Joy’s law) tilt owing to the strong tension-force-driven longitudinal converging flows and the subsequent Coriolis force acting on them.

Fig. A.1. Longitudinal component of the magnetic pressure and magnetic tension forces, $F_{\varphi }^{\mathrm{MP}}$ and $F_{\varphi }^{\mathrm{MT}}$. (a,b) $F_{\varphi }^{\mathrm{MP}}$ and $F_{\varphi }^{\mathrm{MT}}$ for the initial BMR for the case 1, respectively. The cross-sections at θ = θ* are shown. (c,d) The same as panels (a) and (b) for the case 3. (e) Ratio between the maximum amplitudes of $F_{\varphi }^{\mathrm{MP}}$ and $F_{\varphi }^{\mathrm{MT}}$ at the surface (r = 0.985 R⊙) with different combinations of the model parameters, Δrbmr 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 $F_r^{\mathrm{MP}}$ 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 $F_r^{\mathrm{MT}}$ 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 $F_r^{\mathrm{MP}}$ and $F_r^{\mathrm{MT}}$ shown in the same way as in Fig. A.1. (e) Ratio between the maximum amplitudes of $F_r^{\mathrm{MP}}$ and $F_r^{\mathrm{MT}}$ at the polarity-inversion 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₀ = 50 km s⁻¹

We carried out the simulation with weaker Babcock-Leighton α-effect source, in which a₀ is decreased to 50 km s⁻¹ 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 longitudinally-averaged 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 volume-integrated 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 Babcock-Leighton α-effect (a0 = 50 km s−1). (a) Magnetic butterfly diagram at the surface r = 0.985 R⊙. (b) Temporal evolution of the volume-integrated 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 non-axisymmetric 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 Br at the surface (top panels) and the radial velocity vr near the surface (bottom panels) for the dynamo simulation with weak Babcock-Leighton α-effect source term (a0 = 50 km s−1). The same notation as Fig. 12 is used.

We find that in the simulation with weak Babcock-Leighton α-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 newly-emerged 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 time-latitude 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 high-latitude branches, as shown in Fig. B.3c. We confirm that the longitudinally-averaged 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 longitudinally-averaged magnetic fields and horizontal velocities for the dynamo simulation with deeper Babcock-Leighton α-effect (Δrbmr = 0.12 R⊙). (a) Azimuthally-averaged radial field ⟨Br⟩ at the surface r = 0.985 R⊙. (b) Azimuthally-averaged 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) Azimuthally-averaged 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 Babcock-Leighton α-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 Babcock-Leighton α-effect has an effective α_ϕϕ which is positive (negative) in the northern (southern) hemisphere. According to the Parker-Yoshimura sign rule (Parker 1955; Yoshimura 1975), the propagation direction of the αΩ dynamo wave s_αΩ is given by:

$$\begin{array}{c}\hfill {\mathit{s}}_{\alpha \mathrm{\Omega }}\propto {\alpha }_{\varphi \varphi }\mathrm{\nabla }⟨\mathrm{\Omega }⟩\times {\mathit{e}}_{\varphi }.\end{array}$$(B.1)

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 small-scale helical turbulence (e.g., Parker 1955; Moffatt 1978).

Figure B.4 shows the temporal snapshot of the radial magnetic field B_r and the non-axisymmetric flows at the surface from the simulation with deep Babcock-Leighton α-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 Br at the surface (top panels) and the radial velocity vr near the surface (bottom panels) for the dynamo simulation with deeper Babcock-Leighton α-effect (Δrbmr = 0.12 R⊙). The same notation as Fig. 12 is used.

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 ρ0vm = ∇ × (Ψeϕ) where vm = ⟨vr⟩er + ⟨vθ⟩eθ. The meridional circulation is counter-clockwise (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 Br at the surface r = 0.985 R⊙. Thick black solid curves show the contours at |Br| = 0.3 kG. The blue and red cross marks represent the locations of the flux-weighted 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, vr, 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 anti-Joy’s law tilt and the morphological asymmetry of the field strengths the BMRs. (a) Cross-section 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) Cross-section at the fixed latitude seen from the equator to the north pole. (c) Three-dimensional 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 |BL/BF|. 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 Babcock-Leighton α-effect model. Radial field at the solar surface is shown where red (blue) points represent positive (negative) Br. The Solid black arrows denote the direction of the electro-motive-force, ℰ, 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 log-normal 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 longitudinally-averaged magnetic fields and horizontal velocities. (a) Azimuthal mean of the radial field ⟨Br⟩ 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 δ⟨s1⟩=⟨s1⟩−⟨s1⟩t at the surface.

In the text

Fig. 10. Temporal evolution of the volume-integrated kinetic and magnetic energies. The red, green, and black lines denote the kinetic energies of the differential rotation KEDR, the meridional circulation KEMC, and the non-axisymmetric flows KEm ≠ 0. The blue, purple, and orange lines denote the magnetic energies of the mean toroidal field KEtor, the mean poloidal field MEpol, and the non-axisymmetric fields MEm ≠ 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 Br 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 Br at the surface (top panels) and the radial velocity vr 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 zoom-in 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, $F_{\varphi }^{\mathrm{MP}}$ and $F_{\varphi }^{\mathrm{MT}}$. (a,b) $F_{\varphi }^{\mathrm{MP}}$ and $F_{\varphi }^{\mathrm{MT}}$ for the initial BMR for the case 1, respectively. The cross-sections at θ = θ* are shown. (c,d) The same as panels (a) and (b) for the case 3. (e) Ratio between the maximum amplitudes of $F_{\varphi }^{\mathrm{MP}}$ and $F_{\varphi }^{\mathrm{MT}}$ at the surface (r = 0.985 R⊙) with different combinations of the model parameters, Δrbmr 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 $F_r^{\mathrm{MP}}$ and $F_r^{\mathrm{MT}}$ shown in the same way as in Fig. A.1. (e) Ratio between the maximum amplitudes of $F_r^{\mathrm{MP}}$ and $F_r^{\mathrm{MT}}$ at the polarity-inversion line (ϕ = 0).
In the text

	Fig. B.1. Temporal evolution of the cyclic dynamo simulation with weak Babcock-Leighton α-effect (a0 = 50 km s−1). (a) Magnetic butterfly diagram at the surface r = 0.985 R⊙. (b) Temporal evolution of the volume-integrated 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 Br at the surface (top panels) and the radial velocity vr near the surface (bottom panels) for the dynamo simulation with weak Babcock-Leighton α-effect source term (a0 = 50 km s−1). The same notation as Fig. 12 is used.
In the text

Fig. B.3. Temporal evolution of the longitudinally-averaged magnetic fields and horizontal velocities for the dynamo simulation with deeper Babcock-Leighton α-effect (Δrbmr = 0.12 R⊙). (a) Azimuthally-averaged radial field ⟨Br⟩ at the surface r = 0.985 R⊙. (b) Azimuthally-averaged 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) Azimuthally-averaged latitudinal velocity ⟨vθ⟩ near the base of the convection zone.

In the text

	Fig. B.4. Snapshot of the radial field Br at the surface (top panels) and the radial velocity vr near the surface (bottom panels) for the dynamo simulation with deeper Babcock-Leighton α-effect (Δrbmr = 0.12 R⊙). The same notation as Fig. 12 is used.
In the text