A Babcock-Leighton dynamo model of the Sun incorporating toroidal flux loss and the helioseismically-inferred meridional flow

We investigate whether the Babcock-Leighton flux-transport dynamo model remains in agreement with observations if the meridional flow profile is taken from helioseismic inversions. Additionally, we investigate the effect of the loss of toroidal flux through the solar surface. We employ the 2D flux-transport BL dynamo framework. We use the helioseismically-inferred meridional flow profile, and include toroidal flux loss in a way that is consistent with the amount of poloidal flux generated by Joy's law. Our model does not impose a preference for emergences at low latitudes, but we do require that the model produces such a preference. We can find solutions in general agreement with observations, including the equatorward drift of the butterfly wings and the cycle's 11 year period. The most important free parameters in the model are the depth to which the radial turbulent pumping extends and the turbulent diffusivity in the lower half of the convection zone. We find that the pumping needs to extend to depths of about $0.80R_{\odot}$ and the bulk turbulent diffusivity needs to be around 10 km$^2$/s or less. We find that the emergences are restricted to low latitudes without the need to impose such a preference. The flux-transport BL model, incorporating the helioseismically inferred meridional flow and toroidal field loss term, is compatible with the properties of the observed butterfly diagram and with the observed toroidal loss rate. Reasonably tight constraints are placed on the remaining free parameters. The pumping needs to be to just below the depth corresponding to the location where the meridional flow changes direction. Our linear model does not however reproduce the observed"rush to the poles"of the diffuse surface radial field resulting from the decay of sunspots -- reproducing this might require the imposition of a preference for flux to emerge near the equator.


