Toroidal Magnetic Flux Budget in Mean-field Dynamo Model of Solar Cycles 23 and 24

We study the toroidal magnetic flux budget of the axisymmetric part of a data-driven 3D mean-field dynamo model of Solar Cycles 23 and 24. The model simulates the global solar dynamo that includes the effects of the formation and evolution of bipolar magnetic regions (BMRs) emerging on the solar surface. By applying Stokes’s theorem to the dynamo induction equation, we show that the hemispheric magnitude of the net axisymmetric toroidal magnetic field generation rate in the bulk of the convection zone can only partially be estimated from the surface parameters of the differential rotation and the axisymmetric radial magnetic field. The contribution of the radial integral along the equator, which is mostly due to the rotational radial shear at the bottom of the convection zone, has the same magnitude and is nearly in phase with the effect of the surface latitudinal differential rotation. Also, the toroidal field generation rate estimate strongly depends on the latitudinal profile of the surface radial magnetic field near the poles. This profile in our dynamo models significantly deviates from the polar magnetic field distribution observed during the minima of Solar Cycles 22, 23, and 24. The cause of this discrepancy requires further observational and theoretical studies. Comparing the 2D axisymmetric and the 3D nonaxisymmetric dynamo models, we find an increase in the toroidal field generation rate in the 3D model due to the surface effects of BMRs, resulting in an increase in the axisymmetric poloidal magnetic field magnitude.


INTRODUCTION
The total unsigned radial magnetic flux observed in the photosphere during the Sun's 11-year activity cycle peaks at about 10 24 Mx (Schrijver & Harvey 1984;Schrijver & Harvey 1994).Some portion of this magnetic flux is likely produced by the large-scale dynamo in the convection zone.But how and where this flux is produced is under debate.A small fraction of it comes from active regions.Depending on the complexity of active regions, the total unsigned radial magnetic field flux associated with the bipolar sunspot groups (hereafter, bipolar magnetic regions, BMRs) peaks at about 4-10 × 10 22 Mx (Nagovitsyn et al. 2016;Abramenko et al. 2023).A part of this flux must participate in the surface flux transport and turbulent processes, e.g., via helical convection motion, contributing to the global dynamo.Following the basic ideas of Parker (1955), the solar dynamo cyclically transforms the large-scale poloidal magnetic field into the toroidal field and back by means of the differential rotation and the turbulent α effect caused by the cyclonic convective motions.Cameron & Schüssler (2015) (hereafter CS15) showed that the hemispheric magnitude of the net axisymmetric toroidal flux production in the bulk of the convection zone is largely determined by the surface parameters of the differential rotation and the radial magnetic field.Using solar observations, they found that the toroidal flux generation rate by the latitudinal rotational shear is about 10 23 Mx/yr or 10 15 Mx/s.Therefore, this effect is capable of producing the toroidal flux of 10 24 Mx during the solar cycles.It is interesting to compare this number with the typical unsigned radial magnetic field flux emergence rate of BMRs, which is about 10 18 Mx/s (Golovko 1998).Following the arguments of CS15, it is tempting to conclude that the poloidal magnetic flux produced by the magnetic flux of bipolar magnetic regions (BMRs) is sufficient for maintaining the solar dynamo.Indeed, this idea has been implemented in the Babcock-Leighton flux transport dynamo models (Babcock 1961;Leighton 1969;Hazra et al. 2023;Cameron & Schüssler 2023).It is worth noting that the evidences found by CS15 reflect effects of the global axisymmetric dynamo inside the Sun.Our aim is to analyze these effects using the 3D mean-field dynamo model developed by Pipin (2022); Pipin et al. (2023) (hereafter,P22;PKT23).
Despite the simplicity and visual appeal of the magnetic 'butterfly' diagram reproduced by the Babcock-Leighton flux transport dynamo models, this dynamo scenario is not supported by the global-Sun 3D MHD simulations (e.g., (Guerrero et al. 2016;Käpylä et al. 2016;Warnecke et al. 2018Warnecke et al. , 2021)).The simulations favor the dynamo scenario suggested by Parker (1955), in which a substantial portion of the poloidal magnetic flux is generated by helical turbulence deep in the convection zone.The Parker's dynamo scenario suggests that the solar magnetic activity represents itself in the form of the dynamo waves propagating in the convection zone along the isorotation surfaces.
The helioseismic observations support Parker's idea as well, showing the migrating zonal flows ('torsional oscillations'), which reveal a dynamical pattern corresponding to the dynamo waves predicted by Parker's theory (Kosovichev & Pipin 2019).The dynamical mean-field model developed by Pipin (2018) and Pipin & Kosovichev (2019, 2020) reproduced the observed properties of the torsional oscillations and its' extended' 20-year mode.The 'extended' cycle is observed in the evolution of the ephemeral active regions, the coronal green-line emission, the surface torsional oscillation (Wilson et al. 1988), and the evolution of the low-degree harmonics of the global magnetic fields (Obridko et al. 2023).The 'extended' cycle can be a result of the extended mode of the axisymmetric dynamo waves, and the modulation of the heat and angular momentum transport in the solar convection zone by the dynamo-generated magnetic fields (Stenflo 1994;Pipin & Kosovichev 2019).The model also showed the extended 22-year cycle in the time-latitude variations of the meridional circulation, which were found by helioseismic observations (Komm et al. 2018;Getling et al. 2021).
Further developments of our mean-field model included modeling the formation and emergence of BMRs driven by magnetic buoyancy instability (P22; PKT23).Comparing results of the 3D dynamo model runs with the 2D axisymmetric dynamo model, P22 found that the BMRs' activity may increase of the total toroidal magnetic field flux in the convection zone by about 10 percent.Also, the results of P22 and PKT23 showed that the surface unsigned radial magnetic field flux reaches a magnitude of 5 − 7 × 10 23 Mx in the magnetic cycle maximum.Most of this flux originates from the emergence and evolution of the BMRs.While the role of the surface BMR activity in dynamo is subdominant, the 3D model combines the essential elements of the Parker and Babcock-Leighton dynamo scenarios.
In this paper, we investigate the toroidal magnetic flux budget for the axisymmetric part of the dynamo model using the results of the 3D mean-field data-driven dynamo model of Solar Cycles 23 and 24 (PKT23).We compare the modeling results with the corresponding synoptic observations of magnetic fields from the Kitt Peak Observatory (KPO) and the Solar Dynamics Observatory (SDO/HMI).Section 2 briefly describes the model.Section 3 presents the results for the toroidal magnetic flux budgets and comparison with the synoptic observations.Finally, a summary of the main results is presented in Section 4. Some additional details of the BMR model, the basic dynamo model parameters, and the magnetic field evolution are described in Appendixes A and B.

