Dynamics of Large-Scale Solar Flows

The Sun’s axisymmetric large-scale flows, differential rotation and meridional circulation, are thought to be maintained by the influence of rotation on the thermal-convective motions in the solar convection zone. These large-scale flows are crucial for maintaining the Sun’s global magnetic field. Over the last several decades, our understanding of large-scale motions in the Sun has significantly improved, both through observational and theoretical efforts. Helioseismology has constrained the flow topology in the solar interior, and the growth of supercomputers has enabled simulations that can self-consistently generate large-scale flows in rotating spherical convective shells. In this article, we review our current understanding of solar convection and the large-scale flows present in the Sun, including those associated with the recently discovered inertial modes of oscillation. We discuss some issues still outstanding, and provide an outline of future efforts needed to address these.


Solar Differential Rotation and Meridional Circulation
Differential Rotation Due to the solar rotation, axisymmetry (about its rotational axis) of the large-scale flows is established to some degree in the Sun's interior.Differential rotation denotes the longitudinal component of the axisymmetric (longitudinally-averaged) flow which varies with radius and latitude.It arises from the nonlinear interaction of the rotationallyinfluenced solar magneto-convection (e.g., Miesch 2005).Differential rotation represents a shear in the rotation rate and is thought to play a significant role in the solar dynamo by stretching and amplifying the magnetic field lines (Charbonneau 2020).
The solar differential rotation profile can be measured by global helioseismology, which analyzes small frequency splittings of resonant acoustic oscillations (global standing acoustic modes) (Duvall et al. 1984;Thompson et al. 1996;Schou et al. 1998;Howe et al. 2000).Figure 1 shows the observationally-inferred profile of the internal differential rotation of the Sun (Howe et al. 2005;Larson and Schou 2018).We summarize striking features of the solar differential rotation as follows: Solar and Stellar Dynamos: A New Era Edited by Manfred Schüssler, Robert H. Cameron, Paul Charbonneau, Mausumi Dikpati, Hideyuki Hotta and Leonid Kitchatinov Extended author information available on the last page of the article • The radiative interior rotates almost rigidly.
• In the convection zone, the equator rotates about 30% faster than the poles.
• The transition from uniformly-rotating radiation zone to differentially-rotating convection zone occurs in a thin layer from 0.68 R to 0.73 R .This layer is called the tachocline.• In the bulk of the convection zone (0.73 R < r < 0.96 R ), the rotation rate varies strongly with latitude and much more weakly with radius.• In a shallow surface layer (r 0.96 R ), the rotation rate decreases by about 5% at all latitudes.This layer is called the near surface shear layer.• The contours of constant angular velocity are inclined by about 25 • with respect to the rotational axis over a wide range of latitude.In other words, the differential rotation does not follow the Taylor-Proudman theorem.
These observational facts need to be explained by theoretical and numerical models of rotating solar magneto-convection.
Meridional Circulation Meridional circulation represents radial and latitudinal components of the large-scale axisymmetric flow in the Sun, i.e., a poloidal flow in a meridional plane.Meridional circulation, as well as the differential rotation, is believed to play a significant role in the solar dynamo by advecting the magnetic flux in both radial and latitudinal directions (e.g., Charbonneau 2020).The meridional circulation is much weaker than the differential rotation (two orders of magnitudes smaller in flow amplitude) and cannot be inferred using conventional globalmode helioseismology.Therefore, it is an extremely difficult task to measure the meridional flow in the Sun.Near the solar surface, the meridional flow is poleward in both hemispheres with typical amplitudes of ∼10-20 m s −1 .This was first measured by Duvall (1979) using Doppler measurements and then robustly confirmed in follow-up studies by a variety of methods (e.g., Patron et al. 1995;Giles et al. 1997;Hathaway 1996;Braun and Fan 1998;Haber et al. 2002;Ulrich 2010;Basu and Antia 2010).
Local helioseismology can extend these measurements into the deeper convection zone (e.g., Gizon and Birch 2005).In particular, time-distance helioseismology (Duvall et al. 1993) and ring-diagram analysis (Hill 1988) can be used to measure the effects of the meridional flow on north-south propagating waves.In time-distance helioseismology, wave travel times are extracted from the cross-covariance of the Doppler signals between points along meridians.The technique was applied first by Giles et al. (1997) to the Michelson Doppler Imager (MDI) data onboard SOHO (Scherrer et al. 1995).However, this is an extremely difficult measurement to make because the deep meridional circulation is very weak (no greater than 3-5 m s −1 ) and the sensitivity of the travel times to flows also decreases with depth.The analysis of very long time series is required to reduce noise.Furthermore, for accurate measurements it is critically important to apply corrections to the measurements, especially a center-to-limb correction (Duvall and Hanasoge 2009;Zhao et al. 2012), and corrections for the P -angle and the Carrington elements (e.g.Giles 2000;Hathaway and Rightmire 2010;Liang et al. 2017).Often pixels that are in regions of very strong magnetic fields (e.g. in sunspots) are excluded from the cross-covariances (e.g.Liang and Chou 2015).The center-to-limb effect is, in general, very significant and depends strongly on time and instrument (Liang et al. 2018).It may have both an instrumental and physical component.When Giles et al. (1997) made their measurements, no center-to-limb correction was applied to the MDI travel times, as it was very small during the year 1996 (Liang et al. 2018).We refer the reader to the review by Hanasoge (2022) for additional information.
Many different inferences of the solar meridional circulation using time-distance helioseismology have been published.They are not consistent.Using travel time measurements obtained from SDO/HMI in 2010-2012, Zhao et al. (2013) reported an equatorward return flow in the middle of the convection zone (0.82-0.91 R ) and a poleward flow below 0.82 R , indicating a double-cell structure in radius (Fig. 2a).A similar result was obtained by Chen and Zhao (2017) who used seven years of data from HMI and active regions were masked out to remove the contribution from pixels with strong magnetic fields.Note, however, that Zhao et al. (2013) and Chen and Zhao (2017) did not invert the radial component of the meridional flow which is required by the local mass conservation.An attempt to test a mass conservation constraint was made by Jackiewicz et al. (2015) who found that the equation of continuity is poorly satisfied for the inverted flows from GONG and HMI (only two years of data were used).Under the constraint of mass conservation in terms of the stream function, Rajaguru and Antia (2015), using four years of HMI data report a singlecell meridional circulation in each hemisphere, with an equatorial return flow near the base of the convection zone (below 0.77 R ).A similar result was reported by Mandal et al. (2018).Liang et al. (2017) have compared the north-south travel-time data from SOHO/MDI and those from SDO/HMI over the 1-year overlap period of 2010.After correcting for the center-to-limb effect (Zhao et al. 2012), it is found that the travel-times are consistent within the error bars, although an overlap period of only 288 days is far from enough for such consistency test.When comparing the travel-times with different forward models of meridional flows, Liang et al. (2018) found that the MDI/Cycle 23 data point to a single cell meridional flow in both hemispheres, while the HMI/Cycle 24 data point to a double cell in the north and a single cell in the south.
Following this work, Gizon et al. (2020a) carried out comprehensive measurements of the north-south travel-time measurements using all available data sets from GONG  Gizon et al. (2020a) using GONG data (2008-2019), reprinted with permission from AAAS (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019), SOHO/MDI (1996-2011), and SDO/HMI (2010-2019).After correcting for all known systematic errors (CCD orientation, center-to-limb effect, Carrington elements, pixels in sunspots and active regions), it was found that the north-south travel times measured from GONG and MDI are in good agreement over the overlap period of 1996-2011.However, a small travel-time offset is apparent in the HMI data compared to the GONG data over the period 2010-2019.No explanation was found for this HMI offset and the data were set aside: The HMI measurements imply a strikingly different meridional flow pattern between the northern and southern hemispheres, as hinted at by Liang et al. (2018) and confirmed by Gizon et al. (2020a) and Braun et al. (2021).Gizon et al. (2020a) find that the mass-conserving meridional flows during cycle 23 (from MDI/GONG data) and cycle 24 (from GONG data) display a single-cell pattern in each hemisphere, as shown in Fig. 2b.
Work is needed to resolve the remaining issues that affect helioseismic inferences of the deep meridional flow.Efforts are ongoing to identify the source of the HMI travel-time offset, which is extremely small (0.2 s) but significant.The other pressing issue is to gain an understanding of the origin of the center-to-limb effect (e.g.Chen and Zhao 2018).The analysis of travel-times in different frequency bands is a promising avenue to make progress in this area (Chen and Zhao 2018;Rajaguru and Antia 2020).
Observational determination of the solar meridional circulation is crucial, not only for constraining solar dynamo models (e.g., Hazra et al. 2014), but also for properly understanding the angular momentum flux balance in the Sun's convection zone (e.g., Featherstone and Miesch 2015).
Cycle Variations of Large-Scale Mean Flows Solar Doppler observations and global helioseismology have revealed that the Sun's differential rotation exhibits a temporal variation pattern consisting of multiple bands of faster-and slower-than-average zonal flows which migrate in latitude with the phase of the solar dynamo cycle (e.g., Howe 2009).The so-called torsional oscillation (coined by Howard and Labonte 1980) has two distinct branches, one at active-region latitudes and one at latitudes above 50 • .The low-latitude branch migrates equatorward together with the magnetic activity belts and shows a clear 11-yr periodicity.On the other hand, the poleward-migrating high-latitude branch tends to show a rather irregular behavior: it was clearly seen in Cycle 23 but did not appear in Cycle 24 (e.g., Howe et al. 2013).For more details on the observational aspects of the Sun's torsional oscillations, see Norton et al. (2023) (their Sect. 6.2).The amplitude of the solar torsional oscillation is ≈ 3-5 nHz and it is important to understand its physical origin as it likely reflects the nonlinear feedback of the dynamo-generated magnetic fields onto the large-scale flows (e.g., Rempel 2007;Beaudoin et al. 2013;Pipin and Kosovichev 2019;Brun et al. 2022).
As is the case of differential rotation, the poleward meridional circulation shows a temporal variation associated with the solar magnetic cycle (Beck et al. 2002;Gizon 2004;Gizon et al. 2010;Mahajan et al. 2023); see also Sect. 6.3 in Norton et al. (2023).This variation is of order ≈5 m s −1 or more, which is a significant fraction of the maximum meridional flow amplitude.It may be explained, at least in part, by the north-south component of the inflows around active regions (Gizon et al. 2010;Mahajan et al. 2023).In addition to this component, Mahajan et al. (2023) find that there is a residual solar-cycle component as small as ≈2 m s −1 , which is seen around cycle minima.