Introduction
The solar cycle is driven by a self-excited fluid dynamo which is induced by the interaction between the large-scale magnetic field and flows within the convection zone of the Sun (Ossendrijver 2003;Charbonneau 2014).In the first part of the dynamo loop, differential rotation winds up poloidal magnetic field generating toroidal field.This is the so-called Ω-effect.The Ω-effect is both well understood and constrained by observations.
In the second part of this loop, the toroidal field generates new poloidal magnetic field.This new poloidal field has the opposite polarity to the original poloidal field.Each 11year sunspot (or Schwabe) cycle is half of the 22-year magnetic (or Hale) cycle required to revert to the original polarity.Nonaxisymmetric flows and fields play a critical role during this phase of the cycle.The non-axisymmetric processes involved in the second phase, however, are far from being either well understood or constrained.
A major success of helioseismology was the determination of the sub-surface solar rotation profile, which challenged the dynamo-wave paradigm (Gough & Toomre 1991).This challenge lead to the flux-transport dynamo (FTD) model (Wang et al. 1991) where the deep meridional circulation causes the emergence locations of sunspots to drift equatorwards during a solar cycle.Observational and theoretical studies (see eg. Dasi-Espuig et al. 2010;Kitchatinov & Olemskoy 2011;Cameron & Schüssler 2015) provide strong support for the Babcock-Leighton mechanism (BL -Babcock 1961;Leighton 1964Leighton , 1969) ) to be the dominant mechanism in the second part of the dynamo loop, as opposed to the turbulent α-effect (Parker 1955;Steenbeck et al. 1966).At the core of the BL mechanism is the role played by the surface field.Sunspots emerge in bipolar magnetic pairs, with an east-west orientation usually in accordance with Hale's law.There is also a statistical tendency called Joy's law where the following spots to emerge closer to the poles and the leading spots to emerge closer to the equator.Sunspots decay within a few days to months, after which the field is dispersed by small-scale convective motions and transported poleward by the meridional flow.The flux cancellation of leading sunspot fields across the equator allows for the net buildup of a polar field by trailing sunspot fields.
Transport processes are required at the surface in order to transport the radial field from the equator to the poles, and to transport the subsurface toroidal field equatorwards to account for the equatorial migration of the butterfly wings (Spörer's law).
The surface part of the required transport has been established by observations of the surface meridional flow and the success of the surface flux transport model.The helioseismically inferred subsurface meridional flow is a relatively new constraint for these models.The use of the helioseismically-inferred meridional flow profile removes a number of free parameters from Babcock-Leighton type models.This makes a comparison with the observations a tighter test of the model and allows us to better constrain the remaining free parameters.
An additional recent constraint is that toroidal flux is lost through flux emergence, with a timescale estimated to be around 12 years by Cameron & Schüssler (2020, hereafter CS20).This paper will investigate if the BL FTD model, using the meridional flow inferred by Gizon et al. (2020, hereafter G20), and including the toroidal flux loss associated with flux emergence, is consistent with observations.To this end we introduce a loss term in the evolution equation for the toroidal field that is consistent with the evolution of the poloidal flux associated with Joy's law.
In the Babcock-Leighton type of model considered in this paper, the turbulent convective motions are not explicitly simulated, instead their affect on the magnetic field is parameterized (eg.Charbonneau 2014).Mean-field theory (Moffatt 1978;Krause & Rädler 1980) actually shows that including the effect of turbulence introduces a large number of parameters.In the Babcock-Leighton model only a few are kept, the most important of which include the α-effect, an increased turbulent diffusion, and downward diamagnetic pumping.Of these, the Babcock-Leighton α-effect is poorly understood but well constrained by observations, while turbulent pumping and turbulent diffusion are largely unconstrained.
In this paper we will see if the FTD model, with the observed meridional flow and flux loss, is compatible with the Babcock-Leighton model, and what constraints it places on the other processes of the model.

Dynamo equations
In mean-field theory, the axisymmetric large-scale magnetic and velocity fields are decomposed into poloidal and toroidal components as: where A is the ϕ-component of the vector potential field and B is the toroidal component of the large-scale magnetic field, u m is the meridional circulation, and Ω is the differential rotation.For an isotropic turbulent diffusivity that has only a r-dependence, the kinematic mean-field dynamo equations are: where ϖ = r sin θ, u p = u m + γ, η is the turbulent diffusivity, γ the turbulent pumping, S the BL source term, and L the toroidal field loss term due to flux emergence.The source and loss terms will be discussed in Section 2.4.Belvedere et al. (2000) given by equation 5 (left) and cycle-averaged and symmetrized stream function of the helioseismic meridional flow inversions of G20 (right).For the latter, positive values represent clockwise circulation and negative anticlockwise.The dash-dotted and dotted lines represent the approximate locations of the tachocline at 0.7R ⊙ and reversal of the meridional flow direction at 0.8R ⊙ , respectively.

Differential rotation and meridional circulation
The large-scale flows need to be prescribed in kinematic models.For the differential rotation, we use the simple model provided by Belvedere et al. (2000) and shown in the left panel of Figure 1: where the coefficients c i j can be found in that paper.This fit is an approximation of the helioseismologically inferred rotation rate of Schou et al. (1998).
As mentioned in the introduction, we will use the meridional circulation inferred from observations.In this study, we use the inversions of G20.The authors furnish the meridional flow for cycles 23 and 24.In order to keep the parameter space study manageable, we take the average of the two cycles.In addition, since we are not here interested in the asymmetry between both hemispheres, we also symmetrize the flow across the equator.The right panel of Figure 1 shows the meridional circulation we use in all our models.Note that the flow switches from poleward to equatorward at a radius of about 0.785R ⊙ , which we will call the meridional flow turnover depth r t , and the region beneath it the lower or deep convection zone, and the one above the upper or shallow convection zone.

Turbulent parameterizations
We choose a turbulent diffusivity profile which is written as a double step as in Muñoz-Jaramillo et al. (2011): where η RZ = 0.1 km 2 /s, η R ⊙ = 350 km 2 /s, and η CZ are respectively the radiative core, surface, and bulk values of the turbu-lent diffusivity.η CZ is a free parameter and η R ⊙ has been chosen to be consistent with estimates from observations (eg.Komm et al. 1995), the surface flux transport model of Lemerle et al. (2015), andMLT (eg. Muñoz-Jaramillo et al. 2011).The values in the first error function have been chosen so that the drop in diffusivity occurs mostly before the helioseismically determined position of the tachocline (Charbonneau et al. 1999), roughly coinciding with the overshoot region (Christensen-Dalsgaard et al. 2011).
For the turbulent pumping, we adopt the profile given by Karak & Cameron (2016, hereafter KC16 -see also the discussion in their section 2): where here we take r γ = r t = 0.785R ⊙ .This choice will be discussed in Section 4.2.3.

Flux emergence
The emergence of bipolar magnetic regions both removes toroidal flux (see CS20) and creates poloidal flux (because of Joy's law).These two processes are clearly linked as are the respective loss and source terms in Equations 4 and 3.
We take the amount of flux emerging to be proportional to the toroidal flux density where the integration is over the depth of the convection zone.This prescription corresponds to a dynamo where the toroidal field is not necessarily stored near the tachocline but can be distributed throughout the convection zone (KC16, Zhang & Jiang 2022).It is in part motivated by observations of dynamos in fully convective stars (Wright & Drake 2016), and by cyclic dynamo action in 3D MHD simulations of spherical shells without a tachocline (eg.Brown et al. 2010;Nelson et al. 2013Nelson et al. , 2014)).
The flux emergence rate R can be written in general as a function of latitude: where τ 0 is a timescale and f θ (θ) is its latitudinal dependence, which we take to be corresponding to an emergence probability that is constant per unit length of the toroidal field lines.
In general the timescale, τ 0 in Equation 9 depends on the dynamics associated with flux emergence.If these dynamics are dominated by the large-scale field than the buoyant rise time and τ 0 to depend inversely on the mean-field value of B 2 (Kichatinov & Pipin 1993).If however small-scale magnetic fields remain coherent over timescales longer than the correlation time, then the B filling factor which can be far from 1 becomes important.Some previous studies (eg.Schmitt & Schüssler 1989;Moss et al. 1990a,b;Jennings & Weiss 1991) assume τ ∼ B −2 so that the loss term scales like B 3 .In this paper we consider the linear case where τ 0 is a constant, and which would correspond to a case where the field is composed of filamentary structures with lifetimes longer than the turnover timescale of the turbulence and with local field strengths drawn from some distribution which is independent of flux.We stress that the aim in this paper is to consider a simple linear system.We defer the nonlinear case to future work.
The orientation of the flux emergence is governed by Joy's law which states that the leading polarity flux emerges on average closer to the equator than the trailing polarity one.We take the form of Joy's law used in Leighton (1969), sin δ = 1 2 cos θ, where δ is the angle between the solar equator and the line joining the two polarities.Then the rate at which toroidal flux density is lost due to flux emergence is where the subscript L indicates the contribution from the loss term.The tilting of a toroidal flux tube as it emerges gives rise to a θ-component of the same polarity.The rate at which this θ-component of the flux density is lost is where the subscript S indicates the contribution from the source term S .This is what gives rise to the Babcock-Leighton mechanism.From the definition of the poloidal field (Equation 1), We choose a depth R b sufficiently below the base of the convection zone so that the 11-year cyclic component of the field is negligible there.Then multiplying both sides by r and integrating from the base of the convection zone to the surface, we obtain Therefore, in terms of the ϕ-component of the poloidal vector potential A, Equation 12becomes Next, we need to prescribe the radial structure of the source and loss terms.For the source term S , we follow KC16 and assume where r S is the depth to which the source extends.For most of the calculations, we choose r S = 0.85R ⊙ (as in Muñoz-Jaramillo et al. 2010).There are indications that this disconnection should happen much deeper than the usually assumed shallow location of 0.95R ⊙ (Longcope & Choudhuri 2002).We will nevertheless vary this parameter in order to study its impact on the solutions.
For the loss term L, we first distribute the toroidal flux density loss (equation 11) over radius according to the amount of flux initially present, so that where In order for the differential form of the source and loss terms to be valid, the emergence timescale τ e must be considered infinitesimal with respect to the timescale over which the magnetic configuration of the large-scale field changes appreciably.It follows that a timescale separation must hold: In the case of the Sun, τ e ∼ 1 day, and P = 11 years, so the timescale separation is reasonable.This formulation implicitly ignores the effect of the meridional flow on the the emergence process.Nevertheless, this formulation is sufficient to study the general features of the solar cycle.

Numerical procedure
Equations 3 and 4 are nondimensionalized and numerically solved in the meridional plane with 0 ≤ θ ≤ π and 0.65R ⊙ ≤ r ≤ R ⊙ .We use a spatial resolution of 241 × 301 evenly spaced in radius and colatitute grid points and a time step of 5 × 10 −6 R 2 ⊙ /η t , where η t = 10 km 2 /s.The inner boundary matches to a perfect conductor, so that: and the outer boundary condition is radial: The latter is necessary for FTD models to match surface flux transport models (Cameron et al. 2012).A second-order centered finite difference discretization is used for the spatial variables and the solution is forwarded in time with the ADI scheme (Press et al. 1986).We use the code initially developed by D. Schmitt in Göttingen (as also used by Cameron et al. 2012).The linearity of equations 3 and 4 allows us to choose τ 0 and γ 0 such that the dynamo is approximately critical (σ ≤ 5 × 10 −5 per year) with a cycle period of 12 years (within 0.1%), roughly the average period of cycles 23 and 24.This way we reduce our parameter space to only one dimension (η CZ ).Our model being linear also means we can arbitrarily scale A and B. In order to facilitate comparisons with observations, we scale the fields so that the maximum of the surface radial field is 10 G, which is consistent with the observed polar field strengths at cycle minimum (eg.Hathaway 2015).

Toroidal flux loss timescale
The general expression for the toroidal flux decay timescale is where Φ is understood as the net subsurface toroidal flux in the northern hemisphere: and dΦ dt is its decay.In order to calculate the latter, we first need to specify the toroidal flux loss mechanism.For the loss L due to flux emergence, we have and its corresponding timescale will be denoted by τ L .Toroidal flux is also lost due to the explicit diffusion across the solar surface: with an associated loss timescale of τ η .The observational constraint on the flux loss timescale is given by CS20, and corresponds to τ L ≈ 12 years at solar maximum.Since these timescales vary across a cycle, we will calculate them at cycle maximum (to be defined in Section 3.2).It is possible to estimate the range of values τ 0 should take.To do so we assume the tilt angle δ associated with Joy's law is small, so that cos δ ∼ 1, and that the toroidal flux density can be approximated by b(θ, t) = b 0 (t) sin m θ, where b 0 (t) is a time-dependent scalar, and m determines how closely the field is concentrated near the equator.With these approximations, we obtain: where τ 0 is the timescale for flux loss in the model, defined in equation 9, and the limits correspond to m = 0 and m = ∞.
The model parameter τ 0 should thus be comparable to (not more than a factor of 2 smaller than) the observed toroidal flux loss timescale τ L .

Polar cap and activity belt flux densities
An important observational constrain is that the maxima of azimuthally averaged polar flux densities should be around the same strength as maximum flux densities in the activity belt.
Since the evolution equations we are using are linear, we have nominally set the maximum of the polar field to be 10 G.This implies that the azimuthaly averaged radial field in the butterfly wings in the model should also be around 10 G.This is a constraint on the model which we will return to when evaluating whether the model can produce solar-like cycles.

Cycle phase of polar maxima
An important constraint we will take into account is when the maximum of polar flux occurs.As it is observed to happen quite close to the activity minimum, its corresponding phase shift of about 90 • with respect to cycle maximum.To measure this shift in the simulations, we need a definition for when the cycle maxima and the maxima of polar fields occur.
We begin with the surface radial flux of the polar cap and of the activity belt We then define cycle maximum times, T a , as the times when the activity belt flux Φ a is maximum.The times of the maximum polar flux T p is similarly defined.T a and T p are both defined by the signed fluxes, and hence each have one maxima per magnetic cycle (of about 24 years).
For each cycle, i, we then calculate the phase shift between the polar maximum times and the activity maximum times ∆ϕ by where P is the (activity) cycle period.

Results
We first present our reference model including the Babcock-Leighton loss term in Section 4.1, then examine parameter sensitivity (Sections 4.2 to 4.5), and finally gauge the importance of the dynamo wave (Section 4.6).

Reference model
Our reference model has a bulk diffusivity of η CZ = 10 km 2 /s, which places the simulation in the advection-dominated regime (an explanation of why such a low diffusivity is required in our setup is given in Section 4.2.2).Values of τ 0 = 12.5 years and γ 0 = 43.8m/s were required to achieve a critical 12-year (activity) cycle period dynamo.The reference model is presented in Figure 2. It can be seen that we are able to produce a reasonably solar-like butterfly diagram when using the helioseismicallyinferred meridional flow of G20.

Comparison with observational constraints
This model has an emergence loss timescale of τ L = 17.2 years, somewhat longer than the 12 years inferred by CS20.The diffusive loss timescale, on the other hand, is τ η = 96.8years.
The estimate of CS20 is based on all toroidal flux escaping through the photosphere.The combined rate at which flux is lost through the photosphere in the model can be estimated to be 1/(1/τ L + 1/τ η ) = 14.6 years and hence close to the inferred 12-year timescale.The value of ∆ϕ from the model is noticeably larger than the observed value at ∆ϕ = 134°.The maximum value of the surface radial field in the butterfly wings is around 5 G, so about half of the maximum polar field strength (The observed average polar field is similar to the average in the butterfly wings, see e.g. the butterfly diagrams in Hathaway 2015).
Considering the polar field strength is somewhat uncertain, our results are not inconsistent with observations.Our simulations do not have the problem of very large polar fields typical of FTD models (Charbonneau 2020).
The net toroidal flux Φ (lower panel of Figure 3) shows that it reaches a maximum value of about 5 × 10 23 Mx close to cycle maximum, which is in rather good agreement with the estimates of Cameron & Schüssler (2015) for cycles 22 and 23.
An important result of our model is that it achieves confinement of emergences to the low observed latitudes without the need for an emergence probability decreasing faster than sin θ with latitude.This can be seen in the upper panel of Figure 2, where we see that the toroidal flux density is mainly strong near the equator.The left panel of Figure 3 shows the toroidal field is mainly stored deep in the convection zone.This confinement of the toroidal field to deep in the convection zone is a consequence of the imposed radial pumping.The confinement to near the equator is then due to the meridional flow which advects the material from high latitudes towards the equator.The combination of radial pumping and equatorward meridional flow in the lower half of the convection zone leads to a stagnation point near the equator in the lower half of the convection zone where the field builds up until it is removed through emergence (also see Cameron & Schüssler 2017;Jiang et al. 2013).
Our simulated butterfly diagram differs from the observed one in that it lacks a distinct "rush to the poles" of the trailing diffuse field of the decayed sunspots.

Parameter dependence
Our model has 5 free parameters: the source and loss terms timescale τ 0 and the depth where sunspots are anchored r S , the turbulent pumping amplitude γ 0 and the depth it reaches down to r γ , and the turbulent diffusivity in the bulk of the convection zone η CZ .In this section we will first provide a qualitative description of what different choices of the parameters produce.
Importantly, the results and constraints we find are for under the assumption that f θ (θ) = sin θ, i.e. that there is no imposed preference for emergences to occur at low latitudes.We also performed simulations with f θ (θ) = sin 12 θ (as in KC16), and as expected were are able to find critical dynamo solutions which match the observations for a much wider range of parameters.

Influence of the source depth
We here investigate how the choice of r S affects the solutions.We used the same value of η CZ as in reference case, and varied r S .The values of τ 0 and γ 0 were then chosen so that the growth rate is zero and the cycle period is 12 years.We found that the solutions with flux loss are not very dependent on r S for r S extending from just above 0.78 to the surface.This is because the solutions with flux loss require strong pumping which rapidly stretches the poloidal field so that they extend radially to the depth at which the pumping stops.This makes the model insensitive to the initial depth of the poloidal source term in (at least in the region of parameter space near the reference case).

The bulk diffusivity
We only find critical 12-year periodic solutions when the bulk diffusivity is of the order of 10 km 2 /s.This is a consequence of the strong radial shear of the equatorward component of the helioseismically-inferred meridional flow u θ in the lower third of the convection zone.The strong radial shear leads to toroidal flux at different depths being advected in latitude at very different rates.This implies that flux originally concentrated at one latitude over a range of depths will be quickly spread out in latitude.
To understand the role of the radial shear of the latitudinal flow, we can imagine toroidal field initially at one latitude but spread out in radius from r = 0.766R ⊙ to 0.785R ⊙ .These depths were chosen so that the meridional flow will vary from almost 1 m/s equatorwards to almost 0 m/s.Over 5 years this will spread the flux over a latitudinal band of 157 Mm.This spreading out would be similar to a diffusivity of (157) 2 /5 Mm 2 /year = 156 km 2 /s.In the context of Babcock-Leighton FTD models, this is a large value.As a comparison, the 1D model of Cameron & Schüssler (2017) which also assumes f θ (θ) = sin θ, requires the latitudinal diffusivity in the bulk of the convection zone to be lower than about 100 km 2 /s.This effective diffusivity is, however, in agreement with the estimate of of 150 − 450 km 2 /s of Cameron & Schüssler (2016) based on the properties of the declining phase of solar cycle.The essential point of the above is that, unless the toroidal field is confined to a narrow range of depths, the latitudinal shear in the meridional velocity quickly spreads the toroidal field out in latitude.Consequently, if there is no imposed preference for emerging near the equator then the butterfly diagram ceases to be solar-like.The requirement for confinement in latitude is what imposes the constraint that η CZ ≈ 10 km 2 /s.We consider this fixed for the rest of this paper.We also comment that if the radial shear in the differential rotation was weaker, then this constraint would be much weaker.

Turbulent magnetic pumping
With our chosen f θ (θ) = sin θ (Eq.10), we find growing dynamo solutions only for values of r γ not far away from 0.785R ⊙ .This depth is where the helioseismically-inferred meridional flow profile changes direction from poleward to equatorward,  and roughly where numerical simulations suggest the convection zone might be weakly subadiabatic (Hotta 2017).A slightly broader range of pumping depths can be achieved when the loss term is not included, but shallower depths make the butterfly wings broader.
Our conclusion from this is that the pumping depth is fairly tightly constrained if the appearance of spots to low latitudes is only caused by the equatorward meridional flow leading to a build up at low latitudes.The depth of the pumping is poorly constrained if the preference for low latitude emergence is imposed.Observations do not provide estimates for the amplitude of turbulent pumping at depth, and so it is interesting to compare our results with those from global MHD simulations.Shimada et al. (2022) find that in the outer half of the convection zone has a turbulent diffusivity of ∼ 10 km 2 /s, similar to the values we needed in our model without an imposed preference for emergence at low latitudes, they also find that γ r peaks at ∼ 10 m/s, while Simard et al. (2016) and Warnecke et al. (2018) find amplitudes of the order of 1 − 2 m/s or half the root-mean-square velocity.The latter is of the order of 40 m/s according to mixinglength (Vitense 1953;Böhm-Vitense 1958) estimates.It thus appears the pumping velocities required in the reference case are too large by a factor of 2 to 3. We defer a discussion of this to Sections 4.4 and 4.5.

Role of the toroidal flux loss term L
In order to gauge the importance of emergence loss term, in this subsection we switch off the term by setting L = 0 in Equation 4. We first use the same parameters as in the reference case.The resulting cycle period is shorter at 11.3 years.However, the solution is rapidly growing, with a growth rate of about 66% per cycle.This growth is not unexpected, as emergences are no longer able to remove the subsurface toroidal flux and it must now be removed either through its "unwinding" by the new cycle flux, or by diffusive cancellation across the equator.Because emergences no longer deplete the subsurface toroidal flux, more poloidal field is generated so that the polar fields are reversed much faster, explaining the shorter period.
We also investigated the 12 year period critical solutions when L = 0. Doing so required τ 0 = 9.3 years (as against the reference case where τ 0 = 12.5 years), and a turbulent pumping γ 0 = 6.41 m/s.The time-latitude diagrams and meridional cuts for this case are shown in figures 4 and 5.The surface dif-fusion loss timescale is reduced to τ η = 25.6 years, only a factor of two larger than the observed value of the toroidal flux loss timescale.This is due to the lowered pumping amplitude, which makes the diffusion of the toroidal field through the surface less difficult than in the case with strong pumping.Even with such low pumping, the poloidal field near the surface is still almost radial (Figure 5).
Looking at the butterfly diagram, the most apparent difference is the large decrease of the polar field strengths compared to the fields in the butterfly wings.The maximum value of the latter goes up from about 5 to 7.5 G, which is relatively close to the observed value of around 10 G.In this case, we find maxima phase between the polar field maxima and active region maxima is 101 • , similar to that which is observed.
Note that the pumping amplitude of γ 0 = 6.41 m/s in this model is in much better agreement with estimates from global MHD simulations (cf.Section 4.2.3).This is because models without the loss term achieve shorter periods much more easily.In principle, then, very large pumping velocities are not necessary to obtain a functioning dynamo for this class of models.

Sensitivity to the meridional flow and differential rotation
Here we investigate how sensitive the simulations are to the meridional flow and differential rotation.We do this by considering the inferred meridional flow for cycles 23 and 24 separately, and a differential rotation profile which differs significantly from the one of the reference model at high latitudes.As is also the case for the reference solution, we do not impose a preference for emergences at low latitudes (if we impose a preference for emergences at low latitudes, then the parameter space where the model has similar properties to the observations becomes much larger).As all solutions mentioned in this section have qualitatively the same butterfly diagrams as the reference case they are not shown.
First, we consider individually the symmetrized (across the equator) meridional flow profiles of cycles 23 and 24.Using the meridional flow from cycle 23, a 12-year periodic critical dynamo requires γ 0 = 16.3 m/s.This is a substantial reduction from the reference case.Increasing the period to 13.3 years and keeping the criticality requirement led to pumping speeds of γ 0 = 10 m/s.Using the meridional flow from cycle 24, we were unable to find critical solutions with periods shorter than 12 years.A critical solution with γ 0 = 10 m/s required a period of 16.6 years.
Clearly the model, where emergence is not restricted to low latitudes, is very sensitive to the meridional flow.This is because emergences at high latitudes are inefficient at getting flux across the equator, which is what eventually reverses the polar fields.The observational constraint, that the cycle period is similar to the timescale for which toroidal field is lost through the surface due to flux emergence, implies that the cycle period involves a balance between the flux transport to low latitudes and the loss through emergence.
In this context, both mean-field theory (eg.Kichatinov 1991;Kitchatinov & Nepomnyashchikh 2016) and global numerical models (eg.Shimada et al. 2022, and references therein) indicate equatorial latitudinal turbulent pumping could also be substantial.From the BL-FTD modelling this would correspond to an increase in the meridional return flow, and would lead to a reduction in the strength of the required radial pumping.
Second, we consider the sensitivity to differential rotation.Again we consider critical 12-year periodic solutions, using the average meridional flow profile used in the reference case.We now use the differential rotation profile of Larson & Schou (2018), which differs from that of the reference at high latitudes (rather than the analytic fit of Belvedere et al. (2000) often used in dynamo studies and used in the reference case).The required parameter values are τ 0 = 11.8 yrs and γ 0 = 28.5 m/s.τ 0 is again of the order of cycle period and the 12-year estimate for the toroidal flux loss timescale of Cameron & Schüssler (2020).But the pumping velocity is much more reasonable.

Sensitivity to assumption that growth rate is zero and period is 12 years
We have thus far concentrated on the kinematic regime with zero growth rates.The Sun is certainly in a statistically saturated state.
The kinematic case with growth rate zero is relevant if the system is weakly nonlinear.Whether or not this is the case for the Sun is open (for arguments in favour of this see van Saders et al. 2016;Metcalfe et al. 2016;Kitchatinov & Nepomnyashchikh 2017).In the strongly nonlinear case, the period will be substantially affected by the choice of the nonlinearity and the growth rate in the linear regime is no longer a constraint.The observations are thus less constraining in the strongly nonlinear case.For this reason, we have focused on the weakly nonlinear case and have looked for zero-growth rate solutions to the linear problem.The addition of a weak nonlinearity will slightly modify both the growth rate (in the saturated state it will be zero) and period.Hence in this section we consider the sensitivity of the growth rate and periods to τ 0 and γ 0 .
Figure 6 shows the cycle period and the growth of the toroidal flux per cycle as a function of the timescale parameter τ 0 .As in KC16, we observe that increasing the source term amplitude τ −1 0 causes the growth rate to increase until the cycle period becomes too short for the meridional flow to transport the field (see Section 4.1 of KC16).Eventually, the dynamo shuts down completely.Growing solutions can nonetheless be reached by further increasing the source term amplitude τ −1 0 .However, the resulting cycle periods are very short (≲ 3 years) and the dynamo is now driven by a dynamo wave propagating equatorwards in the high-latitude tachocline.
The effect of the pumping amplitude on the growth rate and cycle period is shown in Figure 7.The growth rate is very sensitive to the pumping amplitude at lower values, where the operating threshold is not yet met, as flux emergence then quickly removes the toroidal flux at high latitudes.But the effect of  pumping saturates as its amplitude increases.At some point the time required for the poloidal field to reach the lower convection zone is essentially instantaneous.Note that we have concentrated on solutions near the bifurcation point where dynamo action switches on.This very likely makes the model more sensitive to the different parameters than would be the case if we were considering a non-linear, saturated dynamo.

Role of the dynamo wave
To investigate what role the subsurface meridional return flow is playing, we apply the same procedure as KC16 to our reference model, namely we switch off the equatorward component of the meridional flow (see Section 3.2 of KC16 for a discussion).Figure 8 shows the resulting magnetic field butterfly diagram.For the reference parameters this mode is decaying (and so is not a dynamo) with fields almost entirely located above 45 • .Since there is no equatorward component of the meridional flow, the equatorward migration of the field is due to the negative radial rotation shear in the high-latitude tachocline, and the direction of propagation is in accordance with the Parker-Yoshimura sign rule.We hence, not surprisingly, conclude that the subsurface meridional flow is essential in this model.

Conclusion
Using the helioseismically-inferred meridional flow of G20, we have shown that the Babcock-Leighton FTD model remains gen-erally consistent with observations.We have also shown that the long-standing problem of the latitudinal distribution of sunspots can be solved if turbulent pumping reaches depths just under 0.80R ⊙ , but not much deeper, where the meridional flow's direction switches from poleward to equatorward.High turbulent pumping velocities are necessary to essentially store the toroidal flux under this location, in agreement with the results of G20 (see also Parker 1987).There, the meridional flow, in conjunction with the Ω-effect through the latitudinal shear present in the bulk of the convection zone (which is maximal at mid-latitudes), causes an accumulation of toroidal flux at equatorial latitudes.Turbulent pumping effectively short-circuits the meridional circulation, preventing significant generation of toroidal field at high latitudes.No additional restriction of emergences to low latitudes is required.
Our model using the helioseismically inferred meridional flow, and including the observed toroidal flux loss associated with flux emergence in a way that is consistent with the Babcock-Leighton source term, is able to reproduce the observed properties of the solar cycle, including the latitudinal migration of the sunspot wings and the approximately 11 year period.Our reference model predicts a toroidal flux loss timescale of 14.8 years at cycle maximum, compared to the estimate of 12 years of CS20.

Fig. 1 .
Fig. 1.Rotation profile ofBelvedere et al. (2000) given by equation 5 (left) and cycle-averaged and symmetrized stream function of the helioseismic meridional flow inversions of G20 (right).For the latter, positive values represent clockwise circulation and negative anticlockwise.The dash-dotted and dotted lines represent the approximate locations of the tachocline at 0.7R ⊙ and reversal of the meridional flow direction at 0.8R ⊙ , respectively.

Fig. 2 .
Fig. 2. Time-latitude diagrams of the toroidal flux density b (top), surface radial field B r (R ⊙ ) (middle), and the net toroidal flux in the northern hemisphere Φ (bottom) for our reference model.The vertical dotted lines indicate the times where the snapshots of Figure 3 were taken.

Fig. 3 .
Fig. 3. Meridional cuts of the North hemisphere toroidal field (left column) and poloidal field (as r sin θA, right column) of the reference model for specific times indicated by the vertical dotted lines in Figure 2. The dotted lines are located at radii of 0.95, 0.85, and 0.80R ⊙ , the approximate locations of the bottom of the near-surface shear layer, r S , and r γ respectively.

Fig. 4 .
Fig. 4. Time-latitude diagrams of the toroidal flux density b (top), surface radial field B r (R ⊙ ) (middle), and the net toroidal flux in the northern hemisphere Φ (bottom) for the case without the flux loss associated with flux emergence (and with a period of 12 years and zero growth rate).The vertical dotted lines indicate the times where the snapshots of Figure 5 were taken.

Fig. 5 .
Fig. 5. Meridional cuts of the North hemisphere toroidal field (left column) and poloidal field (as r sin θA, right column) of the the case without the flux loss associated with flux emergence (and with a period of 12 years and zero growth rate) for specific times indicated by the vertical dotted lines in Figure 4.The dotted lines are located at radii of 0.95, 0.85, and 0.80R ⊙ , the approximate locations of the bottom of the near-surface shear layer, r S , and r γ respectively.

Fig. 6 .
Fig. 6.Percental per-cycle growth of the toroidal flux (solid line, left axis) and cycle period (dashed line, right axis) as a function of the timescale parameter τ 0 .

Fig. 7 .Fig. 8 .
Fig. 7. Percental per-cycle growth of the toroidal flux (solid line, left axis) and cycle period (dashed line, right axis) as a function of the pumping amplitude γ 0 .