DYNAMO MODEL
We employ the dynamo model recently developed by P22.In the model, we decompose the mean magnetic field, ⟨B⟩, (associated with the statistical ensemble average) into the sum of the axisymmetric, B, and the nonaxisymmetric parts, B, (1) These parts are further decomposed into a sum of the poloidal and toroidal components, as follows (Moss & Brandenburg 1992): where φ is the azimuthal unit vector, r is the radius vector, and θ is the polar angle.The gauge transformation for potentials T and S involves a sum with arbitrary r-dependent functions (Krause & Rädler 1980).The gauge dependence is removed if the integrals of T and S over the solid angle are zero.We seek the solution for T and S using the spherical harmonic decomposition.In this case, by definition, the m = 0 modes in T and S are excluded, and the representation of Eq.( 3) is gauge invariant (also, see Berger & Hornig 2018).
The evolution of the large-scale magnetic field is governed by the mean-field induction equation, where E = ⟨u × b⟩ is the mean electromotive force; u and b are randomly fluctuating velocity and magnetic field components.Their evolution equations are solved separately, either analytically or numerically using the so-called test-field method (Schrinner 2011).Typically, such a solution includes nonlinear effects that stem from effects of the global rotation and fluctuations of ⟨B⟩ caused by background turbulent motions and small-scale dynamo (Käpylä et al. 2022;Brandenburg et al. 2023).We employ the phenomenological term E (BMR) to describe that emergence of bipolar magnetic regions on the solar surface.
In our model, we assume that the large-scale flow is axisymmetric, ⟨U⟩ ≡ U.The same assumption is applied to the mean entropy and the other thermodynamic parameters.However, the magnetic field is three-dimensional.To get the dynamo equations for the axisymmetric toroidal magnetic field evolution, we take the scalar product of Eq.( 4) with the unit vector φ.We do the same for the uncurled version of Eq.( 4) to get the equation for the vector potential, A. To obtain the evolution equations for the nonaxisymmetric magnetic field, we apply the curl and double curl operations to Eq.( 4).Then, we take the scalar product of these equations with vector r (see details in Krause &Rädler 1980, andMoss &Brandenburg 1992).
We define E using the mean-field MHD magnetohydrodynamics framework of Krause & Rädler (1980); it reads where α ij describes the turbulent generation by kinetic and magnetic helicity, γ ij describes the turbulent pumping, and η ijk is the eddy magnetic diffusivity tensor.The α effect tensor includes effect of the magnetic helicity conservation (Kleeorin & Ruzmaikin 1982;Kleeorin & Rogachevskii 1999), Here, the pseudoscalar ⟨a • b⟩ (where b = ∇ × a) is the small-scale magnetic helicity density, C α is the dynamo parameter characterizing the magnitude of the hydrodynamic α effect, and α K ij and α M ij are the anisotropic versions of the kinetic and magnetic α effects (Pipin 2008;Pipin & Kosovichev 2019;Brandenburg et al. 2023).The radial profiles of α H ij and α M ij depend on the mean density stratification and the spatial profiles of the convective velocity u c , and on the Coriolis number, Co = 2Ω 0 τ c , where Ω 0 is the global angular velocity of the star and τ c is the convective turnover time.The magnetic quenching function ψ α (β) depends on the parameter β = |⟨B⟩|/ 4πρu 2 c (Pipin & Kosovichev 2019).The evolution of magnetic helicity is governed by the equation of the integral balance of the helicity density, where ⟨χ⟩ (tot) = ⟨a • b⟩+⟨A⟩•⟨B⟩ is the total magnetic helicity, η 0 is the constant microscopic diffusivity, ⟨J⟩ = ∇×⟨B⟩ is the mean current density, and parameter R m = u c ℓ c /η 0 = 10 6 is the magnetic Reynolds number, where u c and ℓ c are the RMS convective velocity and the mixing length of the turbulent convection, respectively.We choose η 0 to be small enough so that its amplitude does not affect the evolution of the small-scale helicity density.Following the results of Mitra et al. (2010), we use the mixing-length estimation of the magnetic diffusivity η χ = 1 10 η T , where, η T = u c ℓ c /3. Similarly to α ij , the structure and radial profiles of γ ij and η ijk are determined by the effects of the global rotation and large-scale magnetic field on turbulent flows (Pipin 2008;Pipin & Kosovichev 2019 and P22). Figure 6 in Appendix B shows the convection zone profiles for α K ij , η ijk and the effective drift velocity, which is defined as a sum of the turbulent pumping and meridional circulation (Warnecke et al. 2018).
The calculation of the small-scale helicity density evolution from Eq.( 8) depends on the gauge of the large-scale magnetic field vector potential, ⟨A⟩.The axisymmetric part of the vector-potential satisfies the Coulomb and Poincare gauges, by default.The nonaxisymmetric part of ⟨A⟩ depends on gauges of scalars S and T .In our spherical harmonic decomposition, we exclude the m = 0 modes of T and S. In this case, the representation of Eq.( 3) is gauge invariant (Krause & Rädler 1980;Berger & Hornig 2018).
The electromotive force induced by the bipolar magnetic regions is determined by Parameter α β determines the magnetic field generation effect, which is related to the BMRs tilt (P22).The second term of this equation models the emergence of the nonaxisymmetric magnetic field in the form of bipolar magnetic regions.Parameter V β is related to the buoyancy velocity of BMRs in the convection zone.The reader can find further details about E (BMR) in Appendix and P22; PKT23.
The reference profiles of the mean thermodynamic parameters, such as entropy, density, temperature, and the parameters of the background turbulent convection, such as the convective turnover time, τ c , and the mixing length, ℓ c , are determined from the solar interior model MESA (Paxton et al. 2013).For this computation, we use the default choice of the MESA input parameters, which includes the information on the solar mass (1M ⊙ ), metalicity (z=0.02), and age (4.6 Gyr).The convective RMS velocity is determined from the mixing-length approximation, where ℓ c = α MLT H p is the mixing length, α MLT = 1.9 is the mixing length parameter, and H p is the pressure scale height.Equation ( 10) determines the reference profiles for the eddy heat conductivity, χ T , eddy viscosity, ν T , and eddy diffusivity, η T , as follows, where Pr T = 3/4 is the turbulent Prandtl number.We put the magnetic turbulent Prandtl number Pm T = 10.This choice allows us to match the solar cycle period in our dynamo model.The model includes the effects of the dynamo-generated magnetic field on the evolution of global axisymmetric flows and the heat transport in the solar convection zone.We take into account the effect of the nonaxisymmetric magnetic field using the longitudinal averaging of the Lorentz force and the Maxwell stresses.The reader can find further details in our papers (Pipin & Kosovichev 2019;Pipin 2022;Pipin et al. 2023).The large-scale flow and turbulence parameters profiles of our 3D model as well as the modeled surface evolution of the radial magnetic field are presented in Appendix B.