Solar Convective Flows
Convection on the Sun occurs over a wide range of spatial scales, and while the spectrum is continuous, apparent characteristic scales are commonly cited: granulation, mesogranulation, supergranulation, and giant cells.Granulation (Herschel 1801) is readily apparent in high-resolution images of the solar photosphere, as a pattern of bright upflowing regions separated by darker downflowing lanes.The characteristic upflow cells have diameters of ∼1000 km, lifetimes of about 0.2 hr, and vertical flow speeds of ∼1 km s −1 .The upflow velocity often peaks near the granular boundaries (e.g., Nesis et al. 1992;Rast 1995;Hirzberger 2002;Nordlund et al. 2009;Falco et al. 2017, and reference therein).These properties reflect the compressible flow dynamics of a strongly cooled radiative boundary layer, with observations confirming the convective nature of the flow via measurement of the correlation between the vertical velocity and plasma temperature (e.g., Canfield and Mehltretter 1973).Granulation is well observed and robustly modeled (se e.g., Nordlund et al. 2009), even in quite shallow domains, by codes that capture the rapid change in radiative opacity in the solar photosphere and implement an open lower boundary condition to minimize bottom-up influences on the top-down dynamics of the radiative boundary layer.
Mesogranulation (November et al. 1981), on the other hand, is observationally elusive.With a reported length scale of about 5-10 Mm, ∼60 m s −1 vertical flow speeds, and ∼2-3 hr lifetime, its identification as a convective feature is still debated.Most recent studies suggest that no distinct mesogranular scale is present in the broad range of convective scales observed (e.g., Rincon and Rieutord 2018, and references therein).One possibility is that there is weak advective self-organization of the granular flows, a process first proposed in the context of supergranulation (Rieutord et al. 2000;Rast 2003) but likely more relevant on mesogranular scales (Cattaneo et al. 2001;Berrilli et al. 2005;Leitzinger et al. 2005;Duvall and Birch 2010).However, the absence of a mesogranular scale in the clustering of magnetic elements in high resolution magnetograms suggests that this mechanism too leads to a continuous exponential distribution of scales between 2 and 10 Mm, with no distinctive characteristic peak (Berrilli et al. 2013).Supergranulation (Hart 1954;Leighton et al. 1962) is the largest likely-convective scale of motion readily visible in the solar photosphere.It is observed directly in spectral Doppler shifts away from disk center (due to horizontal motions) and is traced by network magnetic elements which are prominent in magnetograms and in emission in low chromospheric lines such as Ca II K.There is good correlation between Ca II K emission and magnetic flux density (Ortiz and Rast 2005, and references therein).Supergranular cells have diameter of ∼30 Mm, horizontal flow velocities of ∼100 m s −1 , and lifetimes of ∼20 hr.After the intensity contributions of the small scale magnetic field elements has been removed, they show an average continuum intensity contrast across the cells of about 0.1%, corresponding to about one degree Kelvin in brightness temperature (Goldbaum et al. 2009).
The origin of the supergranular motions has been widely debated (see Rincon and Rieutord 2018, and references therein).It has recently been proposed that the scale of supergranulation reflects not a selected convective scale, but is instead defined by the scale above which convective power declines (Lord et al. 2014;Cossette and Rast 2016).This interpretation, and the reasons underlying the power reduction, links the well observed phenomenon of supergranulation to the convective conundrum, an outstanding discrepancy between models and observations (see this section below, and Sects.2.2 and 3.1).We note that the early suggestion that helium ionization plays a role in determining the mesogranular and supergranular scales (Leighton et al. 1962;Simon and Leighton 1964) is not supported by numerical simulations or simplified models base on them (Rast and Toomre 1993;Lord et al. 2014).Additionally, the presence of the network magnetic elements themselves (Crouch et al. 2007;Thibault et al. 2012) or the enhanced radiative loses through them (Rast 2003) do not seem to play a role in scale selection role (Lord 2014).Finally, it is important to note that supergranulation shows peculiar unexplained wave-like properties (Gizon et al. 2003;Schou 2003;Langfellner et al. 2018).
In contrast with supergranulation, which is readily observed but not captured by any local-area or global spherical-shell simulation, solar giant cells (Simon and Weiss 1968), motions on the scale of the solar convection zone depth (∼200 Mm), dominate global spherical-shell simulations but are very difficult to observe (Hathaway et al. 2013; Hathaway and Upton 2021, and references therein).If, in the Sun, giant-cells had the amplitude they do in simulations, they would be easily observed in the solar photosphere.This is the simplest manifestation of the convective conundrum: that supergranulation, rather than giant-cell scale motions, are the largest readily observed motions in the solar photosphere.The implication for solar differential rotation is fundamental.The enhanced amplitude of the large-scale convective motions in global numerical simulations tend to place those simulations in a Rossby-number regime that favors anti-solar differential rotation profiles.
These issues are critical to our understanding of large-scale motions on the Sun.As can be seen in Fig. 3 (from Hathaway et al. 2015), only two of the components described above are evident as distinct features in the observed spectrum of motions in the solar photosphere.Granulation is responsible for the most pronounced peak at high spherical-harmonic degree and supergranulation for the smaller peak near spherical-harmonic degree 120.Added to the plot are vertical fiducial lines indicating the approximate scale of supergranulation and giant cells.Additionally, a blue dotted line has been added to schematically indicate the monotonic increase of power to low wavenumbers seen in all realistic numerical simulations up until the most recent of Hotta and Kusano (2021).In the Hotta and Kusano (2021) simulations in the power rolls over at spherical-harmonic degree ∼10.It is the discrepancy in low spatial-frequency power between simulations and observations that has come to be known as the convective conundrum.
It is important to note that the spectrum plotted in Fig. 3 is a composite, with vertical velocities dominating at high spatial wavenumbers (granular scales) and horizontal motions most important at supergranular scales.The vertical velocity contribution decreases from the granular peak towards lower wavenumbers, with horizontal velocity contribution increasing to spherical-harmonic degree ∼120 before rolling over beyond that.The supergranular peak results from this decrease in the power beyond spherical-harmonic degree ∼120.Thus, with respect to photospheric flow observations, the convective conundrum refers to the scale and amplitude of the horizontal-flows in the photosphere.No global-spherical-shell or localarea simulation of solar convection yet captures the supergranular scale maximum in photospheric horizontal-flow power.

