Turbulent Processes and Mean-Field Dynamo

Mean-field dynamo theory has important applications in solar physics and galactic magnetism. We discuss some of the many turbulence effects relevant to the generation of large-scale magnetic fields in the solar convection zone. The mean-field description is then used to illustrate the physics of the α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha $\end{document} effect, turbulent pumping, turbulent magnetic diffusivity, and other effects on a modern solar dynamo model. We also discuss how turbulence transport coefficients are derived from local simulations of convection and then used in mean-field models.


Introduction
The problem of solar and stellar dynamos is still an open one. In spite of tremendous progress over recent decades, we still do not understand with any degree of certainty the reason behind the equatorward migration of solar activity belts, the dependence of cycle frequency on rotation frequency, or the level of magnetic activity.1 All models of solar and stellar magnetism rely on some assumptions. Even the most realistic simulations suffer from finite resolution and the compromises in the physics that are made. The crucial question is then, when and where we are allowed to make compromises and when not. Among those approximations is the secondorder correlation approximation (SOCA), also known as the first-order smoothing approximation. These are nowadays either replaced by other approximations or by numerical techniques such as the test-field method (TFM), as will be explained later in this review.
The Sun's magnetic field exhibits a clear mean field with spatio-temporal order: antisymmetry of radial and toroidal fields about the equator and the 11-yr cycle. This mean field can well be described by an azimuthal average. The radial component of such an azimuthally averaged mean field has a typical strength of ±10 G. This is not much compared with the peak strength of ±2 kG in sunspots, but much of this is "lost" in the process of averaging. Of course, whatever is lost corresponds to fluctuations, which actually play crucial parts and correlations between different fluctuations lead to various mean-field effects.
Mathematically, once an averaging procedure has been defined, we have the mean field , indicated by an overbar. Then, the difference between the actual and the mean field, and , gives the fluctuating field as ≡ − . The same procedure also applies to all other quantities. This formal distinction between mean and fluctuating fields, which are sometimes also called large-scale and small-scale fields, is important in discussions with observers. Coronal mass ejections, for example, are superficially reported as being part of a large-scale field, but this may not be true anymore when we think of averaging over the full solar circumference. Thus, paradoxically, even if something is large by some standards, it may not qualify as large-scale under this formal definition of an azimuthal average.
Azimuthal averaging is not always a good recipe. Some stars have nonaxisymmetric magnetic fields, and even the Sun is believed to have what is known as active longitudes -a weak nonaxisymmetric magnetic field on top of a predominantly axisymmetric one. Those nonaxisymmetric fields might best be described through low-order Fourier mode filtering. This is probably completely fine, but slightly problematic at the formal level, because then the average of the product of mean and fluctuating fields is no longer vanishing, as it would be in the case of an azimuthal average. This mathematical property is one of several rules that are called the Reynolds rules. However, as alluded to above, the violation of this particular Reynolds rule is probably just a technicality that makes mean-field predictions less accurate. We refer here to the work of Zhou et al (2018) for a detailed investigation. There are a number of other limitations in mean-field theories that will be discussed below.
The purpose of defining mean fields is twofold. On the one hand, they allow us to quantify large-scale magnetic, velocity, and other fields that are observed or that are present in a simulation. On the other hand, they allow us to develop predictive theories for these averages. In these theories, mean fields can sometimes emerge because of instabilities and/or because of suitable boundary conditions. This is possible because of certain mean-field effects, by which one usually means the relations between correlations of fluctuations and various mean fields. Discussing those effects is an important purpose of this review. The ultimate goal of mean-field dynamo theory is to understand and model the Sun and other stars. We therefore also discuss in this review the status of such attempts. For a basic introduction to mean-field theory, which is not the subject of this review, we refer to standard textbooks (Moffatt, 1978;Krause and Rädler, 1980;Zeldovich et al, 1983) and other reviews (Brandenburg and Subramanian, 2005a;Kulsrud and Zweibel, 2008;Miesch and Toomre, 2009;Charbonneau, 2014Charbonneau, , 2020Brandenburg, 2018a;Tobias, 2021;Brandenburg and Ntormousi, 2023;Karak, 2023).

The golden years of dynamo theory
The first mean-field model was constructed by Parker (1955). In his model, the toroidal magnetic field is generated from the dipole field by nonuniform rotation. To overcome the restrictions of Cowling's theorem (Cowling, 1933), Parker suggested that the dipole magnetic field can be regenerated by cyclonic convective motions which transform emerging toroidal magnetic loops into poloidal magnetic field. The coalescing loops can amplify the dipole magnetic field. Studying the combined action of differential rotation and cyclonic motions, he found a solution in the form of a dynamo wave and formulated conditions for the equatorward propagation of dynamo waves. Steenbeck et al (1966) and Steenbeck and Krause (1969) constructed the theoretical basis of mean-field theory, introduced the notion of the mean electromotive force (MEMF) of the turbulence and showed that the Parker effect of the cyclonic convective motions is equivalent to the effective MEMF along the large-scale field. The 1970s can be considered the golden years of mean-field dynamo theory. Back then, Schuessler (1983) stated: "dynamo theory reached the textbook state", mentioning the famous monographs by Moffatt (1978), Parker (1979), Krause and Rädler (1980), Vainshtein et al (1980), andZeldovich et al (1983).
Indeed, intensive theoretical and observational studies led to the establishment of the basic solar dynamo scenario, identification of key dynamo parameters, and the formulation of a general paradigm of the nature of solar and stellar magnetism. Schuessler (1983) summarized that mean-field dynamo models can reproduce the "physics of solar activity to a great extent" including: • Hale's polarity rule of sunspots groups, • the time-latitude evolution of sunspot activity ("butterfly diagram"), • reversals of the polar magnetic field, • the phase relationship between the evolution of poloidal and toroidal magnetic fields and their consistence with the observed butterfly diagrams (Stix, 1976), • rigid rotation of the magnetic sector structure and coronal holes (Stix, 1974(Stix, , 1977, • chaotic variations of dynamo activity either due to a random effect or the dynamo nonlinearity from the Lorentz force (Leighton, 1969;Yoshimura, 1978;Ruzmaikin, 1981), and • a quantitative understanding of the solar torsional oscillations (Schuessler, 1981;Yoshimura, 1981).
We note that the first and second items are based on the assumption that sunspot groups are formed from the large-scale toroidal magnetic field. Already at the time it was recognized that the mean-field models need to take into account the fibril state of the magnetic field which we observed at the solar surface. We return to this point later.
Classical mean-field dynamo models utilize the Ω scenario using the differential rotation (Ω effect) as the source of the toroidal magnetic flux production and the effect for the poloidal magnetic field generation. Since the seminal work of Pouquet et al (1976), it started to become clear that the magnetic helicity results in an important nonlinear contribution to the effect and turbulent magnetic field generation (Kleeorin and Ruzmaikin, 1982).

Mean-field theory and avoiding some of its limitations
We can never expect a mean-field theory to produce an accurate representation of reality. One reason is the fact that the underlying turbulence has stochastic aspects, so each realization with slightly different initial conditions would result in a somewhat different outcome. However, there could be other reasons for discrepancies that we discuss next. Some of those discrepancies can nowadays be avoided.

Mean-field electrodynamics
In mean-field theory, one derives evolution equations for the averaged fields, namely the mean magnetic field , the mean velocity , and the mean thermodynamic variables such as mean specific entropy and the mean density . Often, one neglects the evolution of , , and , which is then already an important limitation.
If one focuses on the evolution of the mean magnetic field only, one often talks about the mean-fields electrodynamics or quasi-kinematic mean-field theory, which can still be nonlinear if the various mean-field transport coefficients depend on the mean fields. If they are unaffected, one talks about kinematic mean-field theory, which is linear. Of course, once there is a dynamo, we have an exponentially growing solution, so the magnetic field would grow without limit, i.e., it would not saturate within kinematic mean-field theory. Obviously, a correct mean-field theory must be nonlinear, but even within the realm of linear theory, there are important lessons to be learnt. Below, we discuss the aspects of nonlocality, which were often omitted out of ignorance, but nowadays we know that this is often not possible and this restriction can easily be relaxed.

Nonlocality
The mean magnetic field is governed by the mean induction equation, which is sometimes also referred to as the mean-field dynamo equation. The most important term here is the mean electromotive force, i.e., the averaged cross product of velocity and magnetic fluctuations. In mean-field electrodynamics, it is often expressed as where the ellipsis denotes higher derivative terms, of which there should be infinitely many, and there should also be time derivatives. The term E 0 is a contribution that can exist already in the absence of a mean field; see Brandenburg and Rädler (2013) for details and numerical experiments. Including only a finite number of derivatives in Eq.
(2) and ignoring time derivatives is another important approximation. In fact, it is usually easier to express E as a convolution between an integral kernel and the mean field. Furthermore, it is instructive to split the integral kernel into two pieces and write E = E 0 +ˆ * +ˆ * / , where the asterisks mean a convolution in space and time, and the hats denote integration kernels. In principle, the spatial derivative can be absorbed as being part of the integral kernel, but separating the kernel intoˆandˆhas conceptual advantages, because they preserve the similarity to Eq. (2). Note also that, unlike Eq.
(2), where we allowed for arbitrarily many derivatives, here, we have no other terms, because all even derivatives are already absorbed inˆand all odd derivatives are absorbed inˆ. Time derivatives can also be absorbed in both of them if the convolution with the kernels is also over time.
For the benefit of better interpretation, both and (and analogously also forˆandˆ) can be broken down into further pieces. The tensor can be split into a symmetric and an antisymmetric tensor. The latter is characterized by a vector, = − 1 2 , which corresponds to a pumping velocity. Having in mind that the magnetic gradient tensor can also be split into symmetric and antisymmetric parts, where the latter is the mean current density, , with = − 1 2 / , we can separate the rank-3 tensor, , into a rank-2 tensor operating only on and the rest operating on the symmetric part of / . The convolution can only be replaced by a multiplication, as in Eq. (2), if the mean field is constant in time (which is normally never the case!) and if it varies at most linearly in space (which is normally also not the case). We return to this point further below.