The boundary conditions and parameters
Our models include the convection zone and the convective overshoot region.The bottom of the integration domain is fixed at the bottom of the overshoot zone, r i = 0.67R, where R ≡ R ⊙ is the radius of the Sun,.The convection zone extends from r b = 0.728R to r t = 0.99R.The top layers of the solar convection zone are not included because of the strong density variations in the very near surface layer, which can be problematic because of the numerical resolution issues.At the bottom boundary, r i = 0.67R, we set the magnetic field induction vector to zero.At the top boundary, we use the condition, which allows penetration of the toroidal magnetic field to the surface (Pipin & Kosovichev 2011).Parameter δ controls the magnitude of the axisymmetric toroidal field at the top boundary.For the set of parameters: δ = 0.999 and B esq = 5G, the surface toroidal field magnitude is around 1.5 G.The poloidal magnetic field is potential outside the dynamo domain.For the nonaxisymmetric magnetic field, we put potentials T and S to zero at the bottom boundary, Note-Model T0 is the axisymmetric baseline model without BMRs.Model S2 is the data-driven model with the BMR emergence.The forth column shows the magnitude of the unsigned toroidal flux in the convection zone and its generation rate; the fourth column shows the magnitude of the unsigned surface axisymmetric radial magnetic field flux; the last column shows the duration of the activity cycles (half dynamo periods of the magnetic cycles).
r i .At the top, we use the standard vacuum boundary conditions, i.e, T = 0, and S is defined by the potential magnetic field solution outside of the dynamo domain.The dynamo model equations are solved using the spherical harmonics with angular degree ℓ ≤ 72.For the numerical solution, we employ the fortran version of the shtns library of Schaeffer (2013).To integrate over the radius, we employ the finite differences using the uniform mesh of 80 points over r b − r t interval and 30 mesh points for the region below the convection zone, r i − r b .The model employs the second-order Crank-Nicholson scheme with alternate directions for integration in time.
The full dynamo cycle of about 20 years is reproduced if the turbulent magnetic Prandtl number Pm T = 10, and the α-effect parameter C α = 0.045.The critical value of this parameter C cr α ≈ 0.04.The magnitude of the alpha effect in our model is about 0.5 m/s.These parameters are adopted in the reference axisymmetric model T0.In the range C α = 0.04 − 0.045, the dynamo solution shows the antisymmetric magnetic field about the equator.The magnetic flux loss due to the mean-field magnetic buoyancy results in a weakly nonlinear dynamo regime with the maximum toroidal magnetic field strength of about 2.5 kG and parameter β = |⟨B⟩| / 4πρu 2 c ≤ 0.2 in the solar convection zone.We evolve the axisymmetric model to a quasi-stationary stage and choose the profile of the magnetic field at the minimum of a magnetic cycle with the sign of the polar magnetic field, which corresponds to the minimum state of Solar Cycle 22 on Jan 1, 1996.This magnetic field profile is considered as the initial state of the data-driven non-axisymmetric model S2.In this model, a part of the axisymmetric poloidal magnetic field flux is generated due to the evolution of the surface BMRs.Therefore, in this model, we reduced the baseline magnitude of the mean-field α effect to the threshold value of the axisymmetric dynamo model (C α = C cr α ).To reproduce the magnitude and duration of Solar Cycles 23 and 24, we further decreased C α = 0.036 after five years from the start of Cycle 23, and increased it back to C α = C cr α at the end of the cycle in 2007, see the further details of the data-driven model S2 in PKT23.To evaluate the effects of BMRs in the magnetic flux budget, we compare the non-axisymmetric data-driven model S2 with the axisymmetric reference model T0 calculated without BMRs.
Some output parameters of our models are given in Table 1, where we use the same model notations as in PKT23.The model results were extensively discussed in our previous papers (Pipin & Kosovichev 2019;Pipin 2022;Pipin et al. 2023).For convenience, we included some of those results in the Appendices.