Solar Inertial Modes
On scales much larger than the supergranules, the non-axisymmetric flows have long timescales and are thus strongly influenced by solar rotation via the Coriolis force.A component of these flows has recently been observed as retrograde waves of radial vorticity at the surface and are known as inertial modes.In this section, we give an overview of their characteristics and properties.
Inertial modes owe their existence to rotation (Greenspan et al. 1968).Their restoring force is the Coriolis force.In a uniformly rotating sphere, the frequencies of inertial modes are limited to a range of |ω| < 2 0 in the co-rotating frame, where 0 denotes the angular velocity.The traditional Rossby modes (or r modes) correspond to a variety of inertial modes that have quasi-toroidal motions.Although these modes have been expected to exist in the Sun and stars since the late 1970's (e.g., Papaloizou and Pringle 1978;Saio 1982;Unno et al. 1989), they were not observed on the Sun until very recently.Inertial modes on the Sun have very long oscillation periods (of the order of months) and very small velocity amplitudes (of the order of 1 m s −1 ).Therefore, long-term and high-precision observations of horizontal flows over many years are required to detect them.Löptien et al. (2018) discovered the solar equatorial Rossby modes using both a granulation-tracking method and a ring-diagram analysis applied to six years of SDO/HMI data.They detected excess power in the radial vorticity along the dispersion relation of the sectoral Rossby modes, ω = −2 eq /(m + 1) in the Carrington frame, for azimuthal orders 3 ≤ m ≤ 15.In this formula, eq is the equatorial rotation rate at the Sun's surface.For m 5, the latitudinal eigenfunctions of the equatorial Rossby modes significantly deviate from the sectoral spherical harmonics: the radial vorticity peaks at the equator but changes sign at middle latitudes (this is due to differential rotation, see below).The detection of solar equatorial Rossby modes has been confirmed in follow-up studies using various other observational datasets and methods (Liang et al. 2019;Proxauf et al. 2020;Mandal and Hanasoge 2020;Hanson et al. 2020;Hathaway and Upton 2021).
The observed equatorial Rossby modes exist only at very large scales with m ≤ m crit ≈ 15.Löptien et al. (2018) speculated that this critical azimuthal order m crit might correspond to the Rhines scale l Rhines = √ R v c / above which rotation strongly affects turbulent convection (Rhines 1975).Assuming m crit ≈ R /l Rhines , a typical speed of turbulent convection can be roughly estimated as v c ≈ 10 m s −1 , which is about one order magnitude smaller than the typical mixing-length estimate (Böhm-Vitense 1958;Stix 2002).This may add another piece of information with regard to the solar convective conundrum (see Sects.2.2 and 3.1 below).
Recently, Gizon et al. (2021) analyzed more than 10 years of data from both SDO/HMI and GONG and detected many additional quasi-toroidal inertial modes with 1 ≤ m ≤ 10.With the help of a 2.5D linear eigenvalue solver applied to a model of the differentially rotating convection zone (Bekki et al. 2022b), they identified modes at middle and high latitudes, in addition to the equatorial Rossby modes.The observational power spectra and the measured eigenfunctions of three selected inertial modes with m = 1, 2, and 3 are shown in Fig. 4.These modes are very sensitive to the solar latitudinal differential rotation (see Gizon et al. 2020b, for a discussion in the β plane).For many of the inertial modes, there exists latitudes at which the phase speed is equal the local differential rotation speed; such latitudes are called critical latitudes (see Fig. 4, middle column).The m = 1 high-latitude mode has a large amplitude (v φ ≈ 15 m s −1 ) at high latitudes (above ≈ 50 • ) and its surface eigenfunction exhibits a spiral pattern in the polar regions.The horizontal flow associated with the m = 1 high-latitude mode was first observed by Hathaway et al. (2013), but misidentified at the time as giant cell convection.It had also been observed by Bogart et al. (2015) and Howe et al. (2015).A linear model suggests that the m = 1 high-latitude mode is self-excited when a large enough latitudinal entropy gradient exists in the convection zone (Bekki et al. 2022b).
More recently, Hanson et al. (2022) reported the detection of another class of inertial modes with north-south anti-symmetric radial vorticity across the equator.These modes propagate in a retrograde direction with the phase speed roughly three times faster than the that of the equatorial Rossby modes.Compared with the equatorial Rossby modes and the high-latitude modes reported in Gizon et al. (2021), these modes have lower velocity amplitudes and are thus much harder to distinguish from the background noise in the power spectra.According to simplified linear eigenmode calculations of a uniformly-rotating Sun (Triana et al. 2022;Bhattacharya and Hanasoge 2023), these modes are not quasi-toroidal, i.e., substantial radial motions are involved.
Though it remains unclear how important the solar inertial modes are to the overall convection zone dynamics, they play an important diagnostic role.Gizon et al. (2021) and Bekki et al. (2022b) have shown that the properties of some inertial modes (i.e., frequencies, linewidths, and surface eigenfunctions) are sensitive to the turbulent viscous diffusivity ν t and to the superadiabaticity δ of the convection zone.Gizon et al. (2021) inferred that, on average, ν t ≤ 10 12 cm 2 s −1 and δ < 2 × 10 −7 .These values are an order of magnitude smaller than the theoretical estimates from the local mixing length model (Christensen-Dalsgaard et al. 1996;Muñoz-Jaramillo et al. 2011).It is noteworthy that both of these parameters cannot be constrained by conventional p-mode helioseismology, and are important in discussions of the convective conundrum (Sects.2.2 and 3.1 below) and the solar The same power spectra but normalized at each latitude by their average value over the frequency range between the orange bars.Excess power is seen at a specific frequency at all latitudes in each of the three cases.The red arrows point to critical latitudes at the surface for the case m = 2. Bottom row: Observed horizontal velocity eigenfunctions at the surface for the m = 1 high-latitude mode, the m = 2 critical-latitude mode, and the m = 3 equatorial Rossby mode.These figures are courtesy of Gizon et al. (2021) dynamo.The amplitudes of the linearly-stable modes might provide additional constraints on the turbulent viscosity (Philidet and Gizon 2023).On the other hand, to better understand the amplitudes of the linearly-unstable inertial modes, nonlinear numerical simulations will be required (Bekki et al. 2022a;Matilsky et al. 2022;Bekki and Cameron 2023).