Beyond SOCA and scale separation limits
An important question concerns the calculation of the and coefficients or kernels. A problem arises from the fact that the differential equations for these expressions are nonlinear and therefore hard to solve analytically. A commonly used approximation is SOCA. It neglects triple (and higher) correlations in the evolution equations for the fluctuating velocity and magnetic fields. This closure can be applied when either the magnetic Reynolds number, Re M , or the Strouhal number, St, are much smaller than unity. These limits are rather restrictive for astrophysical conditions. For example, the convection zones (CZs) of the Sun and stars are in a turbulent state with huge values of the fluid Reynolds number (Re ≳ 10 12 ), the magnetic Reynolds number (Re M ≳ 10 8 ), Rayleigh number (Ra ≳ 10 20 ), and an extremely low Prandtl number (Pr ∼ 10 −4 -10 −7 ); see, e.g., Ossendrijver (2003).
Results of Schrinner et al (2005) showed that SOCA does not work well when Re M exceeds unity and St is not small. There are analytical approaches, e.g., different variants of the so-called approximation (Kleeorin et al, 1996;Field and Blackman, 2002;Brandenburg and Subramanian, 2005a), which can be applied in the high Reynolds number limit. The restrictions inherent to SOCA or the approximation no longer apply when calculating solutions of the underlying differential equations numerically. This is done in the TFM (Schrinner et al, 2005(Schrinner et al, , 2007. Using a set of mean magnetic fields, the TFM allows one to determine the turbulent transport coefficients for arbitrary velocity fields, provided they can be computed or otherwise represented on the computer. The velocity field can be determined either as a solution of the nonlinear Navier-Stokes equations for a forced turbulent flow or it can be obtained as a results of global convective dynamo (GCD) simulations.
To compute E, the solution for the induction equation for the fluctuating magnetic field is needed as well. The original TFM of Schrinner et al (2005Schrinner et al ( , 2007 adopts the scale-separation assumption. It was shown that the TFM describes the dynamo processes for GCD simulations at moderate Reynolds numbers of around 50 rather well (Schrinner, 2011;Warnecke et al, 2018;Viviani et al, 2018). The calculations within TFM have some technical restrictions and are currently unable to meet the very high astrophysical limits of Re, Re M > 10 6 . Nevertheless, the current applications of the TFM concern cases with Re, Re M ≫ 1, which is well beyond the SOCA limits. In recent developments of the TFM, Käpylä et al (2022) took compressibility effects into account. They also studied the effects of the small-scale dynamo on the turbulent electromotive force; see also Rempel et al (2023).
An alternative way of extracting the coefficients of the mean-electromotive force employs a multi-dimensional regression method (Brandenburg and Sokoloff, 2002;Racine et al, 2011;Augustson et al, 2015;Simard et al, 2016). In this approach, instead of solving the equations for the fluctuations in the presence of different mean magnetic fields, as it is done in the TFM, the regression methods try to exploit the form of Eq. (2) for the dynamo-generated large-scale magnetic field. Detailed comparisons of the above method with TFM were done by Warnecke et al (2018). It was found that the TFM gives a more accurate representation of the mean-field coefficients than the multidimensional regression method. We encourage the reader to consult this paper for further details. We return to the problem of extracting turbulent transport coefficients from GCDs in Sect. 7.
The limitations discussed so far are in principle all avoidable: (i) The evolution equations for , , and can be (and have been) included , in addition to that for , but in practice, even this is still an approximation in the sense that the full set of coefficients for the equations is not (or only approximately) known. (ii) The electromotive force can be (and has been) expressed as a convolution, which can most effectively be solved by rewriting the equations as a differential equation, as will be described below. (iii) Numerical solutions can be employed to find specific values for and ; see Warnecke (2018); Warnecke et al (2021) for doing this for solar simulations using the TFM. It often turns out that analytical closure techniques are very useful as a first orientation and they are often also accurate enough for a qualitatively useful model. In special cases, when a more accurate solution is required, the answer may well be obtained numerically using the TFM. The problem is then only that numerical solutions themselves are limited in just the same way as those for a full numerical solution in the solar and stellar dynamo problems. Figure 1 shows results for˜( ) and˜t ( ) with / = 1. Both˜and˜t decrease monotonously with increasing | |. The functions˜( ) and˜t ( ) are well represented by Lorentzian fits of the form (4) Also shown in the lower part of Fig. 1 are the kernelsˆ( ) andˆt ( ), obtained through Fourier transforms of the Lorentzian fits, We see that the profile ofˆt is half as wide as that ofˆ, but it is not known whether this is a general property. It is important to realize that the suggested mean-field modifications employing the Lorentzian forms of the integral convolution kernels are based on empirical results. Nevertheless, they are much more accurate than the approximation of replacing the kernels by functions, which is done in conventional approaches.
Under suitable conditions, the accuracy of the TFM can be so high that discrepancies become apparent that are solely the result of having made unjustified approximations in the comparison. An example is the memory effect. Comparing the growth rate for a supercritical dynamo with that obtained theoretically from the coefficient obtained from the TFM can give noticeable discrepancies if the memory effect is neglected; see Fig. 1 of . The combined Fourier transformed integral kernel is of the form where is well approximated by the turbulent turnover time. Even for stationary flows, the memory effect can be dramatically important (Rädler et al, 2011). In practice, it is cumbersome to solve the integral equation in time. However, as alluded to above, it is possible to approximate this integral equation by a differential equation for E with respect to space and time of the form This has been done in several papers (Rheinhardt and Brandenburg, 2012;Rheinhardt et al, 2014;Brandenburg and Chatterjee, 2018;Pipin, 2023). We return to this in Sect. 5.3.

The use of mean-field theory
If mean-field theory cannot reliably be applied to a regime outside that of the direct numerical simulations (DNS), one must ask what is then the use of mean-field theory. The answer lies in the fact that mean-field theory provides us with an excellent diagnostic "tool" for approaching the problem. Particular features of a solution can usually be attributed to particular terms in the mean-field equations. This would then allow us a more informed answer by saying that the main dynamo mechanism is, for example, of Ω type, or of a specific type of a shear flow dynamo, for example. Thus, mean-field theory may be regarded as a convenient tool for understanding what is going on rather than predicting what might be going on.

The catastrophic quenching problem
Since the 1990s, a problem emerged in that numerical dynamo solutions were found to depend on the value of the microphysical magnetic diffusivity. Typically, the strength of the mean-fields then decreases with increasing magnetic Reynolds number. This is unusual and does not have any correspondence with ordinary hydrodynamics where the large-scale dynamics is usually already captured at moderate fluid Reynolds numbers. In its original form, the catastrophic quenching problem refers to the finding that the volume-averaged electromotive force scales with the microphysical magnetic diffusivity, and thus goes to zero when → 0. To some extent, this is a problem related to the use of periodic boundary conditions. However, even for astrophysically more realistic boundary conditions, numerical simulations reveal that there is still a problem. Plasma relaxation experiments have identified the role of magnetic helicity as the main culprit in causing -dependent large-scale dynamics and catastrophic quenching. We therefore begin by briefly reviewing the essential findings.

Lessons from plasma relaxation experiments
The magnetic field is divergence-free and can therefore be written as = ∇ × , where is the magnetic vector potential. The magnetic helicity density is defined as · . Its evolution equation follows directly from the uncurled induction equation, / = − − ∇Φ, or, using Ohm's law, − = × − 0 , so It yields the evolution equation for the magnetic helicity density, It must here be emphasized that there is an important difference to the equation for the magnetic energy density, While both Eqs. (9) and (10) have analogous terms such as dissipation ∝ · versus ∝ 2 , respectively, and flux terms × versus Poynting vector × / 0 , respectively, there is the work against the Lorentz force, − · ( × ) in Eq. (10), which would be · ( × ) in Eq. (9), but it obviously vanishes. In statistical equilibrium, ⟨2 0 2 ⟩ must balance −⟨ · ( × )⟩, which implies that the current density diverges like | | ∼ −1/2 . By contrast, no magnetic helicity is being produced, and also its dissipation converges to zero like ∝ | · | → 1/2 as → 0.
Already since the 1970s, we know of the conjecture of J. B. Taylor (1974Taylor ( , 1986) that the magnetic field relaxes under the constraint of magnetic helicity conservation to a nearly force-free field to minimize dissipation. The approximate conservation of magnetic helicity has been verified experimentally in plasma relaxation experiments; see, e.g., Ji et al (1995). There are obviously some differences between the solar convection zone and laboratory plasmas, for example, the role of the electron pressure in the generalized Ohm's law could play an important role in explaining why magnetic helicity changes are observed to be faster in plasma experiments than what is predicted by Ohm's law (Ji et al, 1995). Furthermore, the effect has been identified as the main agent for converting magnetic helicity from the turbulent field to the mean field (Ji, 1999).
In the context of plasma relaxation experiments, it is useful to distinguish between electromagnetic and electrostatic turbulence. This distinction refers to the curl-free and divergence-free parts of the electric field written as = −∇Φ − / . In plasma relaxation experiments, turbulence is mostly electrostatic. It can be affected by the electron pressure gradient ( e ) −1 ∇ e in the generalized Ohm's law, where is the unit charge, and e and e are the electron density and electron pressure, respectively; see Ji (1999) for details. This leads to a battery term ∝ ∇ e × ∇ e in the equation for / and to a magnetic helicity flux, which transports magnetic helicity across physical space, as opposed to wavenumber space (Ji, 1999).
There is the possibility that the divergence of a magnetic helicity flux, F f , itself can constitute an effect. This corresponds to = − 1 2 −2 ∇ · F f ; see Vishniac and Cho (2001), who have derived a specific form for such a flux. The subscript 'f' indicates that the flux originates from correlations of the fluctuating magnetic field. Mean-field models of the type described below have shown that a dynamo can operate even without kinetic helicity, i.e., it is based only on shear and current helicity fluxes, provided a nondimensional scaling factor in front of the magnetic helicity flux exceeds a certain critical value (Brandenburg and Subramanian, 2005c). However, there are so far no DNS that have supported this kind of behavior, nor has the proposed flux been confirmed (Hubbard and Brandenburg, 2012). Nevertheless, the idea of an effect being related to the magnetic helicity flux divergence is certainly consistent with the laboratory experiment presented in Fig. 1 of Ji (1999).
The effect reflects the physics of the inverse cascade of magnetic helicity (Pouquet et al, 1976). In the absence of energy input, this is known to lead to a slower turbulent decay of magnetic energy ∝ −2/3 (Hatori, 1984;Biskamp and Müller, 1999). In hydrodynamics, by comparison, the kinetic energy density decays like −10/7 or −6/5 , depending on the initial subinertial range energy spectrum (Davidson, 2000); see also Brandenburg and Larsson (2023) for a comparison with the magnetic case.
In the presence of magnetic driving by applying a voltage drop along the magnetic field, small-scale instabilities such as the tearing instability develop. This leads to a sawtooth-like time dependence in the mean toroidal magnetic flux; see Ji and Prager (2002) for a review. This can be associated with the resulting development of a mean electromotive force, E = × , along with an effect that accomplishes the helicity transport (Ji, 1999).
Unlike astrophysical dynamos, which are generally understood as self-excited ones, the plasma experiments operate in a regime where a magnetic field is always present, but it is then redistributed by the effect. The conceptional similarities and differences have been discussed in detail by Blackman and Ji (2006). In the following, we discuss in more detail the consequences imposed by magnetic helicity evolution in astrophysical dynamos. It is important to emphasize, however, that the same physics that is used to explain the catastrophic quenching phenomenology also applies to plasma experiments such as the reversed field pinch, as was shown in corresponding mean-field simulations by Kemel et al (2011).

Mean fields in periodic domains
Under astrophysical conditions of interest, is so small that the volume-averaged electromotive force would be negligibly small. If this result was actually astrophysically relevant, it would be a "catastrophe," i.e., it would not be possible to understand astrophysical magnetic fields as mean-field dynamos. The solution to this particular problem turned out to be that relating the volume-averaged electromotive force to the volume-averaged mean magnetic field is only of limited relevance to the problem of effect dynamos. Any dynamo would produce a non-uniform field. Especially in a periodic domain, the mean magnetic flux through any of the faces of the periodic domain is constant in time, so if it was zero to begin with, it would always remain zero. A dynamo problem can therefore not be formulated in that case.
A proper dynamo problem should always allow for the possibility of the magnetic field to decay to zero if there is sufficient magnetic diffusivity. Simple examples of nontrivial mean fields in a periodic domain are Beltrami fields of the form which can be solutions of the simple 2 dynamo problem, / = ∇× + T ∇ 2 . Nevertheless, there is still a problem of catastrophic nature because it turned out that the time required to reach the final solution scales inversely with . This is demonstrated in Fig. 2, where we show the evolution of one of the three planar averages. In the beginning, all three mean fields grow in a similar fashion, but at some point, only one of the three reaches a significant amplitude. Note, however, that the ultimate saturation takes a resistive time, res = 1/(2 2 1 ).

Quenching phenomenology
To understand the reason for the catastrophically slow saturation, it suffices to consider Eq. (9) for the magnetic helicity density. For periodic domains, we just have This equation is gauge-independent, because the gauge transformation → ′ +∇Λ yields ⟨ · ⟩ = ⟨ ′ · ⟩, with ⟨ · ∇Λ⟩ = ⟨∇ · (Λ )⟩ − ⟨Λ∇ · ⟩ = 0, using ∇ · = 0 and the fact that the domain is periodic, so the average of a divergence vanishes. In numerical calculations, it is often convenient to adopt the Weyl gauge, which implies that / = − , i.e., the scalar potential term drops out. For fully helical large-scale and small-scale magnetic fields of opposite magnetic helicity, Eq. (12) becomes (Brandenburg, 2001) with the solution This agrees with the slow saturation behavior seen first in the simulations of Brandenburg (2001); see Fig. 2. Here, sat is the time when the slow saturation phase commences; see the crossing of the green dashed line with the abscissa. Interestingly, instead of waiting until full saturation is accomplished, one can obtain the saturation value already much earlier simply by differentiating the simulation data to compute (Candelaresi and Brandenburg, 2013) 2 Since res involves the microphysical magnetic diffusivity, the quenching is still in that sense catastrophic.

The quenching formula
A more complete description is in terms of kinetic and magnetic effects, i.e., and observing the fact that the magnetic helicity evolution of averages and fluctuations is given by Equation (17) allows for the possibility that magnetic helicity can be produced by the mean electromotive force, because, in general, E · ≡ × · ≠ 0. (By contrast, of course, ( × ) · = 0.) In particular, if E = − t 0 , then, E · = 2 − t 0 · , which produces positive (negative) magnetic helicity of the mean field when > 0 ( < 0) Equation (18) is constructed such that the sum of Eqs. (17) and (18) yields Eq. (12). Given that ⟨ · ⟩ is related to ⟨ · ⟩, which, in turn, is related to a magnetic contribution to the effect (Pouquet et al, 1976), Eq. (18) can be rewritten as an evolution equation for the total (Brandenburg, 2008), which can also be expressed in the form where "extra terms" Note that the last term is here a time derivative. Equation (20) resembles the catastrophic quenching formula of Vainshtein and Cattaneo (1992), but it also shows that it needs to be extended in several important ways: when the mean field is no longer defined as a volume average, extra terms emerge that are of the same order as those in the denominator. They can therefore potentially offset the catastrophic quenching.
In practice, this is only partially true, because there are also other terms, for example the aforementioned time derivative term. It is responsible for the fact that a strong field state is only reached after a resistively long time.

Magnetic helicity fluxes and helicity reversals
Magnetic helicity fluxes could in principle remove the catastrophic quenching problem, but only if preferentially small-scale magnetic helicity is being removed (Kleeorin et al, 2000). To see this, let us first consider the problem of an 2 dynamo in insulating boundaries using the Weyl gauge, i.e., The boundary condition implies that = = 0, and is therefore also referred to as the vertical field condition. In this 1-D problem, however, this boundary condition is equivalent to a proper vacuum boundary condition.
The 2 dynamo with this boundary condition was first considered by Gruzinov and Diamond (1994), who found that the saturation field strength of such a dynamo decreases with Re M . This was later confirmed by Brandenburg and Dobler (2001). In Fig. 3, we show the profiles of magnetic helicity, current helicity, and the magnetic helicity fluxes for Runs A of Brandenburg (2018b) with Re M = 180. The computational domain is 0 ≤ ≤ /2 with a perfect conductor boundary condition on = 0 and a vertical field condition on = /2; see Brandenburg (2017) for the relevant mean-field models. For normalization purposes, he defined the reference values He emphasized that the largest contribution to the magnetic helicity density comes from the large-scale field. Near the surface ( = /2), the (negative) magnetic helicity flux from small-scale fields is only about 0.02 m0 , which explains why they are not efficient enough to alleviate the catastrophic dependence of the resulting mean magnetic field (Del Sordo et al, 2013;Rincon, 2021). Subsequent simulations with an outside corona indicated that the magnetic helicity changes sign at or near the outer surface ). This was just a speculation and needs to be reconsidered with the help of global models of the type considered by Warnecke et al (2011Warnecke et al ( , 2012 and Brandenburg et al (2017a). This is shown in Fig. 4, where we present the line-of-sight averaged current helicity density, ⟨ · ⟩ in the plane of the sky using a simulation of Brandenburg et al (2017a). The quantity ⟨ · ⟩ is a proxy of magnetic helicity at small scales and shows clearly the reversal of sign between the dynamo interior and the exterior.

Radial magnetic helicity reversal in the solar wind
If the idea of alleviating catastrophic quenching by magnetic helicity fluxes is to make sense, we would expect to see signs of the expelled magnetic helicity in the solar wind. The magnetic helicity spectrum can be measured in the solar wind by determining the parity-odd contribution to the magnetic correlation tensor, which, in Fourier space, takes the form This would allow one to compute ( ) = Im(˜˜ * ) and ( ) = 1 2 (|˜| 2 + |˜| 2 ), which also obeys the realizability condition | ( )| ≤ 2 ( ).
The Ulysses spacecraft was the only one to cover high heliographic latitudes, where a non-vanishing sign of magnetic helicity can be expected. It turned out that ( ) has, as expected from dynamo theory, different signs in the northern and southern hemispheres. It also has different signs at small and large wavenumbers. This, in itself, is also expected from an 2 dynamo, because the effect produces no net magnetic helicity, but it separates magnetic helicity in wavenumber space. However, the signs are opposite to what is seen at the solar surface, where the helicity in the north is negative at small length scales. In the solar wind, however, it is positive in the north and at small scales. Of course, the meaning of small is here relative and has to be with respect to larger scales, where a sign change in has been seen. If one just assumed a linear expansion of all scales from the solar surface (radius = 700 Mm, to the location of the Earth at 1 AU, we expect a corresponding expansion ratio so that a wavenumber of 1 Mm −1 corresponds to 1/700 AU −1 . In particular, 20 Mm −1 corresponds to 2/70 AU −1 , which is close to the wavenumber where we see a sign-change in Fig. 5. It is unexpected, however, that at the solar surface ( Fig. 5b), the sign in the northern hemisphere changes from minus to plus as increases, while in the solar wind, it changes from plus to minus. This apparent mismatch may not just be a measurement error, but it may actually be a real result and would tell us that the simpleminded picture of expelling magnetic helicity of one sign all the way to infinity may not be accurate.
When the domain is inhomogeneous, and especially when there are boundaries, magnetic helicity fluxes are possible and Eq. (18) takes the form In the steady state, we have In the dynamo interior at the northern hemisphere, > 0, and, assuming 2 to dominate the EMF, we expect −2E · to be negative. However, a negative flux divergence of a negative quantity would eventually make this quantity positive, which is what has been observed. Whether or not this is really the right interpretation remains still an open question. It would clearly be useful to have an independent assessment of this interpretation.

Nonlocal effects of E and catastrophic quenching
Catastrophic quenching in large-scale dynamos is a rather general property. It is a consequence of the build-up of magnetic helicity of the mean magnetic field. It has been conjectured that catastrophic quenching would be prevented if the sources of toroidal field generation are spatially separated from the sources of the poloidal field; see, e.g., Tobias and Weiss (2007). This would be the case in what is known as interface dynamos (Parker, 1993). It could also be through a nonlocal effect. Such a nonlocal effect is an essential ingredient of the Babcock-Leighton and fluxtransport dynamo models; see, Hazra et al (2023). The studies of Brandenburg and Käpylä (2007) and Chatterjee et al (2010) showed that a spatial separation between shear and effects does in general not help to avoid catastrophic quenching for such types of dynamo models. It is interesting, however, that Kitchatinov and Olemskoy (2011) and later also Brandenburg and Hubbard (2015) found that the inclusion of diamagnetic downward pumping of the toroidal magnetic field can alleviate the catastrophic quenching in the Babcock-Leighton dynamo model with a strongly nonlocal .
The catastrophic quenching models are reasonably well reproduced by DNS when the geometries of the setups are sufficiently simple. It would therefore be worthwhile to apply DNS to conditions where turbulent pumping and a strongly nonlocal mean electromotive force can be expected. At present, however, even just the physical reality of a nonlocal effect of Babcock-Leighton type through the decay of active regions rests mainly on the interpretation of solar observations. Turbulence simulations have so far not been able to make contact with such concepts.

Alternative large-scale dynamo effects
Given the difficulties encountered with effect dynamos, there have been various attempts to construct large-scale dynamos that are not based on the effect. A common misconception here is the idea that catastrophic quenching would not apply if just because there is no effect. This is not true, because an M term can always emerge regardless of whether or not there existed an original effect. An example is the shear-current effect. It is due to the presence of shear and boundaries that a helicity can be introduced. Shear of the form = (0, , 0) implies a finite vorticity, ∇× = (0, 0, ) and boundaries would lead to a gradient vector of turbulent intensity near the boundaries. Thus, while there can be hope that catastrophic quenching may not be as strong, this may turn out not to be the case. An example of this was presented in Brandenburg and Subramanian (2005c).

Rädler and shear-current effects
The Rädler effect is another large-scale dynamo effect (Rädler, 1969). In the simplest representation it leads to an EMF proportional to × . It is similar to the shearcurrent effect. In this case it cannot change the magnetic energy of the mean field. Indeed, the energy equation for the mean field is given by In the general case, the generation effects due to global rotation and mean currents can be written as follows (see Krause and Rädler, 1980;Rädler et al, 2003;: where the coefficients 1,2,3 depend on the spatial profiles of the turbulent parameters such as the typical convective turnover time, the convective velocity rms , etc. The last two terms in this equation may lead to an 2 dynamo . For the solar case, the effect can provide an additional non-helical source of poloidal magnetic field generation. Interestingly, Pipin and Seehafer (2009) found that for the solar-type dynamos, i.e., those with equatorward propagation of the dynamo waves, the dynamo effect does not dominate the contributions of the -effect. We will discuss the available scenario in the next section.

Dynamos from negative turbulent magnetic diffusivity
There are two other effects that are noteworthy, although it is not clear that either of them can play a role in stellar convection zones. One is the negative turbulent magnetic diffusivity and the other is the memory effect in conjunction with a pumping effect.
When modeling a negative turbulent magnetic diffusivity dynamo, high wavenumbers must not be destabilized at the same time. Brandenburg and Chen (2020) studied classes of dynamos with a very low critical Re M . The Willis dynamo (Willis, 2012) has a critical Re M of 2.01, which is small compared to 6.3 for the Roberts flow and 17.9 for the ABC flow. In this dynamo, one of the two horizontally averaged field components grows exponentially, because the total magnetic diffusivity in that direction is negative (Brandenburg and Chen, 2020). The other component decays and is not coupled to the former one.
As we see from Fig. 6, t is negative only for < ∼ 1.5. The dependence of the turbulent magnetic diffusivity can be expanded up to second order as where the tildes indicate Fourier transformed quantities. In the proximity of = 1, which corresponds to the largest scale in the computational domain of 2 , we havẽ (0) ≈ −0.233 and˜( 2) ≈ 0.11. In addition, there is still the microphysical magnetic diffusivity, which is positive ( = 0.403). To a first approximation, one can just consider the equation for , which can then be written as We recall that the minus sign in front of the fourth derivative corresponds to positive diffusion if˜( 2) is positive, and so does the plus sign in front of the second derivative, unless the term in squared brackets is negative, which is the case we are considering here.

Dynamos from pumping and memory effects
Pumping effects alone cannot usually lead to interesting dynamo effects, unless there is also a memory effect. This effect means that the mean electromotive force depends not just on the instantaneous mean magnetic field at that time, but also on the mean magnetic field at earlier times. It is therefore described as a convolution between a pumping kernel and the mean magnetic field. This can lead to dynamo action, as has been demonstrated by Rheinhardt et al (2014) for the case of two of the four flow fields studied by Roberts (1972). These are flows II and III with II ( , ) and III ( , ), respectively, and are given by where is the wavenumber of the flow. Both flows have zero kinetic helicity and no effect, but flow II is also pointwise nonhelical. A supercritical three-dimensional magnetic field with growth rate and wavenumber in the direction of the form = 0 ( , ) +i is possible when Re M ≡ 0 / > 4.58 and 2.9 for flows II and III, respectively; see Rheinhardt et al (2014). Here, 0 ( , ) is the eigenfunction.
For both flows, there are -averaged mean fields ( , ) and ( , ), with waves traveling in opposite directions for flow II and in the same direction for flow III; see Figures 6 and 8 of Rheinhardt et al (2014), respectively. These dynamos appear to be atypical, because there is so far no other known example of a flow where pumping produces a memory effect that is strong enough to lead to dynamo action. This may well be due to the absence (until recently) of computational tools for determining the memory effect. Indeed, it was only with the development of the TFM (Schrinner et al, 2005(Schrinner et al, , 2007 that the importance of the memory effect was noticed ) and applied to pumping.
The dispersion relation for a problem with turbulent pumping and turbulent magnetic diffusion t is given by = −i − t 2 . Since Re < 0, the solution can only decay, but it is oscillating with the frequency = Im = . In the presence of a memory effect, is replaced by /(1 − i ), where is the memory time. Then, , and Re can be positive. This is the case for the Roberts flows II and III.
We return to nonlocality and memory effects further below in this article when we discuss concrete solar models; see Pipin (2023). One of the most obvious consequences of the memory effect is a lowering of the critical excitation conditions for the dynamo, which was already reported by Rheinhardt and Brandenburg (2012). Interestingly, for the nonlocal mean electromotive force, the lowering of the critical threshold can be accompanied by multiple instabilities of different dynamo modes that have different frequencies and spatial localization; see Pipin (2023).

Dynamos from cross-helicity
An alignment of velocity and magnetic field, i.e., cross helicity, plays a key role in numerous processes and phenomena of astrophysical plasmas. Krause and Rädler (1980) showed that the saturation stage of the turbulent generation is characterized by an alignment of the turbulent convective velocity and the magnetic field. This consideration does not account for the effects of cross-helicity that take place in the strongly stratified subsurface layers of the stellar convective envelope. For example, the direct numerical simulations of Matthaeus et al (2008) showed a directional alignment of velocity and magnetic field fluctuations in the presence of gradients of either pressure or kinetic energy.
The mean electromotive force in this case is along to the mean vorticity, where, Υ = u (0) · b (0) is the cross helicity pseudoscalar, and is the turbulent turnover time. The superscripts (0) indicate quantities of the background turbulence, which exists in the absence of a mean magnetic field and a mean flow; see our comment after Eq. (2) about the E 0 term. In the standard mean-field framework, it is assumed that Υ = 0; see Krause and Rädler (1980). Yoshizawa and Yokoi (1993) generalized the framework assuming Υ ≠ 0, see the comprehensive review of Yokoi (2013).
Dynamo scenarios based on cross helicity have been suggested in a number of papers (Yoshizawa and Yokoi, 1993;Yoshizawa et al, 2000;Yokoi, 2013). Pipin and Yokoi (2018) showed that the large-scale dynamo instability does not require the existence of a global axisymmetric mean. The mix of axisymmetric and nonaxisymmetric magnetic fields can be produced even in the case Υ = 0, where the overbar means azimuthal averaging. The surface magnetic field of the Sun and other similar stars tends to be organized in sunspots, plagues, ephemeral regions, super-granular magnetic network, etc. These structures tend to demonstrate the alignment of local velocity and magnetic fields (Rüdiger et al, 2011). Therefore, the cross helicity dynamo instability can contribute to dynamo generation effects that operate near the stellar surface. Stellar observations, for example those of Katsova et al (2021), require such dynamo effects to be working in situ at the stellar surface. The solar analogs show an increase of the spottiness with an increase of the rotation rate (Berdyugina, 2005). In that case, cross helicity dynamo effects can be considered as a relevant addition to the standard turbulent generation by means of convective helical motions. Rapidly rotating M-dwarfs show the highest level of magnetic activity (Kochukhov, 2021). There is a population of rapidly rotating M-dwarfs that show a rather strong dipole type magnetic field. These stars show a rather small level of differential rotation. For solid body rotation, an 2 dynamo generates a nonaxisymmetric magnetic field (Chabrier and Küker, 2006;Elstner and Rüdiger, 2007). At high rotation rates, the effect is highly anisotropic (Rüdiger and Kitchatinov, 1993). It cannot employ the component of the large-scale magnetic field along the rotation axis for the generation of an axial electromotive force. Results of Pipin and Yokoi (2018) show that the 2 Υ 2 scenario can produce a strong constant dipole magnetic field. The model predicts the existence of large-scale cross helicity patterns occupying the stellar surface. We hope that this can be tested either in observations or in GCDs.
The nonlinear theory for the cross helicity effect is not yet developed. Sur and Brandenburg (2009) showed that the turbulent generation due to Υ is quenched by large-scale vorticity in a way that is similar to catastrophic quenching given by Eq. (20), i.e., One should remember that for the initialization of the cross-helicity dynamo instability we have to seed both the cross helicity and the magnetic field. The solar type model scenarios based on cross helicity require an effect, which produces poloidal magnetic field and cross helicity at the top of the dynamo domain (Yokoi et al, 2016). Given that cross helicity is an ideal invariant of the MHD equations, it is natural to ask whether systems with strong cross helicity exhibit inverse cascading. The answer seems to be yes; see Brandenburg et al (2014). In Fig. 7 we demonstrate the gradual build-up of magnetic fields in the vertical direction when the system has significant cross helicity owing to the presence of a magnetic field along the direction of gravity (Rüdiger et al, 2011).

Mean-field dynamo models
In general, the effect, as well as any other turbulent generation effect, including the effect (Rädler, 1969), the shear-current effect (Kleeorin et al, 2000), and the crosshelicity effect (Yokoi, 2013) can generate both toroidal and poloidal magnetic fields. Therefore there can be a number of possibilities for solar-types dynamo models Krause and Rädler (1980); Yokoi et al (2016); Pipin and Kosovichev (2018). Some of them skip the effect altogether. For example, Seehafer and Pipin (2009) studied Ω Ω and Ω scenarios, where turbulent generation of the poloidal magnetic field is due to × and shear-current effect, respectively. These scenarios show oscillating solutions and the correct time-latitude diagram of the toroidal magnetic field if the meridional circulation is included. A similar possibility was mentioned earlier by Krause and Rädler (1980) for the Ω scenario. However, the given scenarios result in an incorrect phase relation between activity of the toroidal and poloidal magnetic fields. The aim to search for effect alternatives pursues double benefits. First, the nonhelical source of dynamo generations avoids the above mentioned catastrophic quenching problem. This issue is less important currently. Secondly, and it was already mentioned earlier by Köhler (1973) as well as Steenbeck and Krause (1969), the mixing length estimate of the effect for the solar convection zone parameters results in a very strong effect with a magnitude as strong as the convective velocity rms. Solar observations of the ratio between the typical strength of the toroidal and poloidal fields and the solar cycle period, favor an order of magnitude smaller effect. In addition, the turbulent generation sources in the Ω scenario help reduce the given constraints. We must stress that the GCD simulations of Schrinner (2011), Schrinner et al (2011), andWarnecke et al (2021) showed that the mean-field models need a full spectrum of turbulent effects to describe DNS.
In the case of a solar-like star, i.e., with solar-like stratification, differential rotation, and meridional circulation profiles, the turbulent sources of the poloidal magnetic field generation due to , shear-current, and cross-helicity effects are likely complimentary to the effect.
We thus arrive at the conclusion that the 2 Ω dynamo is, probably, the simplest scenario for the solar dynamo. Also, this scenario seems to fit well with observations of stellar activity of young solar-type stars.

Basic model
We discuss some results of the state-of-the-art mean-field dynamo model of a solar dynamo developed recently by Pipin and Kosovichev (2019) (hereafter PK19). The magnetic field evolution is governed by the mean-field induction equation: The expression for the components of E reads Here, describes the turbulent generation by the effect, represents turbulent pumping, and is the eddy magnetic diffusivity tensor. The effect tensor includes the small-scale magnetic helicity density contribution, i.e., the pseudoscalar ⟨a · b⟩, where is the dynamo parameter characterizing the magnitude of the kinetic effect, and K and M are the anisotropic versions of the kinetic and magnetic effects, as described in PK19. Note that, unlike Eq. (16), where the two contributions are pseudoscalars and have the same dimension, they are here tensorial where only K is a pseudotensor, but M is not, and they have here different dimensions.
The radial profiles of ( ) and ( ) depend on the mean density stratification, the profile of the convective velocity rms , and on the Coriolis number, where Ω 0 is the global angular velocity of the star and is the convective turnover time. The magnetic quenching function ( ) depends on the parameter = | |/( √︁ 4 rms ). In this model the magnetic helicity is governed by the global conservation law for the total magnetic helicity, ⟨A · B⟩ = ⟨a · b⟩ + · (see Hubbard and Brandenburg, 2012;Pipin et al, 2013): where we have used 2 ⟨j · b⟩ = ⟨a · b⟩ /Re M (Kleeorin and Rogachevskii, 1999). Also, we have introduced a diffusive flux of the small-scale magnetic helicity density, F f = − ∇ ⟨a · b⟩, and Re M is the magnetic Reynolds number, for which we employ Re M = 10 6 . Following results of Mitra et al (2010a), we put = 0.1 . Here, the turbulent fluxes of the magnetic helicity are approximated by the only term which is related to the diffusive flux. Besides the diffusive helicity flux, the other turbulent fluxes of the magnetic helicity can be important for the nonlinear dynamo regimes and the catastrophic quenching problem (Kleeorin et al, 2000;Vishniac and Cho, 2001;Chatterjee et al, 2011;Brandenburg and Subramanian, 2005a;Kleeorin and Rogachevskii, 2022;Gopalakrishnan and Subramanian, 2023). The relative importance of different kinds of magnetic helicity fluxes for the dynamo should be studied further.
The above ansatz for the helicity evolution differs from that given by Eq. (18); see also papers by Kleeorin and Ruzmaikin (1982) and Kleeorin and Rogachevskii (1999). Hubbard and Brandenburg (2012) had studied the magnetic helicity evolution for shearing dynamos. They found that employing Eq. (18) in the dynamo problem can result in nonphysical fluxes of magnetic helicity over spatial scales. For the ansatz given by Eq. (18), the nonlinear dynamo models can show sharp magnetic structures inside the dynamo domain. Such structures are connected with concentrations of the magnetic helicity; see, e.g., Chatterjee et al (2011) and Brandenburg and Chatterjee (2018). Even a strong diffusive helicity flux does not seem to correct for those irrelevant features from the numerical solution. The technical point is that the helicity fluxes, which are omitted in Eq. (18), should be consistent with the turbulent effects involved in the mean electromotive force, e.g., the rotationally induced anisotropy of the effect, the magnetic eddy diffusivity, etc. Such calculations are currently absent. Also, we have to take into account the modulation of the magnetic helicity density by the magnetic activity. On the other hand, with the magnetic helicity evolution equation Eq. (38), Pipin et al (2013) found that magnetic helicity density follows the large-scale dynamo wave. This alleviates the catastrophic quenching of the effect. They showed that if we write the Eq. (38) in the form of Eq. (18), we get an additional helicity flux due to the global dynamo, Rewriting Eq. (38) in the form of Eq. (18) we get where the ellipsis refers to additional helicity transport terms due to the large-scale flow. The term E × A consists of the counterparts of the sources of magnetic helicity, which are represented by −2E · B, and the fluxes which result from pumping of the large-scale magnetic fields. The sources of magnetic helicity in the term −2E · are partly compensated in Eq. (39) by the counterparts in E × . This results in the spatially homogeneous quenching of the large-scale magnetic generation and alleviation of the catastrophic quenching problem. The effect of E × was not unambiguously confirmed in DNS because of limited numerical resolution; see Del Sordo et al (2013) and Brandenburg (2018b). The turbulent pumping is expressed by the antisymmetric tensor . The tuning of for the solar-type mean-field dynamo model was discussed by Pipin (2018). We define it as follows, where ( ) = ∇ log , MLT = 1.9 is the mixing-length theory parameter, is the adiabatic law constant. In Eq. (40), the first term takes into account the mean drift of large-scale field due the gradient of the mean density, and the second one does the same for the mean-field magnetic buoyancy effect. The function H ( ) takes into account the effect of the magnetic tensions. It is H ( ) ∼ 2 for small and it saturates as −2 for ≫ 1; see P22.
We employ an anisotropic diffusion tensor following the formulation of Pipin (2008) (hereafter, P08): where the functions ( , ) 1,2 (Ω * ) are determined in P08. Analytical calculations of E in the above cited paper include the effects of a small scale dynamo. In the above expressions for E, we assume equipartition between kinetic energy of the turbulence and magnetic fluctuations which stem from the small-scale dynamo. It was found that for the case of slow rotation (Co ≪ 1), the part of E that depends on the gradients of consists of an isotropic eddy diffusivity and Rädler's × effect due to the small-scale dynamo (see also Rädler et al, 2003). In the case of rapid rotation, the fluctuating magnetic fields from the small-scale dynamo contribute both to isotropic and anisotropic parts of the diffusivity. The effect appears already in the terms of order Ω 2 in the global rotation rate (Rädler et al, 2003). In particular, the part of EMF which corresponds to Eq(42 can be written as It is noteworthy that the full expression for E obtained in P08 is complicated and includes other contributions due to the effects of global rotation , mean shear, mean current density , and the magnetic deformation tensor (∇ ). We skip them in the application to the solar dynamo model. The analytical results about the relations of the specific effects of E and the global rotation rate show qualitative agreement with the DNS of Käpylä et al (2009a) and Brandenburg et al (2012). Yet, a more detailed comparison of the analytical results and the GCD simulations is needed; for further discussions, see Sect. 7. We assume that the large-scale flow is axisymmetric. It is decomposed into the sum of meridional circulation and differential rotation, where is the radial coordinate, is the polar angle,ˆis is the unit vector in the azimuthal direction, and Ω ( , ) is the angular velocity profile. The angular momentum conservation and the equation for the azimuthal component of largescale vorticity, = (∇ × U m ) , determine the distributions of differential rotation and meridional circulation: whereT is the turbulent stress tensor: see the detailed description in Pipin and Kosovichev (2018) and PK19. Also, is the mean density, is the mean entropy; / = cos / −sin / · / is the gradient along the axis of rotation. The mean heat transport equation determines the mean entropy variations from the reference state due to the generation and dissipation of the large-scale magnetic field and large-scale flows: where is the mean temperature, F is the radiative heat flux, F is the anisotropic convective flux; see PK19. The last two terms in Eq. (46) take into account the convective energy gain and loss caused by the generation and dissipation of large-scale magnetic fields and large-scale flows. The reference profiles of mean thermodynamic parameters, such as entropy, density, and temperature are determined from the stellar interior model MESA (Paxton et al, 2015). The radial profile of the typical convective turnover time, , is determined from the MESA code, as well. We assume that does not depend on the magnetic field and global flows. The convective rms velocity is determined from the mixing-length approximation, where ℓ = MLT is the mixing length, MLT = 1.9 is the mixing length parameter, and is the pressure height scale. Equation (47) determines the reference profiles for the eddy heat conductivity , eddy viscosity , and eddy magnetic diffusivity , It should be noted that stellar convection might well have convection zones with slightly subadiabatic stratification in some layers. In those cases, the enthalpy flux can no longer be transported entirely by the mean entropy gradient, but there can be an extra term that is nowadays called the Deardorff term; see Deardorff (1972). Such convection can be driven through the rapid cooling in the surface layers and is therefore sometimes referred to as entropy rain Brandenburg (2016). It is useful to stress that the Deardorff term is distinct from the usual overshoot, because there the enthalpy flux points downward, while entropy rain still produces an outward enthalpy flux. It is instead more similar to semiconvection.
Boundary conditions. At the bottom of the tachocline, i = 0.68 , we assume solid body rotation and perfect conductor boundary conditions. Following to the MESA solar interior model, we put the bottom of the convection zone to b = 0.728 . At this boundary we fix the total heat flux, conv r + rad r = ★ ( )/4 2 . We introduce the decrease by exp (−100 / ) for all turbulent coefficients (except the eddy viscosity and eddy diffusivity), where is the distance from the bottom of the convection zone. The decrease of the eddy viscosity and eddy diffusivity is at most one order of magnitude for numerical stability. The meridional circulation is restricted to the convection zone. Therefore, we put the azimuthal component of the large-scale vorticity to zero, i.e., we set = 0 at b . At the top, t = 0.99 , we employ a stress free and black body radiating boundary. Following ideas of Moss and Brandenburg (1992), we formulate the top boundary condition in a form that allows penetration of the toroidal magnetic field to the surface: Free parameters. The model employs a number of free parameters, including , the turbulent Prandtl numbers Pr T and Pr M,T , , esq , and the global rotation rate Ω 0 . For the solar case we use a period of rotation of the solar tachocline determined from helioseismology, Ω 0 /2 = 434 nHz (Kosovichev et al, 1997). The best agreement of the angular velocity profile with helioseismology results is found for Pr = 3/4. Also, the dynamo model reproduces the solar magnetic cycle period, ∼ 20 years, if Pm = 10. Results of Pipin and Kosovichev (2011) showed that the parameters and esq affect the drift of the equatorial drift of the toroidal magnetic field field in the subsurface shear layer and magnitude of the surface toroidal magnetic field. Solar observations show the magnitude of the surface toroidal field to be about 1-2 G (Vidotto et al, 2018). To reproduce, it we use = 0.99 and esq = 50G. In what follows, we present the results of the solar dynamo model for a slightly supercritical parameter (10% above the threshold). Further details of the dynamo model can be found in PK19. Figure 8 shows profiles of the basic turbulent effects and large-scale flow distributions for the nonmagnetic case. The amplitude of the meridional circulation on the surface is about 13 m s −1 . In the lower part of the convection zone, the equatorward flow is about 1 m s −1 . The angular velocity profile is in agreement with helioseismology data.
Interestingly, the stagnation point of the meridional circulation is near the lower boundary of the subsurface shear layer, i.e., at = 0.9 . This is in agreement with observations of Hathaway (2012) and the helioseismic inversions of Stejko et al (2021). The structure of meridional circulation and turbulent pumping promotes an effective equatorward drift of the toroidal magnetic field below the subsurface shear layer; see Fig. 8(d).

Parker-Yoshimura dynamo waves and extended cycle
The dynamo shown in Fig. 9 demonstrates the numerical solution of the dynamo system including Eqs. (34) and (44)-(46). The time latitude diagrams of the surface radial magnetic field and the toroidal magnetic field in the upper part of the convection zone show agreement with observations of the evolution of the large-scale magnetic field of the Sun (Hathaway, 2015;Vidotto et al, 2018, see also the review of Righmire in this volume). The dynamo waves propagate to the surface-and equatorward. The radial direction of propagation follows the Parker-Yoshimura rule because of the positive sign of the effect in the main part of the convection zone and a positive latitudinal shear. It is noteworthy that at high latitudes, the model shows another dynamo wave family which propagates poleward along the convection zone boundary. This family follows the Parker-Yoshimura rule as well. Further we will see that the latitudinal shear plays a dominant role in this dynamo model and perhaps in the solar dynamo as well (see also Cameron and Schüssler, 2015). The latitudinal drift of the toroidal magnetic field in this model results from turbulent pumping and meridional circulation; see Fig8(d). The GCD simulations of Warnecke, 2018;Warnecke et al, 2021 show the crucial role of turbulent pumping in the solar type dynamo model, as well. The extended mode of the dynamo cycle is another feature of their model. The toroidal magnetic dynamo wave starts at the bottom of the convection zone at around 50 • latitude (see the marks in Fig. 9). It disappears near the solar equator after a full dynamo cycle. On the surface, the extended mode of the solar cycle is seen in the radial magnetic field evolution, in the torsional oscillations of zonal flow, and in the variations of the meridional circulation as well (Getling et al, 2021). The origin of the extended mode of the dynamo cycle is due to the distributed character of the large-scale dynamo and the interaction of the global dynamo modes, where the low order dynamo modes, e.g., dipole and octupole modes, are mainly generated in the deep part of the convection zone. The high order modes are predominantly generated in the near surface level. The phase difference between the models results in a dynamo mode of the extended length (Stenflo, 1992;Obridko et al, 2021).

Torsional oscillations
Solar zonal variations of the angular velocity ("torsional oscillations") were discovered by Howard and Labonte (1980). Since that time it was found that torsional oscillations represent a complicated wave-like pattern which consists of alternating zones of accelerated and decelerated plasma flows (Snodgrass and Howard, 1985;Altrock et al, 2008;Howe et al, 2011). Ulrich (2001) found two oscillatory modes of these variations with periods of 11 and 22 years. Torsional oscillations were linked to ephemeral active regions that emerge at high latitudes during the declining phase of solar cycles, but represent the magnetic field of the subsequent cycle (Wilson et al, 1988). It is interesting that in their original paper, Howard and Labonte (1980) conjectured that the solar torsional oscillation can shear magnetic fields and induce the dynamo cycle. This idea was further elaborated upon in a number of papers. However, the idea looks unreasonable because it conflicts with Cowling's theorem. Also, the magnitude of the torsional oscillations of 3-6 m s −1 is too small in comparison with the magnitude of the magnetic field generated by a dynamo. The first papers by Schuessler (1981) and Yoshimura (1981) suggested that the 11-year solar torsional oscillation can be explained by the mechanical effect of the Lorentz force. The double frequency of the zonal variation results from the 2 modulation of the large-scale flow due to the dynamo activity. On the basis of a flux-tube dynamo model, Schuessler (1981) used the simple estimated of a large-scale Lorentz force and found both the 11 and 22 year mode of the torsional oscillations. This result was further elaborated upon by Kleeorin and Ruzmaikin (1991). The further development of the mean-field theory of the solar differential rotation showed that, in addition to the large-scale Lorentz force, the dynamo induced 2 modulation of the turbulent angular momentum fluxes is also an essential source of the torsional oscillations (Rüdiger and Kitchatinov, 1990;Kleeorin et al, 1996;Küker et al, 1996;Rüdiger et al, 2012). Global convective dynamo simulations (e.g., Beaudoin et al, 2013;Käpylä et al, 2016a;Guerrero et al, 2016) confirmed this conclusion. The strength of the solar torsional oscillations is more than two orders of magnitude less than the differential rotation. It looks like the theory of the torsional oscillations can be constructed using perturbative approximations. Models of this type (see, e.g., Tobias 1996;Covas et al 2000;Bushby and Tobias 2007;Pipin 2015;Hazra and Choudhuri 2017) were inspired by results of Malkus and Proctor (1975). Yet, the constructed models are incomplete because they ignored the Taylor-Proudman balance, which is a key ingredient of solar differential rotation theory (see Kitchatinov 2013, also contribution of Hazra et al, this volume). Complete mean-field dynamo models, which take into accounts the Taylor-Proudman balance (hereafter TPB), were constructed by Brandenburg et al (1992), Rempel (2007), and PK19. Figure 9 shows variations of the zonal acceleration for our mean-field model in following the PK19 line of work. Similar to the results of helioseismology (Howe et al, 2011;Kosovichev and Pipin, 2019) and the results of Rempel (2007), snapshots of the model show that in the main part of the convection zone, the acceleration patterns are elongated along the rotation axis. This is caused by the Taylor-Proudman balance. Near the convection zone boundaries, these patterns deviate in the radial direction, which is in agreement with the above cited helioseismology results, as well. The given observation on the role of TPB shows the importance of the meridional circulation and the dynamo-induced heat transport perturbation (Spruit, 2003;Rempel, 2007) in the theory of torsional oscillations. This fact does not deny the importance of the large-scale Lorentz force and the magnetic modulation of the turbulent angular momentum transport. Results of Figures 9(b) and (c) show that the positive sign of the zonal acceleration propagates from high latitude bottom of the convection zone toward equator sticking to the equatorial edge of the dynamo wave. The torsional oscillation wave is accompanied by corresponding variations of the meridional circulation. These variations are induced by magnetic perturbations of the heat transport (see details in PK19). We emphasize that the given dynamo models also show overlapping magnetic cycles; see Fig. 9(b), similarly to what was originally proposed by Schuessler, 1981]. In this case, the 2 effect of the dynamo on the heat transport and the TPB results in about 4 to 5 meridional circulation cells along latitude. This tracks the zonal variations of angular velocity, which are caused by the mechanical action of the large-scale Lorentz force and magnetic quenching of the turbulent stresses, from polar regions to the equator. PK19 found that the induced zonal acceleration is ∼ (2-4) × 10 −8 m s −2 , which is in agreement with the observational results of Kosovichev and Pipin (2019). However, the individual forces in the angular momentum balance such that the large-scale Lorentz force, the variations of the angular momentum transport due to meridional circulation, the inertial forces, and others are by more than an order of magnitude stronger than their combined action and can reach a magnitude of ∼ 10 −6 m s −2 . Therefore, the resulting pattern of the torsional oscillations forms in nonlinear balance, which includes the forces driving the angular momentum transport, the TPB, and heat perturbations due to magnetic activity in the convection zone (see details in PK19).

Mean-field models based on the EMF obtained from DNS
Here we provide an example of how the mean-field theory is utilized as a tool for understanding what is going on (see Sect. 3.4). We discuss recent studies of meanfield dynamo models constructed based on the electromotive force (EMF) obtained from direct numerical simulation (DNS) of rotating stratified convection, especially focusing on "semi-global" models. The properties of solar and stellar convection, and the various methods for extracting the information of the EMF from DNS are also summarized.

Properties of solar and stellar convection
A quantitative physical description of solar and stellar dynamos, which should be the result of the nonlinear interaction of turbulent flows and magnetic fields, is a great challenge for us and constitutes a significant milestone on the long way to a full understanding of turbulence. Even with state-of-the-art supercomputers, it is impossible to numerically simulate solar and stellar convection and its interaction with the magnetic field and to observe/analyze numerical data in detail with realistic parameters. Therefore, to say with confidence that one has fully understood the solar and stellar dynamo problem, it should be necessary to find a universal law of magnetohydrodynamic (MHD) turbulence, build a reliable sub-grid scale (SGS) turbulence model, and then reproduce the magnetic activities of the Sun and stars quantitatively in an integrated framework by numerical models with incorporating the SGS model. This is because fluid quantities that may be verified in future observations should include the meridional distributions of fluid velocity, vorticity, kinetic helicity, and thus the turbulence model constructed on the basis of these profiles (e.g., Hanasoge et al, 2016). Only when the correctness of the turbulence model is observationally validated should our understanding of the solar and stellar dynamos as a consequence of the turbulent dynamo process be completed. In the near future, a very exciting time may come when we will be able to test and verify various turbulence models under extreme conditions inside the solar and stellar interiors. What physical characteristics should be taken into account when constructing a turbulence model of thermal convection in the Sun and stars? Let us summarize some essential features: 1. Extremely low dissipation: turbulent state with Re ≳ 10 12 , Re M ≳ 10 8 , and a large Pèclet number, Pe ∼ 10 6 − 10 9 (where Pe = Re · Pr). 2. Huge separation of dissipation scales: Pr ∼ 10 −4 -10 −7 , Pr M ∼ 10 4 3. Compressibility: high Mach number O (1) in the upper convection zone makes the convective motion compressible. 4. Anisotropy: spin of stars (i.e., Coriolis force in a rotating system) makes fluid motions anisotropic. 5. Inhomogeneity: density contrast of 10 6 between top and bottom CZs results in multi-scale properties of fluid motion. 6. Non-locality: Radiative energy loss at the CZ surface (open system), allowing the growth of cooling-driven downflow.
In view of these features, it can be seen that the characteristics of thermal convection operating inside the Sun and stars are quite different from those of isotropic turbulence. Those can be considered to some extent in DNS even with the current computing performance, as listed under 3-6, while the others, (items 1 and 2) are unreachable with current grid-based simulations. It should be emphasized, however, that higher resolution simulations using state-of-the-art supercomputers is a classical way forward in turbulence research, and the knowledge obtained from such studies in unexplored low-dissipation regimes will greatly expand the horizon of our understanding of turbulence (e.g., Kaneda et al, 2003;Hotta and Kusano, 2021). Moreover, if sufficient scale separation between the turbulent and mean fields is ensured and the inertial range of the turbulent cascade is captured appropriately, there is the possibility that the evolution of mean-field components, such as largescale flow and large-scale magnetic field, can be approximately reproduced even by simulations with enhanced dissipation compared to the actual solar and stellar values (e.g., Ossendrijver, 2003). It should be remembered, however, that in spite of the rapid increase in computing power, some rather basic questions about the solar dynamo still remain, for example the equatorward migration of the sunspot belts and the formation of sunspots themselves.

Semi-global simulation of rotating stratified convection
On our way toward a reliable SGS turbulence model for solar and stellar interiors, numerical models of convection and its dynamo should be studied, while keeping the characteristic features of solar and stellar convection, as listed under items 3-6 above, in mind. It should be noted that the underlying necessity for numerical modeling is an important component of earlier studies that applied mixing-length type concepts to the dynamo theory, which never successfully explained the magnetic activities of the Sun and stars (e.g., Brandenburg and Tuominen, 1988). In recent years, significant progress has been made in GCD simulations (e.g., Browning et al, 2006;Ghizaru et al, 2010;Käpylä et al, 2012;Masada et al, 2013;Fan and Fang, 2014;Augustson et al, 2015;Hotta et al, 2016;Warnecke, 2018), there is also a growing effort to extract the information of turbulent transport processes from so-called "semi-global" (or local model) MHD convection simulations with the aim of quantifying the dynamo effect of rotational stratified convection (e.g., Brandenburg et al, 1990Brandenburg et al, , 1996Nordlund et al, 1992;Brummell et al, 1998Brummell et al, , 2002Ossendrijver et al, 2001;Käpylä et al, 2006aKäpylä et al, , 2009bSano, 2014b,a, 2016;Bushby et al, 2018;Masada and Sano, 2022). A typical numerical setup of the semi-global model is shown in Fig. 10 schematically. In this setting, the gas is gravitationally stratified in the vertical direction, while periodicity is assumed in the horizontal directions. The governing equations (mostly compressible MHD equations) are solved in a rotating Cartesian frame, and the rotation axis is usually set to be parallel or anti-parallel to the gravity vector. Several studies have simulated the model with the tilt of the rotation axis with respect to the gravity vector, and the latitudinal dependence of the convection has been investigated (e.g., Ossendrijver et al, 2001;Käpylä et al, 2004Käpylä et al, , 2006a.

Extraction of information of dynamo effects
In the semi-global studies, four-types of approaches have been used typically to extract the information of dynamo effects veiled in the convective motion. The starting point of all the four methods is common, the decomposition of the flow field ( ) and magnetic field ( ) into a spatially large-scale, slowly-varying meancomponent, and a small-scale, rapidly varying fluctuating component, as introduced in § 1, i.e., = + and = + , where the lower-case represent the fluctuating component and the overbars denote the mean component. In the case of a semiglobal model, a temporal and horizontal average is often used for deriving the mean component. Then, the equation of mean-field electrodynamics can be derived where E = ′ × ′ is the mean electromotive force (EMF) due to the fluctuation of the flow and the magnetic field. The mean EMF can be described as a power series about the large-scale magnetic component and its derivatives as where represents (tonsorial form of) the -effect, is the turbulent pumping, and denotes the turbulent diffusion.
To obtain the information of the dynamo coefficients, such as , , and , from the MHD convection simulation, there are the following four methods (see general discussion in Sect. 3.3): (i) Methods based on results of analytical theories, e.g., SOCA (or first-order smoothing approximation, FOSA) expressions (ii) Imposed-field method (iii) Test-field method (iv) The multi-dimensional regression method based on the dynamo generated magnetic field Method (i) involves the the estimation of dynamo coefficients based on results of analytical theory. It also exploit the mixing-length approximation in final results. There, the distributions of, for example, the fluctuating components of the convection velocity ( ), vorticity ( = ∇ × ), and the resulting kinetic helicity (H = · ), are directly extracted from the simulation results and used to reconstruct the turbulent and via their analytic forms, derived under SOCA, such as Eq. (16) and = ( /3) 2 , where is the correlation time of the turbulence and is often replaced by the convective turnover time, . Note that anisotropy effects are often neglected in the expressions above, but see Brandenburg and Subramanian (2007), who included them.
Method (ii) is mainly used in the analysis of the numerical results without selfsustained magnetic field. In this method, a uniform external magnetic field is imposed as the mean component to the computational domain, artificially. The effect and turbulent pumping are then inferred from E = × by directly calculating from simulation data Furthermore, one might be tempted to compute t = E · /( 0 2 ), but this would assume that · is vanishing, which is generally not the case for -effect dynamos; see  for details. Method (iii) utilizes a so-called test field, as introduced by Schrinner et al (2005,2007) for the spherical case and Brandenburg et al (2008) for the Cartesian case, allowing for scale dependence. In this method, the evolution equation of ′ T , the fluctuating component of the test field T , which are passive to the velocity field taken from the simulation, is solved additionally to the basic (MHD) equations. From the linear evolution of the test-field, the mean EMF is evaluated and then the full set of turbulent transport coefficients can be obtained. For example, in the case without the large-scale flow, the test-field equation is, for T , with a chosen test field T while taking from the MHD simulation. The original method does not work in cases when magnetic fluctuations are driven through an artificially added electromotive force. In such cases, one needs to use a more general nonlinear method explored by Rheinhardt and Brandenburg (2010) and Käpylä et al (2022). Theses magnetically driven systems are also examples where the imposedfield method gives reliable results in two dimensions with volume averaged mean fields, where no turbulent diffusion can act. Method (iv) can be used only in the analysis of the numerical results with selfsustained magnetic fields. Since the fluctuating and mean components are all known quantities in such simulations, the mean emf, E = × , and the mean magnetic component, , can be directly calculated from the simulation data. Then, the mean profiles of dynamo coefficients are inferred based on a fitting procedure via the relationship, Given E and which are calculated from simulation data, and then find and such that the residual of Eq. (56) is minimized. In the equation above, the contributions from the derivatives of the mean magnetic component to the mean emf are neglected (see, e.g., Racine et al, 2011;Simard et al, 2013Simard et al, , 2016Shimada et al, 2022, for the fitting based analysis of the dynamo coefficient with including the contribution from the first-order derivative of the mean magnetic component). As we have mentioned in Sect. 3.3, the results of Warnecke et al (2018) show that the TFM is more accurate compared to the regression method in obeying Eq. (56); see the detailed comparison in above cited paper.
In all cases, however, the first (and often higher) derivative terms are of the same order as the first term and can therefore not be neglected. This was already done in the work of Brandenburg and Sokoloff (2002), who typically found small diffusion coefficients in the cross-stream direction. This, however, turned out to be a shortcoming of the method and has not been borne out by more advanced measurements (Karak et al, 2014).

Transport coefficients from semi-global turbulence simulations
Here, we briefly review the results of previous semi-global simulations, with a particular focus on the studies that have been dedicated for extracting information about dynamo coefficients. Brandenburg et al (1990), hereafter B90, performed turbulent 3-D magnetoconvection simulations under the influence of the rotation for the semi-global model whose depth is equivalent to about one pressure scale height. They found that, due to the effect of the rotation, a systematic separation of positive and negative values of the kinetic helicity was developed in the vertical direction of the CZ, i.e., in the upper CZ, negative (positive) helicity in the northern (southern) hemisphere, while positive (negative) helicity in the northern (southern) hemisphere. Using the imposed field method, they evaluated the magnitude of the turbulent -effect with anisotropic properties as /( H ) ∼ O (0.1) and /( H ) ∼ O (0.01), where H = · and E = H H + V V . It is interesting to note that these values are about one to two orders of magnitude smaller than ∼ Ω , which is the estimation based on the mixing-length theory. Additionally, it was also suggested that the magnetic helicity showed a similar depth variation, but the sign was opposite to that of the kinetic helicity.
While H had the expected sign (opposite to that of the kinetic helicity), V was found to have the 'wrong' sign (same as that of the kinetic helicity). Such a result was subsequently also obtained by Ferriere (1993). The theoretical possibilities for such effects should be studied further. For example, Rüdiger and Pipin (2000) found that large-scale shear can affect both the sign of the effect and kinetic helicity in magnetically driven compressible turbulence in such a way that they have the same sign, e.g., for Keplerian accretion disks. These ideas were also applied to understanding the finding of a negative effect in stratified accretion disk simulations (Brandenburg, 1998). Ossendrijver et al (2001) also performed the semi-global simulation with a similar model as B90. They showed that, even in the regime where the condition justifying the FOSA (or SOCA) is not satisfied, i.e., in the situation where St = rms / ≳ 1 and > 1, the kinetic helicity was clearly separated into positive and negative values at the lower and upper CZs when taking temporal average of the convective motion over sufficiently long time. Using the imposed field method, they also measured the magnitude of the turbulent -effect and obtained similar values to B90 in terms of and . The rotational dependence of the -effect was also investigated in this work for the first time. They showed that, in the larger Co regime, the underwent a rotational quenching, while the was saturated, where Co is the Coriolis number [see Eq. (37)]. The turnover time was defined, in this work, as = / rms . While the depth-dependence or rotational dependence of the , which was obtained from the simulation, agreed, to some extent, with a theoretical model based on the mixinglength theory (Rüdiger & Kitchatinov 1993), their amplitudes were one to two orders of magnitude smaller than those predicted from the theoretical model. Noteworthy, the critical threshold of the effect parameter in mean-field dynamo models (see Sect. 6.1) is about same magnitude less than the mixing-length models of the solar convection zone predicts; see Sect. 6.
In Käpylä et al (2004Käpylä et al ( , 2006a, additionally to the rotational dependence, the latitudinal dependence of the turbulent -effect was studied in the semi-global convection simulations with varying the inclination of the rotation axis with respect to the gravity vector. With the imposed field method, they found that, for slow and moderate rotation with Co < 4, the latitudinal dependence of the followed cos profile with a peak at the pole (see also, Egorov et al, 2004), while, in the rapid rotation regime with Co ≈ 10, it rather peaked much closer to the equator at ≃ 30 • . Additionally, the vertical profile of the directly evaluated from simulation was found to be qualitatively consistent with analytic expression derived under the FOSA even when changing the latitude. A practical application of these results was the development of a kinematic mean-field solar dynamo model in Käpylä et al (2006b). In it, the rotation profile deduced from the helioseismic observation and the meridional profiles of the -effect and turbulent pumping obtained with the semi-global simulation of Käpylä et al (2006a) are integrated into the framework of the -Ω dynamo, and then the solar dynamo mean-field model was constructed. It is interesting that their kinematic dynamo model correctly reproduced many of the general features of the solar magnetic activity, for example realistic migration patterns and correct phase relation.
The existence of large-scale dynamo, i.e., self-excitation of the mean magnetic component, in rigidly-rotating convection was demonstrated for the first time in the semi-global simulation by Käpylä et al (2009b). By changing the angular velocity, they showed that the large-scale dynamo could be excited only when the rotation is rapid enough, i.e., Co ≳ 60, with Eq. (37) as the definition of Co which is same as that used in Ossendrijver et al (2001) and Käpylä et al (2006a); see, e.g., Tobias et al (2008) and Cattaneo and Hughes (2006), and Favier and Bushby (2013), for unsuccessful large-scale dynamo in rigidly-rotating convection probably due to slow rotation, and/or short integration time. From the measurements of the turbulenteffect and the turbulent diffusivity by the TFM, they also suggested that while the magnitude of the -effect stayed approximately constant as a function of rotation, the turbulent diffusivity decreased monotonically with increasing the angular velocity, resulting in the excitation of the large-scale dynamo in the higher Co. The reliability of the dynamo coefficients extracted with the test-field method from the simulation was validated with the one-dimensional mean-field dynamo model in which the testfield results for and were used as input parameters by studying the excitation of the large-scale magnetic field at the linear stage. Note that the oscillatory properties of the large-scale dynamo in rigidly-rotating convection and its possible relationship with 2 dynamo mode with inhomogeneous profile were also found in Käpylä et al (2013); see, e.g., Baryshnikova and Shukurov (1987) and Mitra et al (2010b) for the oscillatory 2 dynamo.

Weakly-stratified Model
Below we review recent mean-field dynamo models linked with semi-global MHD convection simulations, where the large-scale dynamo is successfully operated; see Masada and Sano (2014a,b, 2016, 2022 for a series of numerical studies.
While Käpylä et al (2009bKäpylä et al ( , 2013 were the first to demonstrate that rigidly-rotating convection can excite the large-scale dynamo as reviewed above, their simulation model was a so-called "three-layer polytrope" consisting of top and bottom stablystratified layers and the CZ in between them. Therefore, it was suspected for a while that the essential factor for the successful large-scale dynamo observed there might be the presence of the stably-stratified layer assumed in their model rather than the rapid rotation (e.g., Favier and Bushby, 2013). To pin down the key requirement for the large-scale dynamo, the impact of the stably-stratified layers on the large-scale dynamo was studied in Masada and Sano (2014b), hereafter MS14a, in which twotypes of semi-global models with and without stably-stratified layers are compared with the same control parameters and the same grid spacing. It was found in this study that a large-scale dynamo was successfully operated even in the model without the stably-stratified layer, and confirmed that the key requirement for it should be a rapid rotation if we evolved the simulation for a sufficiently long time than the ohmic diffusion time. Note that a relatively weak density stratification (the density contrast between the top and bottom CZs is about 10) was assumed in the simulation model employed in this study as well as Käpylä et al (2009bKäpylä et al ( , 2013. With these results, Masada and Sano (2014a), hereafter MS14b, explored the mechanism of the large-scale dynamo operated in the rigidly-rotating stratified convection by linking the mean-field (MF) dynamo model with the DNS. In this study, the FOSA based approach was adopted in the MF modeling. The mean vertical profiles of the kinetic helicity and root-mean square velocity were directly extracted from the simulation data and then the vertical profiles of the turbulent , turbulent pumping ( ) and turbulent diffusivity ( ) were reconstructed according to the analytic expressions of in anisotropic forms of dynamo coefficients under the FOSA (e.g., Käpylä et al, 2006a). Although recent numerical studies indicate that the small-scale current helicity, i.e., · , is important for the -effect when the magnetic field is dynamically force (e.g., Ossendrijver et al. 2002). The coefficients α, γ , and η represent the α-effect, turbulent pumping, and turbulent magnetic diffusivity, respectively. All the terms related to u h and Bz h can be ignored in considering the symmetry of the system. All the variables, except for η0, depend on the time (t) and depth (z). The MF dynamo described by Equation (2) falls into the α 2 -type category. The MF theory predicts that the α 2 mode can generate a large-scale magnetic field with an oscillatory nature (e.g., Baryshnikova & Shukurov 1987;Rädler & Bräuer 1987;Brandenburg et al. 2009). A key ingredient for the oscillatory mode is the nonuniformity of the α-effect, which can arise naturally as an outcome of rotating stratified convection in the stellar interior. Using the rigidly rotating system studied here, the α 2 dynamo wave was excited, which propagates only in the depth direction. However, as shown by Käpylä et al. (2013b), in the global system, it can travel also in the latitudinal direction due to the strong antisymmetry of the α-effect across the equator.
The dynamo-generated MF produces a Lorentz force that tends to "quench" the turbulent motions and control the nonlinear evolution and saturation of the system. Since there is no definitive model to describe the magnetic quenching effect (e.g., Rogachevskii & Kleeorin 2001;Blackman & Brandenburg 2002) as yet, we adopt the prototypical models, which are the dynamical α-quenching, algebraic γ -and η-quenching of the catastrophic-type; (see Brandenburg & Subramanian 2005, for the quenching), where ReM = ηk/η0. The dependence of the MF model on the quenching formula should be discussed in detail in a subsequent paper; however, at least the conclusions of this Letter remain independent from the choice of the quenching models. The characteristic wavenumber kc and the equipartition field strength Beq are given by kc(z) = 2π/Hd and Beq(z) = ρuz 2 h in our model, where Hd = −dz/d ln ρ h is the density scale height.
Here, the subscript "k" refers to the unquenched coefficient, which is calculated from DNS results of the saturated convective turbulence.
In the first-order smoothing approximation (FOSA), the unquenched coefficients αk, γk, and ηk in anisotropic forms are given by (e.g., Käpylä et al. 2006Käpylä et al. , 2009b where τc is the correlation time, Heff is the effective helicity, and urms is the rms velocity. The vertical profiles of Heff and u 2 rms in the reference DNS model are shown in Figure 1 important (Pouquet et al, 1976;Brandenburg and Subramanian, 2005b), its contribution was ignored in this study. As the correlation time , the convective turnover time defined by = ( )/ rms was chosen there ( is the density scale-height as a function of the depth). By solving one-dimensional MF 2 dynamo equation in which these profiles were used as input parameters, i.e., with the time-depth diagram for the mean (horizontal) magnetic component ( ℎ ) was obtained. In Fig. 11, we show ( , ) for the MF model (panels (b)) and its DNS counterpart (panels (a)). Note that, for ensuring the saturation of the magnetic field growth, the quenching effect was also taken into account. Since the DNS results were quantitatively reproduced by the MF 2 dynamo, MS14b concluded that the large-scale magnetic field organized in the rigidly-rotating turbulent convection was a consequence of the oscillatory 2 dynamo.
Reproducing the DNS results with mean-field models using coefficients from the original DNS is an important verification of the whole approach. This has been done on many occasions in the past; see, for example, the work by Gressel (2010) and Warnecke et al (2021).

A strongly stratified model
In MS14a,b, a weakly-stratified model, in which the density contrast between top and bottom CZs is about 10, was adopted. However, the actual Sun has very strong stratification with a density contrast of ∼ 10 6 between top and bottom CZs, resulting in a large segregation of time scales from minutes to months. Aiming at the   Fig. 12 (a) 3-D view of the strongly stratified convection (for the progenitor run without rotation). The black (gray) tone denotes downflows (upflows). (b) Time-depth diagrams for and . The normalization is the equipartition field strength, eq ≡ ( 2 ) 1/2 . In MS16, a one-layer polytrope with a super-adiabaticity of ≡ ∇ − ∇ ad = 1.6 × 10 −3 was used; see MS16 and MS22 for details. application to solar and stellar interiors, Masada and Sano (2016), hereafter MS16, performed a convective dynamo simulation in a strongly stratified atmosphere with a density contrast of 700 in a semi-global setup. Due to the strong solar-and stellar-like density stratification, multi-scale convection with a strong up-down asymmetry, i.e., slower and broader upflow regions surrounded by a network of faster and narrower downflow lanes, was developed in this simulation, as shown in Fig. 12(a). Even in such a situation, a large-scale dynamo was found to operate. As shown in Fig. 12(b), the mean magnetic field components observed there showed a time-depth evolution similar to that seen in the weakly-stratified model (MS14a,b), suggesting that an oscillatory 2 dynamo is responsible for it. It was intriguing that, in addition to the mean horizontal component, the large-scale structures of the vertical magnetic field were spontaneously organized at the CZ surface in the case of the strongly stratified atmosphere, as shown in Fig. 13.
A possible physical origin of such surface magnetic structure formation is the negative magnetic pressure instability (NEMPI; see § 8 for details). NEMPI is a mean-field process in the momentum equation, where the Reynolds and Maxwell stresses attain a component proportional to the square of the mean magnetic field, which acts effectively like a negative pressure by suppressing the turbulent pressure. Since its growth rate becomes larger for stronger density stratification (e.g., Jabbari et al, 2014), one can imagine that it may play an important role in organizing sunspot-like large-scale magnetic field structures in the upper part of the solar CZ. Although its presence has been confirmed numerically in forced MHD turbulence (e.g., Brandenburg et al, 2011;Warnecke et al, 2013), it does not play a significant role in organizing the surface magnetic structure seen in MS16 because of their relatively rapid rotation; Ro = 0.02 was assumed there, while, according to Losada et al (2012), Ro ≳ 5 is required to excite the NEMPI.
The large-scale structure of the vertical magnetic field observed in MS16 is similar to that observed in the large-scale dynamo by forced turbulence in a strongly stratified atmosphere (Mitra et al, 2014;Jabbari et al, 2016). This implies that there would and the local convective turnover time defined by ( ) ( ) t = á ñ r z H z u z cv,l 2 h 1 2 (dashed) are compared as a function of the depth, where r º á ñ r H dz d ln h . In the upper CZ, the condition t t cv,l pk,max is always satisfied. Since, in such a situation, the small-scale convective motion violently disturbs the coherency of the magnetic flux, we would have to say that the Parker instability would not be responsible for the largescale structure formation observed in our simulation.
Next, the large-scale flow and its association with the surface magnetic structure are analyzed. For casting light on the largescale pattern, the small-scale structures with k/k c  8 are eliminated by applying Fourier filtering (e.g., Warnecke et al. 2016;Jabbari et al. 2016). A series of snapshots where á ñ B z k8 and á ñ u z k8 on the horizontal plane at z/d cz = 0.04, are shown in Figure 5(a), where · á ñ k8 denotes Fourier filtering. The overplotted arrows are the velocity vectors composed of á ñ u x k8 and á ñ u y k8 . Additionally, 2D spectra of B h 2 , B z 2 , rv h 2 , and rv z 2 in the upper CZ are also show each depth is spatially av from 0.0 to 0.25 and is tem the corresponding reference It is significant that, in bipolar "band-like" structu the horizontal magnetic flu Although the faster horizon be found in/around the r dynamo-saturation (left col not necessarily associated dynamo-saturated stage ( addition, we can find from t in the large-scale magnetic of the large-scale flow in large-scale flows are not th large-scale magnetic structu Appendix A, hByih also shows a similar cyclic behavior with 348 hBxih, but with a phase difference of ⇡/2, suggesting that an 349 "↵2-type" dynamo is excited in our DNS models (see, e.g.,

454
where He↵ is the effective helicity, and ⌧cor is the corre- In the previous section, we verified that the spin rate of 395 the system is essential for exciting the large-scale dynamo.

396
The critical Rossby number that separates the success or fail-   be an as-yet-unknown mechanism for the self-organization of large-scale magnetic structures, which would be inherent in a strongly stratified atmosphere.
In Masada and Sano (2022), hereafter MS22, with varying angular velocity as a control parameter, the rotational dependence of the large-scale dynamo was explored in a series of DNSs of rigidly-rotating convection. They linked its cause through MF dynamo models with DNSs where a strongly stratified polytrope was adopted as a model of the convective atmosphere, as in MS16.
In Fig. 14(a), DNS results are shown where a time-depth diagram of is depicted for models with different values of Co. While in the slowly rotating model with low Co, the large-scale magnetic component fails to grow, it is found to be spontaneously organized in the rapidly-rotating models with high Co. It was found from DNS that the large-scale dynamo was excited when Co ≳ Co crit , where Co crit is the critical Coriolis number in the range 25 ≲ Co crit ≲ 80, with Eq. (37) as the definition of the Coriolis number. It is remarkable that Co crit , which determines the success or failure of the large-scale dynamo, was almost the same -independent of the density contrast (see Käpylä et al, 2009b) or the geometry of the simulation model (see, e.g., Käpylä et al, 2012;Warnecke, 2018, for Co crit in the global simulations); see MS22 for the quantitative comparison between models.
To explore the underlying mechanism of the rotational dependence of the largescale dynamo, the influence of the rotation on the turbulent transport coefficients was also studied in MS22 with the FOSA-based approach similar to that adopted in MS14b. In Fig. 14(c), the vertical profiles of the turbulent effect and turbulent diffusivity reconstructed with the analytic expressions of Eq. (57) were shown. It can be found that, as the spin rate increases, the turbulent diffusivity decreases significantly, but the profile of the effect remains almost unchanged. This result suggested that the rotational dependence of the large-scale dynamo observed in MS22 may be primarily due to a change in the magnitude of the turbulent diffusion. In fact, this insight was confirmed by the evidence that the MF dynamo model with incorporating the dynamo coefficients shown in Fig. 14 reproduced quantitatively the result of the DNS; see Fig. 14(b) for the time-depth diagram of obtained in the MF models with using different dynamo coefficients extracted from the corresponding DNSs.
MS22 concluded that the independence (dependence) of the turbulent effect (turbulent diffusivity) on the rotation was the essence of the rotational dependence of the dynamo. This is not only the same as the conclusion obtained by Käpylä et al (2009b) from weakly-stratified convective dynamo simulations using the TFM, but also the same as that obtained by Shimada et al (2022) from global solar dynamo simulation with using the "self-sustained field method". Although we don't know whether the independence (dependence) of the -effect (turbulent diffusion) on the rotation, seen in these studies, is universal or not, it may give an important suggestion not only on the turbulence modeling but on the solar dynamo modeling.

Solar-stellar connections and questions beyond standard mean-field theory
The title of our review is "Turbulent processes and mean-field dynamo" Obviously, there are mean-field effects in turbulence that are not just of dynamo type, and the mean-field dynamo is not just related to the solar dynamo problem, but its relevance goes much beyond. Here, we highlight, just some effects, but we refer to the reviews of Käpylä et al (2023) for additional examples.

Origin of sunspots and active regions
An important goal in solar dynamo theory is to compute synthetic butterfly diagrams. The question then emerges from which depth to take the mean toroidal field, for example. The usual argument here is to invoke Parker's theory of sunspot formation and to postulate that the field at some depth translates directly to one at the surface. This is critical because the final result depends on the assumed depth. It is possible that sunspots are not deeply rooted, but are actually a surface phenomenon. No successful and self-consistent model of shallow formation of active regions or sunspots exists as yet. Noteworthy in this context is NEMPI, which is a mean-field theory of the Reynolds and Maxwell stresses. This theory is extremely successful in that its results agree remarkably well with direct numerical simulations (DNS).2 The problem is only that the effect is not strong enough to make real sunspots or active regions. Because of this remarkable agreement between theory and simulations, and because it is an important mean-field process, we shall discuss here a bit more detail.
The essence of the effect is the contribution of the turbulent hydromagnetic pressure to the horizontal force balance. The turbulent pressure is a small-scale effect, but it reacts to the large-scale magnetic field. As the magnetic field increases, it suppresses the turbulence locally, disturbing therefore the horizontal force balance. Although this large-scale magnetic field itself contributes with its own magnetic pressure to the horizontal force balance, the effect from the suppression of the turbulence is often stronger, so the net effect is a negative one. This is why the mean-field effect from a large-scale magnetic field is a negative effective magnetic pressure. This idea goes back to early work of Kleeorin et al (1989Kleeorin et al ( , 1996, who developed the mean-field theory for this effect. In the beginning, it was not clear what kind of numerical experiments one could try to test the a negative effective magnetic pressure effect. The first mean-field simulations were done with a uniform horizontal magnetic field . This led to the development of magnetic flux concentrations near the surface, but those began to sink downward as time went on. A similar effect was soon also seen in DNS . The sinking of such structures was explained by the negative effective magnetic pressure: a positive magnetic pressure would lead to the rise of structures (Parker, 1967) while a negative one leads to a sinking. The sinking of magnetic structures had the side effect that the structures disappeared from the surface and became less prominent.
Subsequent experiments with a vertical field had a more dramatic effect on the general appearance of structures. Because the ambient field was vertical, the downflow had little effect on the magnetic flux concentrations themselves (Brandenburg Fig. 15 Cuts of / eq ( ) in the plane at the top boundary ( / = ) and the plane through the middle of the spot at = 0. In the cut, we also show magnetic field lines and flow vectors obtained by numerically averaging in azimuth around the spot axis. Adapted from . Fig. 15 shows the spontaneous development of a magnetic spot.  found that NEMPI saturates slightly below equipartition value.
Most of the numerical experiments where done with forced turbulence, where one had explicit control over the degree of scale separation. For the results presented here Re M ∼ 19 and Pr M = 0.5. The value of Re M = rms / f is very small because it is based on the forcing wavenumber f , which is chosen to be large. Furthermore, Pr M is chosen to be less than unity, because NEMPI is not expected to work for Pr M ≥ 14. This is rather different in actual stellar convection, where the development of magnetic structures takes different shapes (Stein and Nordlund, 2012;Masada and Sano, 2016;Käpylä et al, 2016b). Following Cameron and Schüssler (2015) (hereafter, CS15, also see the chapter by Cameron & Schüssler) we now estimate the budget of the toroidal magnetic flux in the dynamo region. Using the Stokes theorem and the induction equation Eq. (34), we define the derivative of the toroidal magnetic field flux in the northern hemisphere of the Sun as

Dynamo flux budget and impact of the surface activity on the deep dynamo
where Φ N tor = ∫ Σ dS, Σ is the meridional cut of the northern hemisphere of the solar convection zone, Σ stands for the contour confining the cut and the differential dl is the line element of Σ. The same can be written for the southern hemisphere flux Φ S tor . Similarly to CS15, we use the boundary conditions, and we estimate the RHS of the Eq. (60) in the coordinate system which is co-rotating with angular velocity of the solar equator, 0 = sin Ω 0 , and Ω 0 the surface angular velocity at the equator, here, = 0.99 , = 0.67 , are the radial boundaries of the dynamo region. In compare to CS15 we have additional contributions in the budget equation. Figure 16 shows profiles of the kernels 1−4 for the period of the magnetic cycle minimum. The estimations are based on results and parameters of the mean-field model presented above. Noteworthy, the south hemisphere should show the profiles of the opposite sign (see CS15). The results for 1,4 qualitatively similar to CS15. This is because the mean-field model show the qualitative agreement with solar observations for the time latitude evolution of the surface radial magnetic field. The diffusive decay of the toroidal magnetic flux is captured as well because of the phase shift between evolution of the poloidal and toroidal magnetic field in dynamo model and presumably in the solar dynamo as well. The model show the sharp poleward increase of 1 . This effect produces the winding of the toroidal magnetic field from poloidal component by the latitudinal shear. The effect of the radial shear, 2 , has maximum near the bottom of the convection zone, where its magnitude is less than the 1 . Figure 17(a) shows the time evolution of the RHS contributions of Eq. (61). In our model the We see that 2 is about factor 2 less than 1 . Winding of the toroidal field by the latitudinal shear seems to be the main generation effect in our model and, perhaps, in the solar dynamo, as well. The radial shear is less efficient because it is small in the main part of the convection zone. Also, it has the opposite sign near  the convective zone boundaries. This justifies applications of simple 1-D Babcock-Leighton dynamos to the solar observations as argued by CS15. Together with the fact of the poleward increase of 1 it explain the relative success of correlation of the polar magnetic field strength and the magnitude of the subsequent magnetic cycle for the solar cycle prediction (Choudhuri et al, 2007). Figure 17(b) shows the budget of the toroidal flux generation rate ( 1 + 2 ) and loss rate( 3 + 4 ) for our dynamo model. The parameters of the budget are larger than those deduced by CS15 from solar observations. The difference is because of additional generation and loss terms. The budget which includes only the surface activity contributions (green line in Fig.17b) is less than the full case. Also, the magnitude of the generation rate by the latitudinal shear can be larger than in the solar observations because of difference in the latitudinal profiles of the surface radial magnetic field. We guess that in the dynamo model the radial magnetic field increase poleward steeper than in observations. This issue have to be studied further. Figure 17b shows the budget for another two dynamo models. In one case, we neglect the generation effect of the radial subsurface shear in region r=0.9-0.99R. In compare to the standard case, this model shows the reduction of the generation rate, the amplitude of the generated toroidal flux, and increase of the loss rate. Therefore we conclude the importance of the subsurface shear layer for our dynamo model.
The above analysis shows the importance of surface activity for the dynamo model and perhaps for the solar dynamo as well. Sunspot activity in the form of magnetic bipolar regions is one of the most important aspects of magnetic surface activity. A consistent approach to include it in dynamo models is at present absent. Also, the origin of sunspots and their bipolar magnetic field is not well known; see Sect. 8.1. The Babcock-Leighton type and flux-transport dynamo models use a phenomenological approach. It is also applicable to mean-field models. Pipin (2022) studied the effect of surface activity on convection zone dynamos. Here, we briefly discuss some results of the paper. The emergence of bipolar magnetic regions (BMRs) is modeled using the mean electromotive force which is represented by the and magnetic buoyancy effects acting on the unstable part of the axisymmetric magnetic field as follows: where the first term takes into account the BMR's tilt and the second term models the surface magnetic region in the bipolar form. To produce the bipolar like regions we have to restrict spatially in Eq. (62) to the small scales that are typical for the solar BMR; see details in the above cited paper. Position and emergence time are chosen to be random and modulated by the large-scale magnetic activity. The BMR's -effect parameters are random as well; see details in (Pipin, 2022;Pipin et al, 2023). The given approach could be refined further using the 3D hydrodynamics, effects of the Coriolis force and the theory of the Joy's law developed recently by Kleeorin et al (2020). Figure 18 illustrates the formation of BMR simulated in the dynamo model. It was found that the BMR effects on the dynamo are restricted to the shallow layer below the surface. At the surface, the effect of the BMR on the magnetic field generation is dominant. Compare to the standard axisymmetric mean-field model discussed in the subsections above, the nonaxisymmetric dynamo, which includes the emergence of tilted BMR, can result in additional dynamo generation of the large-scale poloidal magnetic field and to an increase of the polar magnetic field. The red line in Fig. 17(b) shows the budget for this nonaxisymmetric dynamo model. We see an increase of the toroidal flux generation rate in the nonaxisymmetric model because of the surface BMR activity. Similar to Cameron and Schüssler (2015), we can conclude that sunspot surface activity seems to play an important part in the solar large-scale dynamo.

Effect of corona on the dynamo
Usually, dynamo models are limited to the star embedded in a vacuum, which is described by boundary conditions on the stellar surface. However, the boundary conditions have a determining influence on the global solutions, such as the symmetry about the equator. With the assumption of an external vacuum, all induction effects in the corona are neglected. Since the solar surface rotates differentially, the highly conductive plasma in the corona also causes induction effects through shear. Observations of coronal rotation are very scarce. There is evidence from extended coronal holes of rigid rotation in latitude (Timothy et al, 1975;Bagashvili et al, 2017). Kinematic dynamo models involving the corona with various assumptions on its rotation and conductivity give a wide range of solutions (Elstner et al, 2020). A notable influence of the corona on the dynamo in the convection zone was also observed in DNS by Warnecke et al (2016). A too weak density contrast and too strong viscous coupling of the corona to the star in their model probably underestimates the effect of the Lorentz force in the corona. Considering a dynamical situation with dominant Lorentz force in the corona, the solution in the Sun corresponds to that with vacuum boundary condition independent of rotation and conductivity in the corona. The magnetic field in the corona varies in time to a nearly force-free solution. Further investigations of the star-corona coupling are needed to clarify the exchange of magnetic energy and helicity. Noyes et al (1984) developed an early understanding of the observed stellar cycle periods, cyc . In those years, there where just six stars with measured rotation and cycle periods. Remarkably, for the stars of Noyes et al (1984), those values have not changed much with the more accurate data of Baliunas et al (1995); see Table 1 for a comparison of their cycle periods and the more recent data sets Boro Saikia et al, 2018;Bonanno and Corsaro, 2022). However, for the many stars of Baliunas et al (1995), the modern analysis of Olspert et al (2018) and Boro Saikia  Bonanno and Corsaro (2022) that were also included in the sample of Brandenburg et al (2017b). There are refer to the data of stars given in the Table 2. The dotted lines denote the fits determined by Brandenburg et al (2017b) while the upper (lower) solid lines denote fits to the stars of Bonanno and Corsaro (2022) with lowercase (uppercase) letters.

Stellar cycle periods
In recent work of Bonanno and Corsaro (2022), new cycle data were collected for altogether 67 stars. Their new sample includes stars with less accurate data points, so the existence of different branches was no longer a pronounced feature. In addition, many of the new data points are different from the earlier ones of Brandenburg et al (2017b); see Table 2. As in their paper, we denote G and F dwarfs by the same blue italic symbols and K dwarfs by the same red roman symbols.
To see how strong this revision of the data is, we plot in Fig. 19 the ratios rot / cyc versus log⟨ ′ HK ⟩ for all stars of Bonanno and Corsaro (2022) and highlight with lowercase and uppercase letters the stars that were also included in the sample of Brandenburg et al (2017b). We see that the new data are remarkably consistent with the old ones. Out of the eight stars on the branch of active stars, five where listed by Brandenburg et al (2017b) as having two periods. Of the 16 inactive stars, three were listed with two periods, but the case of the Sun was classified by Brandenburg et al (2017b) as somewhat different, because the 80 years Gleissberg cycle does not fit well on the active branch and, unlike all the other stars with two cycle periods, which are all younger than 3.3 Gyr, the Sun is relatively old.
The early data of Noyes et al (1984) did already suggest cyc ∝ Ω 1.25 for the cycle frequency cyc = 2 / cyc versus angular rotation rate Ω. This dependence is reproduced by considering free dynamo waves and assuming axisymmetric mean fields, =ˆ+ ∇ ×ˆ, with ( , ) ∝ i( − ) and writing −i = −i cyc + , where both cyc and are assumed to be real. The mean-field dynamo equations result in traveling wave solutions with a dispersion relation of the form At least up to moderate rotation rates, it is reasonable to assume that and Ω ′ are proportional to Ω. The crucial assumption in arriving at an approximation that matches Eq. (63) is to assume that the relevant wavenumber is selected not by the condition of marginal excitation, but by the assumption that = ( ) is maximized. Thus, has to obey d /d = 0, which yields cyc ∝ ( Ω ′ ) 2/3 ∝ Ω 4/3 . By contrast, if the dynamo is quenched to the being marginally excited, then cyc ≈ T / 2 , which would be either independent of Ω, or perhaps even decreasing with Ω, if T decreases with increasing Ω due to quenching. Of course, nonlinear dynamos must always be quenched to reach a steady state. This led Brandenburg et al (1998) to suggest that Eq. (63) could be obeyed if both and T are antiquenched in such a way that T is antiquenched more slowly than , so that cyc would increase with increasing magnetic field strength, and yet the dynamo would still saturate. Whether this is the only viable solution to this puzzle remains still unclear.
Recently, a number of numerical dynamo models were applied to investigate the relation of the cycle period on the stellar rotation rate in solar analogs (Pipin, 2015;Strugarek et al, 2017;Warnecke, 2018;Viviani et al, 2018;Hazra et al, 2019;Noraz et al, 2022). Figure 20 shows some of these results including the results of observations of Brandenburg et al (2017b) and the survey of Lehtinen et al (2016). It is interesting that the saturation branch of stellar activity for the young solar analogs with rotation periods of less than 10 days is well reproduced in the very different solar-like dynamo models including various GCD simulations (Strugarek et al, 2017;Warnecke, 2018;Viviani et al, 2018), flux transport model of Hazra et al (2019) and mean-field model of Pipin (2021). In Fig.20, this branch is marked by the green line.
The mean cycle period in this branch is almost independent of stellar rotation rate. The nonkinematic nonlinear model of Pipin (2021) shows multiple periods along this line. Pipin (2021) found that saturation of dynamo activity is accompanied by a depression of latitudinal shear, a concentration of magnetic activity to the surface, and changes in the meridional circulations from one cell to multiple cells per hemisphere structure. According to the conclusions of this paper, it is clear that in the saturated state, dynamo waves do not follow the Parker-Yoshimura rule. Their cycle period is determined by turbulent diffusion and meridional circulation. This is why the predictions of flux transport and nonkinematic mean-field dynamo models coincide. The independence of the cycle period on the rotation rate can be typical for the dynamo solutions which show concentration of the magnetic activity toward the boundaries of the dynamo (see Pipin, 2015;Pipin and Kosovichev, 2016 (2017b), green crosses those of Lehtinen et al (2016), orange squares the models of Warnecke (2018), and the asterisks are from the models of Pipin (2021); act/inact marks the active and inactive branches of activity, respectively; while kin/nkin stand for kinematic and nonkinematic models, respectively (adapted by permission from Pipin, 2021).
The inactive branch of the nonkinematic mean-field dynamo models shows a fairly strong positive inclination (see Fig. 20), which is absent in the kinematic models. We see that the dynamo model can reproduce a power law ∼ Co 0.5 , avoiding the antiquenching concept of Brandenburg et al (1998). In fact, the nonkinematic dynamo models show a frequency doubling phenomenon for models in the range from 10 to 15 days rotation period (see Figs. 3 and 8 of Pipin, 2021). The frequency doubling or the second harmonic generation is known from nonlinear optics. This is typical for wave propagation in nonlinear media. In dynamo waves, second harmonics are generated because of 2 effects such as the magnetic feedback on the large-scale flow, magnetic helicity conservation, and magnetic buoyancy effects. The second harmonics can be found in the solar activity as well . For the solar case, they are subdominant. However, they can become dominant for fast rotating stars. This makes the interpretation of magnetic activity cycles difficult (Stepanov et al, 2020).
It is important to note that the GCD simulations of Viviani et al (2018) show an increased level of nonaxisymmetry with an increase of rotation rate and a transition to preferred nonaxisymmetric dynamos for solar-type stars with a rotation period of less than 14 days. In simulations, this transition happens when the rotation profile changes from an antisolar to a solar-like profile with a faster equator and slow poles. Stellar Zeeman Doppler Imaging observations of Donati and Landstreet (2009) and See et al (2016), for example, show that the magnetic topology depends on stellar mass and rotation rate. In a certain interval of rotation periods (below 14 d) and stellar mass, the results of Viviani et al (2018) are compatible with the observational findings mentioned above.
Summarizing, we suggest that the Parker-Yoshimura dynamo regime can work for solar analogs rotating with periods above 14-15 days. In an interval of stellar rotation periods between 10 and 15 days, frequency doubling and a transition to the nonaxisymmetric dynamos occurs. For lower rotation periods, the dynamo transits to a saturation stage. It can be characterized by high magnetic activity and multiple dynamo periods which are independent of stellar rotation rate. This new picture should be further improved by including possible dynamo effects of surface activity in the form of bipolar magnetic regions and star spots.
Another interesting observation is that different types of stars, including partially convective main sequence solar-type stars, fully convective M-dwarfs, and evolved post-main sequence giant stars show similar scaling with the Rossby number for the unsaturated regime (see, e.g., Wright and Drake, 2016;Wright et al, 2018;Lehtinen et al, 2020). In particular, to derive the Rossby number, these study employ a simple parameterization of the convective turnover time suggested by stellar interior models. It was found that for evolved stars, the Rossby-independent parameterizations break down in the rotation-activity relation (Lehtinen et al, 2020). This constitutes a strong argument in favor of the turbulent dynamo paradigm suggesting a common role of the turbulent process in dynamos operating in these stars.

Analogy of mean-field and the chiral magnetic effect
The effect in mean-field dynamo theory is an effect that emerges after averaging over the scale of several turbulent eddies. We know already that turbulent diffusion is somewhat analogous to microphysical diffusion, which also emerges after averaging, but here after averaging over atomic scales. Interestingly, even for the effect there can be an effect on atomic and subatomic scales, because fermions, such as electrons, are chiral. The spin of an electron emerging from the decay of a neutron is anti-aligned with its momentum vector, so their dot product is a negative pseudo-scalar, called the chirality. Positrons have positive chirality. In the presence of an ambient magnetic field, the spins align, but electrons and positrons move in opposite directions, causing therefore an electric current. This constitutes a microscopic effect (Rogachevskii et al, 2017;Brandenburg et al, 2017c), where 5 is the normalized chiral chemical potential (with units of inverse length), is the microscopic magnetic diffusivity, fine ≈ 1/137 is the fine structure constant (quantifying the strength of electromagnetic interaction between charged particles), is the reduced Planck constant, ≈ 3 × 10 10 cm s −1 is the speed of light, B ≈ 10 −16 erg K −1 is the Boltzmann constant, and is the temperature.
The applications of chiral MHD are manifold and range from condensed matter systems and heavy ion collisions to neutron stars and the early Universe; see Kharzeev (2014) for a review. Interestingly, because this microscopic effect produces helical magnetic fields, and because the total chirality is conserved (Boyarsky et al, 2012), this effect does not last forever, but is being quenched in a form analogous to the catastrophic quenching formula, which takes the form (Rogachevskii et al, 2017) where is a coupling constant which, in the catastrophic quenching formalism, is related to 2 t 2 f / 2 eq , and Γ flip is a spin-flipping parameter, which is related to 2 2 f in the catastrophic quenching formalism (see, e.g., Field and Blackman, 2002;Blackman and Brandenburg, 2002).
There is a vast range of recent work in this field, which goes well beyond the scope of the present paper. We just mention here the paper of Masada et al (2018), who studied chiral magnetohydrodynamic turbulence in core-collapse supernovae. They found that the inverse cascade related to the chiral effects impacts the magnetohydrodynamic evolution in the supernova core toward explosion.

Direct Statistical Simulations
An alternative or extension to mean-field theory in the usual sense is to solve the time-dependent system of one-point and two-point correlation functions. This is what is now known as Direct Statistical Simulations (Tobias et al, 2011;Marston, 2013, 2017) and has primarily been applied to two-dimensional turbulent shear flow problems. The dimensionality of the two-point correlation function doubles for those directions over which homogeneity cannot be assumed. On the other hand, the dynamics of the low order statistics is usually slower than that of the original equations. In addition, it is possible to reduce the complexity of the problem by employing Proper Orthogonal Decomposition (Allawala and Marston, 2016;Allawala et al, 2020). This approach has not yet been applied to magnetohydrodynamics and the dynamo problem. Such an approach would be able to incorporate both smallscale and large-scale dynamos at the same time. The small-scale dynamo problem would be solved through correlation functions, as has been done in Brandenburg and Subramanian (2000), for example. But their simulations were isotropic and did not include anisotropic dependencies on position. Nevertheless, this approach has the potential of being a strong competitor in addressing the high Reynolds number dynamics of problems of astrophysical and geophysical relevance; see Tobias (2021) for a recent review touching on these aspects.

Looking forward
In this review, we have provided some insight into recent developments in our understanding of the generation of astrophysical large-scale magnetic fields. The current development of mean-field theory allows us to go beyond some of the original restrictions that were related to the assumption of large scale-separation and the inappropriate neglect of nonlinear effects due to higher order correlations in contributions to the mean turbulent electromotive force. A big portion of the progress comes from the development in the DNS of astrophysical turbulence. Noteworthy, the classical mean-field theory is based on the fundamental equations of electrodynamics and has well-known limits. With the new steps forward, we can take into account results of the DNS, e.g., the spectral kernels, and treat them as experimental facts. The necessity of some phenomenological additions to classical mean-field theory are motivated both by DNS and observations of the magnetic activity in astrophysical systems, such as those in our Sun and other stars. In this way, mean-field models become a valuable tool to understand the real and virtual worlds of the dynamo in stars and in DNS.