MAGNETIC FLUX BUDGET
To estimate variations of the axisymmetric toroidal magnetic flux in the dynamo region, we follow the approach of CS15 and apply Stokes' theorem to the induction equation.The time derivative of the axisymmetric toroidal magnetic field flux in the Northern hemisphere of the Sun is calculated as follows, where Φ N tor = Σ B ϕ dS, Σ is the area of the meridional cut through the solar convection zone in the hemisphere, δΣ stands for the contour line confining the cut, and dl is the line element of δΣ; U, B, and E, E (BM R) are the large-scale axisymmetric flow velocity, magnetic field vector, and the longitudinally averaged components of the mean electromotive force.The integration over contour δΣ is performed clockwise, first, along the top boundary at r t = 0.99R from the equator to the North pole, then along the North axis to the bottom boundary at r i = 0.67R, and, finally, along the radius in the equatorial plane, back to the top boundary.A similar contour integration can be written for the Southern hemisphere flux, Φ S tor .Similarly to CS15, we estimate the RHS of Equation ( 15) in the coordinate system co-rotating with the solar equator with the angular velocity Ω (E) = 462 nHz : Here, U 0ϕ = R sin θΩ (E) is the rotational velocity, R ≡ R ⊙ , superscripts (E) and (N ) refer to the solar equator and North pole; superscripts (t) and (i) denote the top boundary level at r t = 0.99R and the inner boundary level at r i = 0.67R.The integral kernels, I 1 and I 2 , represent magnetic induction U×B calculated along the top boundary and the equatorial radius, respectively.The kernels, I 3 and I 4 , represent the contributions of the turbulent electromotive force.
The nonaxisymmetric components of the large-scale magnetic field and magnetic field of BMRs can contribute to this budget equation via the longitudinally averaged effects of the mean-electromotive force and effects of E (BM R) .
In the solar case, our model shows that the axisymmetric toroidal field is dominant inside the convection zone (P22; PKT23).Therefore, we neglect the direct effects of the nonaxisymmetric magnetic field and BMR's activity in Eq. ( 16).Injections of BMRs in the convection zone and the near-surface activity of the nonaxisymmetric magnetic field generate the poloidal magnetic field.These effects are included in the terms I 1 and I 2 .It is noteworthy that the first two terms of Eq.( 16)could be transformed using integration by part as follows, where A is the axisymmetric vector-potential, see the Eq(2.While the first line (Eq.17) is very useful for the numerical estimation, especially using observational data sets, it could be misinterpreted sometimes.The second line preserves the standard interpretation of the the differential rotation effect.Here, the term, −r t sin θA(r t , θ), represents (with the factor 2π) the net flux of the meridional component of the poloidal magnetic field across the radial section of the convection zone at polar angle θ.The second term rA(r, π 2 ) represents the net flux of the axisymmetric radial magnetic field through the northern hemisphere of radius r.
To estimate I 1 and I 2 , we use the mean angular velocity profile produced in our dynamo models.In this procedure, we neglect the effect of the torsional oscillations on I 1 and I 2 (associated with the large-scale dynamo).We normalize the kernels by dividing them on the solar radius.
Figure 1(a-c) shows the time-latitude diagrams of the integral kernel, I 1 , calculated for our dynamo models and for the observation data.To estimate I 1 from observational data, we use the same surface differential rotation profile as in the model.We use synoptic maps of the radial magnetic field obtained from the data set produced cooperatively by the NASA Goddard Space Flight Center and NOAO Space Environment Laboratory.This work utilizes the SOLIS data obtained by the National Solar Observatory (NSO) Integrated Synoptic Program (NISP).This data set is combined with the synoptic maps of the radial magnetic field from SDO/HMI (Sun et al. 2011).The dynamo models qualitatively agree with the observational data, as well as the results of CS15, regarding the magnitude and sign of I 1 .This agreement comes from the similarity of the axisymmetric magnetic field evolution in the dynamo models and observations.However, in our model, parameter I 1 is confined to a much narrower region closer to the poles than in the observations.This is because of a stronger concentration of the radial magnetic field to the poles in our model in comparison to the solar observations.
Figure 2 shows the latitudinal profiles of I 1 and I 4 and the radial profiles of I 2 and I 3 for model T0 for the magnetic cycle minimum in 1996.The models show a sharp poleward increase of I 1 (Figure 2(a)).This effect contributes to the axisymmetric toroidal magnetic field generation by the latitudinal shear in the bulk of the convection zone.
The term I 2 describes the axisymmetric toroidal magnetic field generation by the radial gradient of the angular velocity.This term has a maximum near the bottom of the convection zone, where it has same sign and nearly the same order of magnitude as I 1 at the northern pole.In the main part of the convection zone, it changes sign.The terms I 3 and I 4 include the effects of the diffusive loss and the loss of the magnetic flux due to magnetic buoyancy.Generally, the sign of these terms is opposite to the sign of the flux generation terms I 1 and I 2 .
Figure 3 illustrates the latitudinal profiles of I 1 for the minima of Cycles 22, 23, and 24 in our dynamo models and the observations.In all cases, the observations show a step-like increase of I 1 above 50 • latitude and almost uniform I 1 near the solar poles.The dynamo models show similar profiles, but the distribution of I 1 is not uniform near the poles.PKT23 found that the polar magnetic field in solar observations does not increase toward the poles as strongly as our models show.This can be a reason for nearly uniform I 1 near the solar poles for the data of observations.16)) for the cycle minimum in 1996.The linear threshold for the y-axis is at 100 G/yr/R.
Figure 4a shows the time evolution of the RHS contributions in Equation ( 16) for models T0 and S2, and for the observational data set.Notably, these variations in the Southern hemisphere have the opposite sign.Model S2 shows a larger magnitude of of I 1,2 than model T0.The magnitude of this increase is about 10 percents of the whole I 1,2 .The larger I 1,2 is due to the effects of BMRs activity on the magnitude of the axisymmetric poloidal magnetic field.Similar effect was found by P22.The models fail to reproduce the relatively low N I 1 dℓ of the solar cycle 24.Both models show that the integrals of I 1 and I 2 are about the same magnitude, and they evolve sin-phase.These results suggest an important role of the radial rotational shear in the axisymmetric mean-field dynamo.
Figure 4b shows the phase diagrams illustrating the toroidal flux generation rate (horizontal axis) and the loss rate (vertical axis) in our dynamo models.The model parameters are similar to those deduced by CS15 from solar observations.Some differences are because of the additional generation and loss terms included in our models.Also, the behavior of the radial magnetic field near the poles affects the magnetic flux budget considerably.
Figure 5 shows the axisymmetric toroidal magnetic field (images in the top panels), the poloidal magnetic field lines (the solid lines show the clockwise oriented field), and the corresponding toroidal field generation rate by the differential rotation, i.e., ∇ (Ω) × ∇ (r sin θA) (bottom panels) in the convection zone in model S2 during Cycle 23.The toroidal field generation rate reaches its maximum of about 2 kG/yr at mid-latitudes near the bottom of the convection zone.These images and the accompanying animation illustrate the dynamo wave propagation and the induction effects of the radial and latitudinal gradients of the angular velocity on the evolving large-scale magnetic field.In the lower part of the convection zone, the radial gradient of the angular velocity stretches the closed parts of the axisymmetric poloidal magnetic field lines.This results in the opposite signs of the toroidal flux generation rate in the lower and upper parts of the convection zone.Simultaneously, it produces the dynamo wave propagating along the lines of constant rotation according to the Parker-Yoshimura rule (Parker 1955;Yoshimura 1975).The latitudinal shear stretches the open poloidal field lines in the upper half of the dynamo domain.The positive sign of the α-effect in the Northern hemisphere results in a coordinated action of the radial and latitudinal shear in the upper part of the convection zone.It follows that the radial gradient of the solar differential rotation provides a substantial part of the toroidal flux inside the convection zone and affects the magnitude of the mean toroidal flux there.It is noteworthy that, similar to results presented in Fig. 2 and Fig. 4a effects of the differential rotation near the surface of the solar poles and near equator below the bottom of the convection zone have the same sign.The open poloidal flux and the latitudinal rotational shear control the magnitude of the subsequent cycle.Their properties affect the correlation of the polar magnetic field strength during the solar minima with the magnitude of the toroidal field dynamo wave of the subsequent magnetic cycle.