Observations of Large-Scale Flows on Solar-Type Stars
It is possible to obtain some information about the rotation of distant stars using observations of photometric variability, radial velocity measurements, and asteroseismology.Such methods are challenging, however recent results are setting new constraints on stellar differential rotation.
One possibility is to measure the effects of rotation on stellar oscillations.Rotation lifts the degeneracy of the frequencies of the modes of oscillation with different azimuthal orders.
Latitudinal and radial differential rotation may leave a signature in the fine structure of the oscillation power spectra.See, e.g., Aerts et al. (2010) and García and Ballot (2019) for reviews of asteroseismology.In the last years it has become possible to probe the deep interior of some evolved stars, including red giants (Beck et al. 2012;Deheuvels et al. 2012) and, more recently, subgiants (Deheuvels et al. 2020).Probing the slowly-rotating Sun-like stars is especially difficult because rotational splitting frequencies are often too small compared to the mode linewidths.Constraints on the mean angular velocity and the inclination angle of the rotation axis have been obtained in a few cases (e.g., Gizon et al. 2013;Chaplin et al. 2013).Strong latitudinal differential rotation has recently been observed on a few selected main-sequence solar-type stars (Benomar et al. 2018).Solar-like differential rotation profile (equator faster than the poles) is detected on 13 stars with a much stronger amplitude than on the Sun (on average twice the solar value).Subsequently, Bazot et al. (2019) reported latitudinal rotation contrasts closer to the solar value for the two solar analogs 16 Cyg A and B. The equator-to-pole difference in angular velocity, , may depend on the mean rotation rate of the stars, with the majority of stars in Benomar et al. (2018)'s study rotating faster than the Sun.Only upper limits have been set on the radial differential rotation of Sunlike stars (Nielsen et al. 2014).To understand differential rotation in stars, it is important to study its relationship to the mean rotation rate * , the stellar mass M * , the stellar age, and the stellar composition.
Surface rotation and differential rotation can also be studied using observations of stellar photometric variability.The presence of starspots and faculae on rotating stars modulate the photometric signal at the rotation period (e.g., Nielsen et al. 2013;Reinhold et al. 2022).Some information about latitudinal differential rotation can be retrieved when magnetic active regions are present over a range of latitudes on the stellar surface (e.g., Reinhold et al. 2013;Nielsen et al. 2019).The stellar differential rotational versus mean rotation rate may be described by a power law, ∝ n * .Figure 3 of Barnes et al. (2005) summarizes the different studies that were available at the time.The long-term monitoring of photometric modulations point toward a weak rotational dependence, n = 0.24 (Henry et al. 1995), while similar observations in H & K bands of Ca II emission suggest a significantly higher value, n = 0.7 (Donahue et al. 1996).The spectroscopic analysis of Reiners and Schmitt (2003) gave a power law index of n = 0.66.Barnes et al. (2005) proposes an index of n = 0.15, but also show that this index is very sensitive to the spectral-type diversity of the target sample considered (see also Reiners and Schmitt 2003).For instance, Reinhold et al. (2013) find a value of n = 0.3 for cool stars and Balona and Abedigamba (2016) report n = 0.2 for G-stars.Asteroseismic studies tend to find values for n between 0.3 and 0.45 (García and Ballot 2019), and underline its sensitivity to the effective temperature T eff range of the stars considered (see Reinhold and Gizon 2015, and references therein).
The dependence of the latitudinal contrast on spectral type also appears to be significant for fast main-sequence rotators.Collier Cameron ( 2007) studies fast M, K, G, and F main-sequence rotators and reports a strong ∝ T 8.6 eff dependence on the effective temperature, consistent with Barnes et al. (2005).Large error bars are associated with the hottest spectral types, however the models of Kueker and Ruediger (2011) tend to confirm a strong dependence (see Reinhold and Gizon 2015, for a detailed discussion).
To take into account both rotational and spectral-type aspects, Saar (2010) proposed to consider the global shear as a function of the stellar Rossby number Ro s = τ c / * (see also Brun et al. 2017 andNoraz et al. 2022a for definitions and prescriptions).Indeed, it is possible to parameterize the convective turnover time τ c as a function of T eff (Cranmer and Saar 2011).In particular, he finds that ∝ Ro −1 s for unsaturated rotators ( ≤ 12 ), pointing then toward n = 1 when fixing T eff and the composition.Numerical studies of the last decade have confirmed that the Rossby number is a major parameter to consider regarding the characterization of large-scale flows along stellar evolution (see Sect. 2.6).However, it has to be mentioned here that the high-Rossby regime still needs observations to constrain theoretical results.A difficulty lies here in the sensitivity of the aforementioned techniques to the rotation rate, which decreases signals significantly when considering slowrotators.
Finally, recent observational studies have started to investigate the impact of the metallicity, for instance with the solar analog HD 173701 studied by Karoff et al. (2018).Its parameters are indeed close to the solar ones, while having a significantly higher metallicity ([Fe/H] = 0.3 ± 0.1).Using different methods previously mentioned, the authors report a solar-like differential rotation (fast equator) with a latitudinal contrast being twice the solar one.Monitoring of its chromospheric and photometric emissions also show a cyclic activity shorter than the one observed on the Sun (P cyc = 7.4 against 11 years), while having a higher amplitude of variation, which underlines the entanglement between composition, large-scale flows, and magnetism (Brun and Browning 2017;See et al. 2021).
To summarize, differential rotation is a characteristic quantity of stellar convective envelopes.Apart from the Sun, the exact quantification of the surface rotational contrast remains difficult because it requires a high degree of precision of the instruments used and long acquisition periods for the different targets.New observations by the upcoming PLATO mission should allow major advances in this direction (Rauer et al. 2014).A particular focus for future prospects will lie on the monitoring and characterization of slow-rotators, which appear to be challenging targets for current techniques (Noraz et al. 2022a;Donati et al. 2023, see also Sect.2.6).In the meantime, theoretical modeling of these flows are hence crucial for the understanding and support of these observations, and in order to guide those to come.

Models of Large-Scale Flows
Large-scale flows in the Sun are generated and maintained by thermal convection.In this section we briefly summarize the current theoretical understanding of how that occurs.

Mixing Length Theory and Energy Budget
Mixing-length theory (MLT; Böhm-Vitense 1958) remains a remarkably useful way to describe the mean energy transport by convection even when the actual dynamics is far from localized eddy motions.In the mixing length formulation, we can relate energy flux to the convective velocity and the stratification through a parameter called the mixing-length parameter α MLT = L MLT /H p , where L MLT and H p = − (d log p/dr) −1 are the mixing length and the pressure scale height in the convecting fluid.The mixing-length parameter, α MLT is a parameter of order unity which mainly affects the stellar radius (Demarque and Percy 1964) when used in a stellar structure model.Here we take α MLT = 1 (i.e., L MLT = H p ) for the simplest mixing-length formulation.The typical temperature perturbation in the convection T is then evaluated as: where δ = (d log T /d log p)−(d log T /d log p) ad is the superadiabaticity, the deviation from the adiabatic stratification.When the speed of sound in the medium is much faster than convection, the sound wave instantaneously relaxes any pressure perturbation p, and the density perturbation ρ can be estimated from the linearized equation of state as (2) A simplified equation of motion is converted with rough dimensional analysis as where v c is the typical convection velocity.Here we can evaluate the typical time for the convection as τ c = H p /v c .Then the convection velocity can be written as where c s is the speed of sound.Together, these approximations yield several important approximations for the superadiabaticity: While some more careful treatments include additional physical effects, such as variation in the mean molecular weight of the plasma or the fluid drag force, in the practical application of MLT, the relationships captured by Equation (5) change little.What is important is that these relations lead to a ratio of the kinetic energy E kin to the internal energy E int that scales with superadiabaticity δ: In the convection zone, heat is mainly transported by the convection, with the enthalpy flux given by Normalizing this equation yields Over the depth of the convection zone, with the superadiabaticity being 10 −6 and 10 −1 at the base and surface respectively (e.g., Ossendrijver 2003), the convective velocity amplitudes should change by a factor of a few hundred.At the base of the solar convection zone, the flow is subsonic (Mach number ∼ 10 −3 ) and the internal energy of the fluid is much larger (10 6 times larger) than the kinetic energy of the flows.

The Effect of Stratification on Observed Horizontal Convective Flow Amplitudes
As is clear from the mixing-length calculation above, one of the most important aspect of stellar envelope convection is the steep stratification of the mean state.The solar convection zone is about 210 Mm deep, over which the density changes by a factor of about one million and the pressure by about 800 million.The density scale height in the photosphere is about 150 km, while at the convection zone base it is equal to nearly half the depth.This has profound influence on the convective dynamics.By mass conservation, only a very small fraction of the upwelling fluid from the deep convection zone makes it into the photosphere.The rest must overturn.Over each scale height, the density decreases by a factor of 1/e, so that 1 − 1/e of the mass must overturn.Similarly, the downwelling fluid must entrain mass at this rate.For simple assumptions about the flow geometry, this implies a characteristic horizontal flow scale at each depth d = 4H ρ , where H ρ is the density scale height (Nordlund et al. 2009).Taking this to be the integral (driving) scale of the motions (Stein et al. 2009) allows a simple two component model that can reproduce the observed spectrum of horizontal motions in the solar photosphere (Lord et al. 2014).
For statistically steady motions and small horizontal density gradients compared to the mean vertical stratification, the equation of mass continuity takes the form This suggest two flow regimes.For ∂v z /∂z v z /H ρ the motions are nearly divergenceless and for ∂v z /∂z v z /H ρ the vertical stratification dominates.In these two limits, Equation (9) can be used to determine the horizontal velocity power spectrum given that of the vertical velocity.For small scale motions, high horizontal wavenumbers k h , the motions are nearly isotropic with where v * h • v h is the power spectrum of the horizontal flows and v * z v z is that of the vertical flows.For low wavenumber components of the flow, on the other hand, stratification is important and The cross over between these two regimes occurs at the integral scale 4H ρ at each depth.Using this balance, a model of the horizontal motions observed in the solar photosphere can be constructed from the vertical velocity spectrum at each depth.For example, the highwavenumber vertical-velocity spectrum can be taken to be Kolmogorov above the integral scale 4H ρ , with a mixing-length approximation to determine the total integrated power.Since no scales larger than that are driven at any given depth, the amplitudes of modes with scales larger than the driving scale at any depth can be determined by their decay with height from the depth at which they were last driven.A potential flow approximation on these scales can be used to model that amplitude decrease with height.With these ingredients, working from the bottom of the convection zone upwards, the power spectrum of the horizontal velocity at each depth can be determined.It matches that seen in three-dimensional radiative hydrodynamic simulations (Lord et al. 2014).The broader and critical take away from this simplified approach is that the horizontal-velocity spectrum in the photosphere depends on the vertical-velocity flow amplitudes at depth.The larger the horizontal flow scale observed in the photosphere, the deeper it originates.This suggests that the dramatic decrease in the observed horizontal-velocity power above the supergranular scale on the Sun reflects weak convective driving at depth (at depths below ∼10 Mm) and that supergranulation represents the largest buoyantly driven scale of motion (Lord et al. 2014;Cossette and Rast 2016).
A number of reasons for the low-convective amplitudes at depth are possible, including highly non-local dynamics (maintenance of small scale downflowing plumes generated in the photosphere with little horizontal diffusion) that ensures that the mean stratification of the solar convection zone is closer to adiabatic than numerical models can achieve (Cossette and Rast 2016; Rast 2020), this possibly due to the presence of small scale magnetic field (O'Mara et al. 2016;Bekki et al. 2017); subadiabatic stratification in the deep solar convection zone due to internal heating by radiation; reduced convective amplitudes due to the stabilizing influence of rotation in the lower convection zone (Featherstone and Hindman 2016;Vasil et al. 2021); or a combination of these.Any mechanism that leads to reduced convective amplitudes in the deep layers of the solar convection zone will be reflected in reduced low wavenumber power in the horizontal-flows at the surface.
In radiative hydrodynamic simulations, a factor of ∼2.5 reduction in convective amplitudes below ∼10 Mm depth is sufficient to resolve the convective conundrum (Lord et al. 2014).Though such a reduction has not been yet achieved in a first principles model of convection, we note that, by simple mixing-length scaling arguments (Equation ( 5)), it is consistent with the reduction in superadiabaticity (δ < 2 × 10 −7 ) suggested by the Gizon et al. (2021) analysis of solar inertial modes.

Gyroscopic Pumping and Thermal Wind Balance
Given the convective motions, it is important to understand how, in combination with solar rotation, they generate and maintain large-scale solar differential rotation and meridional circulation.For that purpose, the concepts of gyroscopic pumping and the thermal wind balance are useful (McIntyre 1999;Miesch et al. 2008;Miesch and Hindman 2011).We discuss those in this section.For simplicity, the magnetic field is ignored, though it may be critical in some cases (Hotta et al. 2022).We use a notation where a quantity Q is divided into a mean (longitudinal average) Q and perturbation Gyroscopic pumping is a reflection of angular momentum conservation.The longitudinally averaged longitudinal equation of motion under the anelastic approximation can be written as where L = v φ and v m = v r e r + v θ e θ are the specific angular momentum and the meridional flow velocity, respectively, with = r sin θ in the spherical geometry (r, θ, φ).We also note ρ 0 the density background profile (spherical average).The flow velocity is then divided into components, as above, v = v + v and, assume a steady state ∂/∂t = 0 balance, the equation for gyroscopic pumping can be written While the angular momentum conservation equation ( 12) determines the temporal evolution of the differential rotation, the gyroscopic pumping balance (Equation ( 13)) mainly determines the meridional flow v m in a steady state.The distribution of the specific angular momentum L in the Sun is known to be mostly cylindrical, i.e., ∂ L /∂z ∼ 0, where z denotes the direction of the rotational axis.Thus, gyroscopic pumping can be rewritten as Since we know that the sign of ∂ L /∂ is greater than zero, the sign of the axial torque T = −∇ • ρ 0 v m v φ directly determines the direction of the meridional flow.
The thermal wind balance equation is derived from the longitudinal vorticity equation, where ζ = ∇ × v is the vorticity.This equation is for the evolution of the meridional flow (v r and v θ ) in terms of the longitudinal vorticity ζ φ .In a steady state (∂/∂t = 0), which reduces to the Taylor-Proudman theorem ∂ /∂z = 0 when advection ∇ × (v × ζ ) and the latitudinal entropy gradient ∂ s /∂θ are ignored.Under the Taylor-Proudman constraint, contour lines of the angular velocity must be parallel to the rotational axis.As shown in Fig. 1 the solar differential rotation does not follow this configuration, showing instead the prominent tachocline and the near surface shear layer at its boundaries and a more conical profile in the interior.The latitudinal entropy gradient plays a dominant role in forcing the rotation profile away from the Taylor-Proudman state in the deep convection zone where the rotational influence is important.The Reynolds stresses arising from vigorous surface convection (and the magnetic field which we ignored here) are also crucial in the near surface layer (e.g., Hotta et al. 2022) (see Sect. 2.5 for more details).

Governing Equations for Numerical Simulations
Simulating the global scale motions seen on the Sun directly requires simulating the convective motions in a spherical domain over many scale heights.This in turn requires running efficient numerical code on high-performance supercomputers.For tractability, a number of approximations must be made during formulation.
The solar convection zone is fully ionized below about 20 Mm, and, for global spherical shell magnetohydrodynamic models of convection below that depth, we can reliably assume the equation of state is close to that of a perfect gas.Simulations of the deep solar convection zone must include rotation along with gravitational stratification, and since the superadiabaticity in the solar convection zone is tiny (δ 10 −6 ), solving an entropy equation is preferable to formulations in terms of the total energy or internal energy.The set of equations to be solved can thus be written as Here, we ignore the viscosity, the magnetic diffusivity, and the thermal conductivity due to the large fluid/magnetic Reynolds and Peclet numbers in the Sun.Moreover, since perturbations in the thermodynamic variables scale with the superadiabicity (Equation ( 2)) a linearized equation of state is appropriate, where γ is the specific heat ratio and subscripts 0 and 1 denote the background and perturbed variables, such that Q = Q 0 + Q 1 , and we typically assume the hydrostatic balance for the background stratification ρ 0 , p 0 , and T 0 : Pressure gradient and the gravitational forces in the momentum equation are then due to perturbations about this background state, Due to a large optical depth in the deep convection zone, the diffusion approximation can be used for radiation energy transfer Q rad , where κ rad is the radiation diffusion coefficient estimated from the local opacity.
Employing an anelastic approximation is important for deep solar convection simulations since the speed of sound c s is much faster than the convection velocity v c .Explicitly solving for sound waves in the domain severely restricts the CFL condition on the time stepping t , and a huge number of the time steps would be required to evolve the convective motions and larger-scale flows over dynamical time scales.To avoid this difficulty, the anelastic approximation (Gough 1969) is widely used (Clune et al. 1999), which simplifies the equation of continuity, filtering out sound waves by taking the sound speed to be infinite.In the context of MHD, the anelastic approximation eliminates the fast magneto-sonic waves, while preserving the Alfvén waves and the slow magneto-sonic waves.
Other sonic-filter formulations have been developed including the Lantz-Braginsky-Roberts (LBR) method, in which a reduce pressure is introduced and interactions between fluctuating pressure and stratification are neglected (Lantz 1992;Braginsky and Roberts 1995).The LBR method has the advantage of conserving energy well in both unstable convective zones and stable radiative interiors, critical for simulating gravity wave excitation and propagation (Brown et al. 2012;Vasil et al. 2013).
Recently, a method that has come to be known as the Reduced Speed of Sound Technique (RSST) has also found extensive use.With this, the continuity equation is altered, so that the effective speed of sound is reduced by a factor of ξ (Rempel 2005;Hotta et al. 2012).An advantage of the RSST method is that it does not require global communications in a parallel computing environment, in contrast to the anelastic approximation which requires frequent global communication to solve the elliptic equation ( 20).Thus, the RSST is very useful when solving the MHD equations with massively parallel supercomputers.
Additionally, an inhomogeneous ξ can be employed, taking ξ large in the deep layers where the sound speed is fast and ξ = 1 in the near-surface layer where the anelastic approximation is not valid.This enables simulation over a continuous domain that covers the whole convection zone (Hotta et al. 2019).

Modeling of the Solar Large-Scale Flows
The last fifty years have brought tremendous advances in the simulation of solar and stellar convection and our understanding of the global flows that result in rotating domains.Initial analytic analysis of convection in a rotating sphere (Chandrasekhar 1961;Roberts 1968;Busse 1970) were extended to numerical linear and nonlinear numerical studies aimed at understanding the behavior of rotating turbulent astrophysical bodies (Gilman 1975(Gilman , 1977)).
In particular, the aim was to understand the importance of nonlinear process in both the solar (Gilman 1979) and terrestrial (Cuong and Busse 1981) contexts.In the solar case, this was strongly motivated by the need to explain solar differential rotation, which had at the time been observed for more than one century (Carrington 1860).
The earliest numerical calculations were done using the Boussinesq approximation, but quickly extended to include stratification using the anelastic approximation (Gilman and Glatzmaier 1981;Glatzmaier and Gilman 1982).At the same time, dynamo calculations were solving for the evolution of a magnetic field in these simulations, adding new theoretical constraints on the understanding of its generation within the Sun (Gilman and Miller 1981;Glatzmaier 1984Glatzmaier , 1985)).In parallel, stellar evolution models became sophisticated enough to model in detail the solar stratification (i.e., density, pressure, temperature, equation-of-state, opacity as functions of depth) and these models were found to be highly consistent with deductions based on helioseismic measurements (Model S: Christensen-Dalsgaard et al. 1996).With the stratification well modeled, notable advances could be made in a more realistic reproduction of thermal convection (e.g., Miesch et al. 2000;Brun and Toomre 2002).Although spatial resolutions were moderate (N θ < 128, max < 85), the solar-like differential rotation profile (with faster equator and slower poles) was reproduced in these anelastic studies.
However, the differential-rotation profiles found in numerical solutions tended to obey the Taylor-Proudman theorem, i.e., ∂ /∂z = 0, in contrast to the solar observation (Fig. 1).Motivated by a thorough assessment of mean-field dynamics (Rempel 2005), Miesch et al. (2006) adopted latitudinal entropy gradient at the bottom boundary and achieved non-Taylor-Proudman differential rotation, as indeed is suggested by Equation ( 16).This somewhat ad-hoc method was adopted by several follow-up studies (e.g.Miesch et al. 2008;Fan and Fang 2014).In particular, Brun et al. (2011) included the overshoot layer between the convection zone and the deeper radiative layer below in a self-consistent model, suggesting that the interaction between the two layers is responsible for the crucial latitudinal entropy gradient, the tachocline, and the conical profile observed in the bulk of the convection zone.
This conclusion has not gone without debate.Earlier, Miesch et al. (2000) concluded that anisotropic entropy transport by the overshooting downflows may not be enough to maintain a solar-like differential rotation profile, though Hotta (2018) point out that an efficient small-scale dynamo can amplify the effect (see also Hotta et al. 2015a), and help produce required non-Taylor-Proudman state.Recently Matilsky et al. (2020) have pointed out that the latitudinal entropy gradient achieved is sensitive to the radial boundary condition imposed.Resolution of these uncertainties is critical to our understanding of the origin of the observed radial profile of the solar differential rotation.
The other non-Taylor-Proudman feature in the solar convection zone is the near-surface shear layer (NSSL).The NSSL is thought to be generated by the radially-inward angular momentum transport by near-surface convection which is not strongly influenced by rotation (e.g., Foukal and Jokipii 1975).Numerical simulations have succeeded in reproducing part of the observed feature (Guerrero et al. 2013;Hotta et al. 2015b;Matilsky et al. 2019), with turbulent viscosity playing an important role (Hotta et al. 2015b), however, the models fail to reproduce key aspects of the NSSL, especially at mid-latitude.Matilsky et al. (2019) argue that the role of large-scale columnar convection (banana cells) and meridional circulation and their balance differ with latitude, but are unable to produce a solution that captures the NSSL in both the equatorial and high-latitude regions.Overall, a consensus has not been reached on the maintenance mechanism for the NSSL.