DISCUSSION AND CONCLUSIONS
The solar dynamo is responsible for the generation and decay of the axisymmetric toroidal and poloidal magnetic field fluxes on the Sun.Motivated by the results of Cameron & Schüssler (2015) (CS15), we analyzed the dynamo magnetic flux budget using solar synoptic observations and results of our mean-field dynamo models (Pipin et al. 2023).One of the hypotheses of CS15 is the dominant role of the surface latitudinal differential rotation in the generation   16) for the instances of the activity minima of Cycles 22, 23, and 24 in our dynamo models and the observations.The linear threshold for the y-axis is set at 100 G/yr/R. of the axisymmetric toroidal magnetic field of the Sun.Our results suggest that the radial rotational shear can be important as well, at least for the mean-field type distributed dynamo models.The study shows the generation rate of the toroidal field by the differential rotation in the 3D dynamo models, which include the emergence of Bipolar Magnetic Regions (BMRs), is higher than in the corresponding axisymmetric 2D models without BMRs.This is because the BMRs activity provides the additional source of the poloidal magnetic field generation for this type of dynamo (P22).
We find that the toroidal flux generation rate by the surface latitudinal differential rotation strongly depends on the radial magnetic field profile in the polar regions.Observations show that the radial magnetic field is almost uniform near the solar poles during epochs of the solar minima.In contrast, the dynamo models show a sharp increase in the toroidal flux generation rate toward the poles.One reason could be a strong effect of the meridional circulation on the near-surface radial magnetic field.The toroidal flux generation rate by the latitudinal differential rotation reaches its maximum at the poles in the observations and dynamo models during the solar activity minima.This property explains the relative success of the correlation between the polar magnetic field strength and the magnitude of the subsequent magnetic cycle for the solar cycle predictions (e.g., Schatten et al. 1978;Choudhuri et al. 2007).Our results suggest that this correlation may depend on the profile of the radial magnetic field near solar poles.The long-term measurements of the large-scale solar magnetic fields in the solar polar regions can help to resolve the issue.
where the first term describes the α-effect caused by the BMR tilt, and the second term models the magnetic buoyancy instability.The magnetic buoyancy velocity, V β , includes the turbulent and mean-field buoyancy effects (Kitchatinov & Rüdiger 1992;Kitchatinov & Pipin 1993;Ruediger & Brandenburg 1995): where, function ξ β , defines the location and formation of the unstable part of the magnetic field, the function H (β) describes effects of magnetic tension, where, β = |⟨B⟩| / 4πρu 2 c , (H (β) ∼ 1 for β ≪ 1) and subscript 'm' marks unstable points.To determine the radial position of the unstable point, r m , we use the maximum of product B I β (r, θ).Here, we restrict r to the upper part of the convection zone (above 0.85R), where I β reads, where B is the strength of the axisymmetric toroidal magnetic field and ρ is the density profile.For the case of ζ = 1, we get Parker's instability condition (Parker 1979).In our model, we use ζ = 1.2.In the unstable region, we must have I β > 0. We define ξ β as follows, where ψ is a kink-type function of radius, The latitudinal and longitudinal coordinates θ m and ϕ m in function ξ β are taken from the active region database.We relate the BMR areas to the parameter m β as follows, where S a is the maximum observed area (in the millionth of the hemisphere) of the bipolar active regions.In our computations, we exclude BMRs with S a < 50.The solar observations show a wide range of variations in the BMR emergence time, δt and growth rate, τ a .Also, there are periods of simultaneous emergence of several BMRs in one hemisphere.To avoid the overlaps, we reformat the emergence initiation time as follows.Firstly, for such cases, we shift the emergence of subsequent BMRs by two time steps after the end of the previous BMR emergence.Secondly, we define the minimal emergence time δt min = 2 days, (Weber et al. 2023), and assume that these BMRs have growth of τ 0 = 1 day.Specifically, we define: where δt a is the total emergence time of the BMRs, τ a is the growth rate.The parameters S a , θ m and ϕ m are taken from the NOAA database of solar active regions (https://www.swpc.noaa.gov/).The α-effect of the E (BMR) is given as follows Here, we put the amplitude of the α-effect to be determined by the local magnetic buoyancy velocity.The ξ α parameter controls the random fluctuation of the BMR tilt.The parameter, C αβ controls the mean amplitude of the BMR α-effect and tilt.To reproduce Joe's law, we put C αβ = 0.5 (PKT23).Note that, for C αβ = 0, the BMRs are oriented along the axisymmetric toroidal field.We model the randomness of the tilt using the parameter ξ α .Similar to Rempel (2005), the ξ α evolution follows the equation, Here, g is a Gaussian random number.It is renewed every time step, τ h .The τ ξ is the relaxation time of ξ α . .The parameters ξ 1,2,3 are introduced to get smooth variations of ξ α .Similarly to the above-cited papers, we choose the parameters of the Gaussian process as follows: g = 0, σ (g) = 1.We choose τ ξ = 2 months, which is roughly equal to the lifetime of big active regions (Norton et al. 2023).The ξ α varies independently in the northern and southern hemispheres.

B. Dynamo model parameters
Figure 6 illustrates distributions of the angular velocity, meridional circulation, the α -effect, and the eddy diffusivity in the nonmagnetic case model.The amplitude of the meridional circulation on the surface is about 13m/s.The angular velocity distribution is in agreement with the helioseismology data.

Figure 1 .
Figure 1.Contour lines in panels (a) and (b) show the time-latitude diagrams of the near-surface toroidal magnetic field in the range of ±1kG for Models T0 and S2; the background color images show the generation rate of the toroidal field by the latitudinal shear , i.e., integral kernel I1 of Equation (16); the green squares in panel (b) mark the latitudinal positions of the BMRs's injections; panel (c) shows the toroidal field generation rate calculated for the axisymmetric magnetic field using the Kitt Peak Observatory (KPO) and SDO/HMI synoptic maps, hereafter R ≡ R⊙.

Figure 2 .
Figure2.Estimation of contributions to the magnetic flux budget (Equation (16)) for the cycle minimum in 1996.The linear threshold for the y-axis is at 100 G/yr/R.

Figure 3 .
Figure 3. Panels (a), (b), and (c) show the latitudinal dependence of the surface term, I1, to the toroidal magnetic flux budget Equation (16) for the instances of the activity minima of Cycles 22, 23, and 24 in our dynamo models and the observations.The linear threshold for the y-axis is set at 100 G/yr/R.