Towards Modeling Other Stars
Using models such as the ones presented previously, numerical simulations are a powerful tool to study the formation and dynamics of large-scale flows in the astrophysical context.In particular, full sphere simulations, resolving a broad range of turbulent convective scales, are suitable tools to probe the large-scale dynamics of distant stars.As examples, numerical studies of solar-type stars recently provided new constraints on the trends we can expect for differential rotation ( ∝ 0.46 * ) on G and K stars (Brun et al. 2022), and highlight the different large-scale flows regimes possible.These regimes can be characterized by the Rossby number (Gastine et al. 2014;Brun and Browning 2017;Hindman et al. 2020).
Thanks to numerical simulations (Matt et al. 2011;Guerrero et al. 2013;Käpylä et al. 2014;Simitev et al. 2015;Brun et al. 2017;Karak et al. 2018), three regimes are currently acknowledged in global rotating models of main-sequence solar-type stars: 1.At low Rossby numbers, differential rotation profiles are highly constrained (Taylor-Proudman) with flows that become cylindrical (Gilman 1975;Elliott et al. 2000;Miesch et al. 2006;Brown et al. 2008).In extreme cases, such profiles show alternating prograde and retrograde zonal jets, also called Jupiter-like jets (Rhines 1975;Heimpel et al. 2016).The meridional circulations under these conditions is typically multicellular in each hemisphere and aligned along the vertical axis (see for example Brun et al. 2017).When magnetic fields are included in the calculation, the Lorentz force feedback can quench the flows, (Brun 2004;Yadav et al. 2015), resulting in a significant decrease of the differential rotation contrast ( -quenching) and the appearance of trans-equatorial meridional-circulation cells (Brun et al. 2022).2. At intermediate Rossby numbers, the differential rotation adopts a typical solar-like conical profile, with a fast equator and slow poles (see for instance Miesch et al. 2006;Hotta et al. 2022).In this regime the power sustaining the differential rotation contrast can reach tens of percent of the stellar luminosity (Brun et al. 2022).3.At high Rossby number (typically over the unity), differential-rotation profiles in simulations become "anti-solar" (Gastine et al. 2014), showing then a slow equator and fast poles (Gilman 1977).
As stars spin down along the main sequence (Skumanich 1972;Gallet and Bouvier 2013;Ahuir et al. 2021), its rotational influence changes and so does by definition its Rossby number which will increase.The resulting transitions between these large-scale flow regimes may then influence the evolution of the star (Metcalfe et al. 2022), and in particular play important roles in global dynamo processes (Strugarek et al. 2018;Warnecke 2018;Guerrero et al. 2019).Large-scale flows indeed strongly influence the nature of such dynamo processes, however the Lorentz force feedback on these flows are likely to play an important role, both in the construction of the solar-type differential rotation profile (Hotta et al. 2022) and magnetic cycle emergence (Gilman 1983;Augustson et al. 2015).Recent studies highlighted that the latter behavior may arise in specific stellar parameter ranges (Strugarek et al. 2018), where the back-reaction of the Maxwell stress modulates the latitudinal differential rotation contrast, and thus the large-scale magnetic field generation.Probing of these modulations (torsional oscillations) have recently been constrained on the Sun (see Sect. 1.1) and quantified for other stars in numerical parametric studies (Brun et al. 2022).
Recent observations show that the composition of the star too plays a role in the convective dynamo (See et al. 2021).One-dimensional stellar models show varying sensitivities to metallicity (Amard and Matt (2020), see also Sect. 4 of Noraz et al. (2022a) for a discussion).Typically, when metallicity is increased at a given stellar mass and age, the opacity also increases, and thus so too does the temperature gradient within the star.The convective zone becomes deeper in proportion to the stellar radius, with longer convective turnover times at its base because of the higher inertia induced by a higher density.This then likely modifies the convective scale distribution and also decrease the Rossby number of the star, which subsequently impacts the dynamics of the large-scale flows (Bessolaz and Brun 2011).Such metallicity effects may lie at the origin of the observed differences in the magnetism of HD 173701 (mentioned in Sect.1.4).Karoff et al. (2018) suggest that the higher metallicity of that star could either enhance the magnetic field generated or the observed facular contrast, yielding then the large amplitude brightness variations observed during the activity cycle.Comparisons between solar twins of different metallicities are needed to either confirm or decipher such mechanisms.To guide those observations, local-area numerical studies of the metalicity impact on photospheric convection have already started (Witzke et al. 2022), and studies using global simulations are currently being investigated (Noraz 2022).
A major point to mention here, is that clear detections of anti-solar rotators are still pending for solar-type stars on the main sequence.In addition to the decrease of sensitivity regarding slow-rotation for current observational techniques, recent results highlight that a change in nature of the dynamo may be induced by a rotational transition toward anti-solar rotation (Karak et al. 2015;Viviani et al. 2019;Noraz et al. 2022b;Brun et al. 2022).In that context, a possible disappearance of star-spots is currently discussed in the community, which would make rotational characterization of the high-Rossby regime even more difficult via photometric techniques.An active search is currently underway to better constrain this regime for solar-type stars (Reiners 2007;Reinhold and Arlt 2015;Benomar et al. 2018;Noraz et al. 2022a), and its implications for the solar case.