Figure 4 .
Figure 4. (a) Time evolution of the integral RHS contributions of Equation (16) for the northern hemisphere for models S2 (solid lines) and T0 (dashed lines) and for the observational data set; (b) the phase diagram of the toroidal flux generation and loss rates (the horizontal axis shows the sum of π/2 0 I1rtdθ+ r t r i I2dr) and the vertical axis shows the sum of r t r i I3dr+ π/2 0 I4dθ) for model T0 (black curve) and model S2 (blue curve).

Figure 5 .
Figure 5. (a) The evolution of the axisymmetric toroidal magnetic field (color images) and the poloidal magnetic field streamlines (contours, solid lines show the clockwise directed magnetic field) during Cycle 23 in model S2, illustrating the dynamo wave propagation in the convection zone; the boundaries of the dynamo domain are shown by solid black lines, the green line marks the bottom of the convection zone; (b) the toroidal field generation rate by the differential rotation, i.e., ∇ (Ω) × ∇ (r sin θA) (The online version shows the animation of this Figure.The left and right columns of the animation illustrate variations of the axisymmetric magnetic field parameters (panels a) and the toroidal field generation rate (panels b) for Cycles 23 and 24 in Model S2 ).

Figure 6 Figure 7
Figure6illustrates distributions of the angular velocity, meridional circulation, the α -effect, and the eddy diffusivity in the nonmagnetic case model.The amplitude of the meridional circulation on the surface is about 13m/s.The angular velocity distribution is in agreement with the helioseismology data.Figure7(a) shows a snapshot of the radial magnetic field at the top of the dynamo domain for model S2 for the epoch of the solar maximum around the year 2000.Figure7(b) shows the time-latitude diagram of the axisymmetric magnetic field in the model S2.Further details can be found inPKT23.