Discrepancies Between Observations and Models and Possible Solutions
Significant difficulties remain when making direct comparison of the observed large-scale motions on the Sun and those produced in numerical simulations still exist.Most fundamentally, until recently (Hotta and Kusano 2021), all global spherical shell models of global convection and circulation produced anti-solar differential rotation profiles (a slow equator and a fast pole) at solar rotation rates and at the solar luminosity, and no local area models of solar convection produces a peak in the photospheric horizontal velocity a supergranular scales, as observed.These discrepancies likely reflect a mismatch in convective amplitudes at depth, that has been come to be known as the 'convective conundrum' (O'Mara et al. 2016).The global implications were first recognized when increased numerical resolution in simulations resulted in lower diffusivity, faster convective flows, and consequent difficulty in achieving a solar-like differential rotation profile (Miesch et al. 2008).

The Convection Conundrum
There is an observation mismatch between flow velocities on the Sun and those found in numerical simulations.In the surface layers this is dramatically illustrated by the absence of supergranulation in local area models, compared to its importance in photospheric observations.That aspect of the convective conundrum is discussed in Sect.1.2 and 2.2 above.In summary, these disparities suggest reduced convective amplitudes with depth, by a factor of about 2.5 below 10 Mm, with several causes for this reduction suggested in the literature, any one or more of which would suffice.
Convective amplitudes in the near solar surface region of the solar convection zone can be measured using a variety of local helioseismic techniques.As illustrated by Fig. 5 these do not agree with each other or with numerical simulations (Hanasoge et al. 2012;Gizon and Birch 2012;Greer et al. 2015).In general models produce flow with significant power at low wavenumbers, monotonically increasing to low wavenumber, while observations indicate reduced power there.A modeling exception is the recent Hotta and Kusano (2021) simulation, in which the horizontal velocity power rolls over at large scales, though at scales somewhat beyond that of supergranulation.An observational exception is the ring-diagram result of Greer et al. (2015) which shows significant power at scales larger than supergranulation at a depth of 0.96 R .
While refinement and revisions of the techniques employed have brought measurements and models closer together (Nagashima et al. 2020;Proxauf 2020), there are significant and important remaining discrepancies.It is imperative to resolve these, not just to understand the Sun and its dynamo, but because the solar case serves as the touchstone for stellar modeling.

Columnar Convective Modes
Numerical simulations of rotating convection have repeatedly found that the giant-cell convection tends to exist as banana cells in a strongly rotationally-constrained regime (e.g., Miesch et al. 2000).These banana cells are located outside the tangential cylinder and can be seen as north-south aligned downflow lanes across the equator.They are known to propagate in a prograde direction with frequencies higher than the local differential rotation rate (Miesch et al. 2008).
The prograde propagation of banana cells can be understood in terms of a special class of inertial modes called columnar convective modes or thermal Rossby modes (e.g., Busse 2002;Miesch et al. 2008;Bekki et al. 2022a).They are z-vorticity waves arising from the compressional β-effect due to the strong background density stratification.A linear dispersion relation of the columnar convective modes was first derived by Glatzmaier and Gilman (1981) using a cylindrical model and later by Hindman and Jain (2022) in Cartesian geometry.The most realistic spherical-shell model has recently been presented by Bekki et al. (2022b).They have shown that the dispersion relation of the columnar convective modes is very sensitive to the superadiabaticity δ, and that the modes become convectively-unstable (exponentially growing) when the background is slightly superadiabatic (δ > 0).These modes have the associated Reynolds stress v r v φ > 0 and v θ v φ > 0 (< 0) in the northern (southern) hemisphere, implying that they transport the angular momentum radially upward and equatorward outside the tangential cylinder.The great importance of banana cells on the establishment of the solar-like differential rotation has been repeatedly appreciated in previous literature (Käpylä et al. 2011;Gastine et al. 2013;Hotta et al. 2015b;Matilsky et al. 2020;Camisassa and Featherstone 2022).
Despite their prominence in numerical simulations, columnar convective modes have not been observed on the surface of the Sun, though very large-scale flows have been deduced from correlation tracking of the solar supergranulation (Hathaway et al. 2013;Hathaway and Upton 2021) and possible indirect evidence may come from the alignment of the solar supergranulation (Lisle et al. 2004;Nagashima et al. 2011).It remains a mystery why the columnar convective modes are not directly detected in the solar surface observations.One possible scenario is that they indeed exist hidden in the deep convection zone with substantial amplitudes but are concealed by surface small-scale convective feature (Guerrero et al. 2013).The other scenario is that they are simply absent in the Sun or too weak to be detected (van Ballegooijen 1986;Lord et al. 2014, Sect. 2.2 above).If the latter is the case, the angular momentum needs to be transported by something other than the large-scale columnar convective modes.For instance, the equatorward angular momentum transport can be achieved in the Sun by the Reynolds and/or Maxwell stresses at much smaller spatial scales.

Differential Rotation
An important aspect of the convective conundrum is its implications for differential rotation.In early phase of the solar differential rotation research, the solar differential rotation profile could be reproduced relatively easily.As supercomputer power grew, allowing simulation of higher resolution, the problem became apparent.Simulations began to fail to reproduce the solar-like differential rotation profile at solar rotation rates.
A well-known feature of differential rotation is that a fast equator (poles) is obtained with strong (weak) rotational influence (e.g., Gastine et al. 2013;Featherstone and Miesch 2015;Karak et al. 2015).The rotational influence is measured by the Rossby number Ro = v/(2 H ), where v is a characteristic convective velocity and H a characteristic spatial scale.Although we believed that the Sun is in a low Rossby-number regime, high resolution simulations most often produce anti-solar differential rotation profiles (fast poles).High resolution simulations introduce the small-scale turbulence which is important for heat transport but is only very weakly rotationally constrained.This has been recently confirmed with scale-dependence analyses of the angular momentum flux (Mori and Hotta 2023).The small-scale turbulence tends to transport the angular momentum radially inward which leads one-cell meridional flow, which transports angular momentum and accelerate the poles (see also Featherstone and Miesch 2015).In other words, high-resolution decreases the convective spatial scale H placing the simulation in a high Rossby number regime, resulting an anti-solar differential rotation.
Three mechanisms have been employed in global simulations in order to maintain a solar-like differential rotation profile, i.e, to maintain the low Rossby number, in the face of vigorous smaller-scale convective flows.
• Reducing the luminosity L (Hotta et al. 2015b), which leads to a reduction of v.
• Increasing viscosity ν and/or thermal conductivity κ (Miesch et al. 2000(Miesch et al. , 2008;;Fan and Fang 2014;Hotta et al. 2016), which leads to a reduction of v and an increase of H . • Increasing radiative conductivity κ rad (Käpylä et al. 2019;Noraz 2022), which leads to a reduction of v while conserving L.
None of these reflect a possible physical mechanism operating on the Sun but not captured in simulations.They cannot provide an answer to the question: Why does the Sun have such a low Rossby number?Recently, Hotta and Kusano (2021) reproduced the solar-like differential rotation in extremely high resolution simulation without using these manipulations.Hotta et al. (2022) showed that in that simulation the strong magnetic field maintained by the convection when the magnetic diffusivity is very low, as it is in their simulation, causes angular momentum to be transported outward.A double-cell meridional flow is generated and that flow transports angular momentum equatorward, resulting in the solar-like differential rotation profile achieved.An important point is that the Hotta and Kusano (2021) simulation convection remains in a high Rossby number regime but the columnar convective modes no longer play an dominant role in the angular momentum transport.The Rossby number range in which a solar-like differential rotation regime is possible is under investigation, and recent studies suggest that this range depends sensitively on the Prandtl number (Käpylä 2023;Noraz 2022).
4 Future Prospects a theoretical view point, some more significant progress can be made in the near future.Hotta and Kusano (2021) show that the high-resolution simulations remain a promising approach to understand the solar large-scale flows.Simulations have not yet reached numerical convergence, where by numerical convergence we mean that when we double the resolution the large-scale structure does not change.The solar turbulent convection in the simulations has a spatial spectrum which has more power at low wavenumbers than is observed on the Sun.With the injection scale at 200 Mm and the dissipation scale around 1 cm, it is impossible to directly resolve all the energy containing scales with numerical simulations.Although higher-resolution numerical simulations are required to better understand the nature of turbulent convection, it is also crucial to construct reliable turbulence models which can mimic the essential features of unresolved flows.
Such models require a deep understanding of the key physical components of the convection.Important progress has been made in understanding the essential role of rotation (e.g., Vasil et al. 2021, and references therein) and preliminary work has begun to characterize the highly nonlocal convective flows that result from radiative cooling of the photosphere which may play an important role in heat transport and allow a mean gradient much closer to isentropic than that achieved by current simulations (Brandenburg 2016;Cossette and Rast 2016).Moreover, radiative heating of the lower convection zone likely plays an important role, one that has only begun to be examined.Brun et al. (2011) have undertaken an important study of the interaction between the convection and radiation zone, but, as Käpylä (2019) pointed out, a fully consistent treatment of the overshoot layer requires huge numerical resources, which are currently not available (see also Hotta 2017).Solving these problems on convection is essential to understanding the large-scale flow dynamics.
Another direct approach is to combine in one simulation the radiative magnetohydrodynamics of the photosphere and the global scale convective motions.Including the photosphere may have a significant impact on the deep large-scale convection (Spruit 1997).Hotta et al. (2019) found only weak (or no) influence from the photosphere on the deep convection, however their result may not reflect a lack of importance of these flows, but the difficulties faced in maintaining them with depth given the diffusivities required in global models.Even when the resolution of the simulation is increased (and numerical or explicit diffusivities are reduced), the problem persists because the horizontal scales of the flow structures are also decreased.What is required to maintain nonlocal transport in a simulation is that the diffusion time across a downflowing plume be greater than the transit time of that fluid across the fluid layer.This is extremely difficult to achieve.Because of the steep stratification, local area radiative magnetohydrodynamic simulations, and ambitious global models which include an upper radiative boundary, may be able to produce granular scale downflow structures but not the low thermal diffusivity required for them to maintain their role in transport at depth.
Over the past fifty years, tremendous progress has been made in understanding and simulating solar convection and the global scale flows that result in a rotating domain.The Sun allows us to directly confront that progress with observations that, with the advent and development of helioseismology over the same time period, have also made previously unimaginable advances.New diagnostics based on the study of the solar inertial modes are expected to provide additional constraints on the physics of the deep convection zone.Some gaps in the observations and some fundamental issues in our understanding however remain.Resolving those will not only advance our understanding of the origin of large-scale global flows, but will also allow us to more robustly model the solar dynamo, perhaps predict solar behavior critical to human activities in space and on Earth, and extend our understanding to other stars for which comparison with data will remain less constraining.Understanding the Sun is thus a touchstone activity.

Fig. 1
Fig. 1 Internal profile of the solar differential rotation deduced from global helioseismology and averaged from April 2010 to February 2021.Panels (a) and (b) show the results obtained by the Global Oscillation Network Group (GONG) (data courtesy of R. Howe, using the method of Howe et al. 2005) and the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO) (Larson and Schou 2018).Panel (c) shows the radial differential rotation at selected latitudes.Grey shades denote the layers of strong radial rotational shear known as the tachocline and the near-surface shear layer (NSSL)

Fig. 2
Fig. 2 Latitudinal component of the meridional flow inferred by time-distance local helioseismology.The red and blue shades correspond to the northward and southward directions respectively.(a) The result obtained by Zhao et al. (2013) using SDO/HMI data (2010-2012), reproduced with permission.(b) The result obtained by Gizon et al. (2020a) using GONG data (2008-2019), reprinted with permission from AAAS

Fig. 3
Fig. 3 From Hathaway et al. (2015): Solar Doppler velocity spectrum as determined from Helioseismic and Magnetic Imager (HMI) observations (red curves, with and without removal of image artifacts), along with a random-phase synthetic spectrum (black curve).Vertical blue dashed fiducial lines have been added indicating the approximate scale of supergranulation and giant-cells.The blue dotted line approximately overplots the spectrum seen in numerical simulations.Figure without added dashed and dotted blue lines used with permission of Hathaway et al. (2015)

Fig. 4
Fig. 4 Observational power spectra in the Carrington frame and eigenfunctions of three selected inertial modes of the Sun.Top row: Power spectra of the longitudinal component of velocity u φ for m = 1 (left column) and m = 2 (middle column), and power spectrum of the colatitudinal component of velocity u θ for m = 3.The blue curves show the differential rotation rate at r = 0.96 R and R .The purple contour indicates the region affected by active-region flows.Middle row:The same power spectra but normalized at each latitude by their average value over the frequency range between the orange bars.Excess power is seen at a specific frequency at all latitudes in each of the three cases.The red arrows point to critical latitudes at the surface for the case m = 2. Bottom row: Observed horizontal velocity eigenfunctions at the surface for the m = 1 high-latitude mode, the m = 2 critical-latitude mode, and the m = 3 equatorial Rossby mode.These figures are courtesy ofGizon et al. (2021) Figure 1 ofLord et al. (2014) confirms that this two-component continuity balance determines the relationship between the vertical and horizontal-velocity power spectra in radiative hydrodynamic simulations.

Fig. 5
Fig. 5 Horizontal velocity power spectra near the solar surface.(a) Comparison of the spectra between numerical simulations and the observations at r = 0.96 R (Birch 2023).Blue: A global full-spherical simulation of rotating magneto-convection by Hotta and Kusano (2021).Navy: A local cartesian box simulation of solar convection by Hotta et al. (2019).Black (gray area): Observational upper limits inferred by deepfocusing time-distance helioseismic measurement of Hanasoge et al. (2012), revised recently by Proxauf (2020).(b) Comparison of the spectra obtained by various observational measurements at various depths.Red: Multi-ridge fitting ring-diagram analysis by Greer et al. (2015), revised recently by Nagashima et al.(2020).Green: Local correlation tracking of surface granulation(Proxauf 2020).Magenta: The SDO/HMI ring-diagram pipeline(Bogart et al. 2011a,b;Proxauf 2020).All observations reported in the above plots are available online(Birch 2023)