Issue 
A&A
Volume 652, August 2021



Article Number  L6  
Number of page(s)  22  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/202141462  
Published online  06 August 2021 
Letter to the Editor
Solar inertial modes: Observations, identification, and diagnostic promise^{⋆}
^{1}
MaxPlanckInstitut für Sonnensystemforschung, 37077 Göttingen, Germany
email: gizon@mps.mpg.de
^{2}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, 37077 Göttingen, Germany
^{3}
Center for Space Science, NYUAD Institute, New York University Abu Dhabi, Abu Dhabi, UAE
^{4}
W. W. Hansen Experimental Physics Laboratory, Stanford University, Stanford, CA 94305, USA
^{5}
AIM, CEA, CNRS, Universités Paris et ParisSaclay, 91191 GifsurYvette Cedex, France
^{6}
Institut Supérieur de l’Aéronautique et de l’Espace (ISAESUPAERO), 31400 Toulouse, France
^{7}
National Solar Observatory, Boulder, CO 80303, USA
Received:
3
June
2021
Accepted:
1
July
2021
The oscillations of a slowly rotating star have long been classified into spheroidal and toroidal modes. The spheroidal modes include the wellknown 5min acoustic modes used in helioseismology. Here we report observations of the Sun’s toroidal modes, for which the restoring force is the Coriolis force and whose periods are on the order of the solar rotation period. By comparing the observations with the normal modes of a differentially rotating spherical shell, we are able to identify many of the observed modes. These are the highlatitude inertial modes, the criticallatitude inertial modes, and the equatorial Rossby modes. In the model, the highlatitude and criticallatitude modes have maximum kinetic energy density at the base of the convection zone, and the highlatitude modes are baroclinically unstable due to the latitudinal entropy gradient. As a first application of inertialmode helioseismology, we constrain the superadiabaticity and the turbulent viscosity in the deep convection zone.
Key words: Sun: rotation / Sun: oscillations / Sun: interior / Sun: helioseismology / Sun: general
Movie associated to Fig. 2 is available at https://www.aanda.org
© L. Gizon et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
The free oscillations of a nonrotating spherical star have zero radial vorticity and are called spheroidal modes: they are the pressure (p), surfacegravity (f), and gravity (g) modes. The p and f modes, discovered on the Sun by Leighton et al. (1962), are used to infer the structure and dynamics of the solar interior (ChristensenDalsgaard 2002). The solar g modes would also have important diagnostic potential regarding the radiative interior of the Sun; however, they evanesce in the convection zone and their amplitudes at the surface are exceedingly small (García et al. 2007; Alvan et al. 2015).
When slow uniform rotation is included in the model, additional modes of oscillation become possible. In particular, quasitoroidal modes that resemble classical Rossby modes, known as r modes, are predicted (Papaloizou & Pringle 1978). They owe their existence to the Coriolis force, have frequencies on the order of the rotation frequency, and propagate in the retrograde direction. Adding the Sun’s differential rotation introduces critical latitudes where the phase speed of a mode is equal to the local rotation velocity. In the inviscid case, the eigenvalue problem is singular at the critical latitudes (Watson 1981; Charbonneau et al. 1999). Adding viscosity changes the eigenvalue problem from second order to fourth order (e.g., Baruteau & Rieutord 2013). The singularity disappears and new quasitoroidal modes appear, which are analogous to those of the plane Poiseuille viscous flow in classical hydrodynamics (Gizon et al. 2020a, and references therein). In the following, we loosely refer to the modes with frequencies on the order of the rotational frequency as inertial modes.
Inertial modes were detected on some rapidly rotating stars (see the review by Aerts 2021). The search for the Sun’s inertial modes requires observations over many times the 27day solar rotation period due to their low frequencies and amplitudes. Equatorial Rossby modes modified by the solar differential rotation have already been reported (Löptien et al. 2018). Here we report observations of a rich spectrum of inertial modes of the Sun over a wide range of latitudes, and we show they can be used to directly probe the superadiabaticity and turbulent viscosity in the deep convection zone. The degree to which the lower half of the convection zone is superadiabatic (or subadiabatic) is important in the context of storing the toroidal magnetic field so that it can build up over the course of the 11year solar cycle (Hotta 2017). The turbulent viscosity is one of the important turbulent transport processes that acts in combination with the observed meridional flow (Gizon et al. 2020b) to explain the equatorward drift of the latitudes at which sunspots emerge (Cameron & Schüssler 2016).
By definition, a normal mode is separable in time and space; it is characterized by a single eigenfrequency that is independent of position and by a displacement eigenfunction that is independent of time. Working in the frequency–latitude domain is key to the observational discovery and the identification of the quasitoroidal normal modes of the Sun.
2. Observations
We use helioseismic maps of horizontal flows near the solar surface provided by the Stanford ringdiagram pipeline (Bogart et al. 2011a,b) applied to continuous highresolution observations from the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO) for the period from 1 May 2010 to 6 September 2020. The two horizontal flow components are standard data products: u_{θ}(θ, ϕ, t) in the colatitudinal direction and u_{ϕ}(θ, ϕ, t) in the longitudinal direction (θ and ϕ increase southward and prograde, respectively; see Proxauf et al. 2020). The flows are measured either with a cadence of dt = 27.28 h and an effective spatial resolution of 15° in both coordinates, or with a more rapid cadence of dt/3 and a finer spatial resolution of 5°; the spatial sampling is half the resolution such that there is a 50% overlap between neighboring measurements. The highest latitude is 67.5° for the lowresolution maps and 80.0° for the highresolution maps. The longitude, ϕ, is defined in the Carrington frame of reference, which rotates at the frequency Ω_{Carr}/2π = 456.0 nHz with respect to an inertial frame (close to the equatorial rotation rate at the surface). The zero and yearly frequencies were removed from the data.
The structure of the Sun and its differential rotation is nearly symmetric with respect to the solar equator. Consequently, modes with a toroidal component can be called either symmetric or antisymmetric depending on the north–south symmetry of the surface radial vorticity. This terminology has been used before in the literature (Charbonneau et al. 1999). A symmetric mode has a symmetric u_{θ} and antisymmetric u_{ϕ}, while an antisymmetric mode has an antisymmetric u_{θ} and symmetric u_{ϕ}. After symmetrizing (superscript “+”) or antisymmetrizing (“−”) the data with respect to the equator, we computed the Fourier transform of the two velocity components in longitude and in time,
where j is either θ or ϕ, ω is the angular frequency, m is the integer longitudinal wavenumber, and the sums were taken over all longitudes and all times. We considered frequencies in the range ω/2π ≤ 400 nHz and m in the range from 1 to 10 to focus on the largescale motions. For each choice of velocity component j, symmetry s, and wavenumber m, the power spectral density is a function of colatitude and frequency.
For illustration purposes, we show the detection of three global modes of oscillation in the inertial frequency range in Fig. 1 (15° resolution). For each mode, there is clear excess power at the same frequency over a range of latitudes. Several types of modes can be seen. The symmetric m = 1 mode at a frequency near −86 nHz is visible at all latitudes in the power spectrum; it has most of its power at latitudes of 50° and above (the 5° observations show that the power keeps increasing with latitude up to at least 80°). It corresponds to the highlatitude velocity features previously reported (Hathaway et al. 2013; Bogart et al. 2015; Hathaway & Upton 2021), although it was not recognized as a normal mode of the whole convection zone. The second example is the symmetric m = 2 mode of oscillation at −73 nHz. This mode is also seen over the entire latitude range of the observations, but it has most of its power concentrated near the critical latitude of 38° (see Fig. 1). The power is strong above the critical latitude, but decreases toward the poles. The third example is the m = 3 equatorial Rossby modes (Löptien et al. 2018) at a frequency of −269 nHz, for which the power is mostly confined to lie between the critical latitudes (±59° for this mode’s frequency) where the mode is trapped (Gizon et al. 2020a).
Fig. 1.
Power spectra showing selected modes of oscillation in the Carrington frame. Each column corresponds to a particular m and velocity component, as indicated at the top. Each row shows a different representation of the power spectrum. Top row: the power spectral density is plotted as a function of frequency and latitude. The two blue curves show m(Ω − Ω_{Carr})/2π at the surface and at r = 0.95 R_{⊙}, where Ω(r, θ) is the solar angular velocity in the inertial frame. The purple contour delineates the region in frequency–latitude space affected by inflows into active regions, m(Ω_{AR} − Ω_{Carr})/2π (see Fig. C.1). Second row: the power at each latitude is normalized by its average value over the frequency range between the red bars; this shows that each mode has excess power over a large range of latitudes. The red arrows point to the critical latitudes of ±38° at the surface for the mode with frequency −73 nHz. Third row: the power is averaged over the selected latitude bands specified on the plots, and the frequency resolution is reduced to 12.24 nHz. The red dots point to modes that are not activityrelated (see text) and they are listed in Table A.1. 
We have detected many tens of normal modes of oscillation at low frequencies, as shown in Figs. C.2–C.11 and reported in Table A.1. These modes are associated with significant (above 95% confidence level) excess power in at least one of three latitude bands (low latitudes below 30°, mid latitudes from 15° to 45°, and high latitudes from 37.5° to 67.5°). While the most striking features in the power spectra are narrow peaks, a closer inspection reveals ranges in frequency and latitude of additional excess power. For example, the m = 8 power spectrum for (Fig. C.13) has excess power at low latitudes at frequencies between −135 nHz and −65 nHz, which can be attributed to the presence of a dense spectrum of modes adjacent to the equatorial Rossby mode.
To avoid misidentifying activeregion inflows (Gizon et al. 2001) as modes of oscillation in the power spectra, we defined a region in frequency–latitude space based on the activeregion rotation rates and latitudes (Kutsenko 2021), as shown in Fig. C.1. A peak in the power spectrum for the entire observation period is not reported in Table A.1 if it is in the activity area and does not have significant power during the quietSun period (February 2018 to September 2020), see Figs. C.2–C.11. We also checked that the reported modes were not misidentified due to leakage from the window function (Liang et al. 2019), and that they are also seen in the Global Oscillation Network Group (GONG) data (Fig. C.12).
For each mode, we extracted the two velocity components of a mode eigenfunction in a narrow frequency range around the mode frequency within one linewidth (Proxauf et al. 2020). Examples of the surface eigenfunctions are shown in the left column of Fig. 2.
Fig. 2.
Observed and model eigenfunctions for the modes shown in Fig. 1. Left column: observed velocity ( for the m = 1 and m = 2 modes, for the m = 3 mode). Middle columns: corresponding eigenfunctions of the 2D model for ν_{t} = 100 km^{2} s^{−1} and δ = 0, at the surface and through the central meridian, together with the kinetic energy density. The thick black curves show the critical latitudes. Rightmost column: eigenfunctions of the 1D model at the surface. The retrograde propagation of these modes in the Carrington frame is illustrated as an online movie. The other velocity components are shown in Fig. C.15 and the radial vorticity is shown in Fig. C.16. 
3. Mode identification
To identify the observed modes of oscillation, we computed the eigenmodes of a spherical shell with 0.710 ≤ r/R_{⊙} ≤ 0.985 that rotates like the Sun. The internal rotation rate is specified by pmode helioseismology (Larson & Schou 2018). For the sake of simplicity, we chose the model to have very few free parameters: a constant fluid viscosity ν_{t} and a constant superadiabaticity δ, which gives the degree of convective instability. Each mode eigenfunction is proportional to exp(imϕ − iσt), where σ is the complex eigenfrequency. For each m, we solved the 2D (r − θ) eigenvalue problem (Appendix B.1). In order to highlight the main physics, we also computed the purely toroidal (horizontal) modes at fixed radius r = R_{⊙}. For each m, we solved the 1D (θ) eigenvalue problem where the only free parameter is the turbulent viscosity (Appendix B.2). The oscillation spectrum of this 1D model is much less cluttered (no radial overtones for the inertial modes and no convective modes).
We then sought a match with the observed modes. We did not tune the parameters of the 2D model to match the observations exactly; we performed a sensitivity study using δ = −10^{−6}, −2 × 10^{−7}, 0, 2 × 10^{−7}, 10^{−6} (Fig. C.17) and ν_{t} = 50, 100, 250, 500 km^{2} s^{−1} (Fig. C.18). We found that δ = 0 and ν_{t} = 100 km^{2} s^{−1} provide a good match (Fig. 2) for the surface eigenfunctions and eigenfrequencies of the three modes of Fig. 1. As part of the identification, we sought modes of the models that have long lifetimes or are growing (Fig. C.14). The identified modes are representatives of three main families of modes: the highlatitude inertial modes (Fig. 2a), the criticallatitude inertial modes (Fig. 2b), and the equatorial Rossby modes (Fig. 2c). This classification is supported by the dispersion relations at small m (Fig. 3).
Fig. 3.
Mode frequencies in the Carrington frame for the observations and the 2D model (Re σ). The symbols show the observed modes (diamonds for symmetric modes and squares for antisymmetric modes). The red symbols show the highlatitude modes, the orange symbols the criticallatitude modes, and the black symbols the equatorial Rossby modes. The rose and grayshaded areas show the observed frequency ranges of excess power (last column of Table A.1). For reference, the blueshaded area gives the range of rotation rates at the equator between the surface and 0.95 R_{⊙}. The curves give the dispersion relations for the modes of the 2D model with ν_{t} = 100 km^{2} s^{−1} and δ = 0. The red curve is the dispersion relation for the highlatitude modes. The solid and dashedblack curves are for the fundamental (n = 0) and first overtone (n = 1) equatorial Rossby modes. 
The highlatitude inertial modes are analogous to the “wall modes” in plane Poiseuille flows (Gizon et al. 2020a). They are seen in both the 1D and 2D eigenvalue problems for m ≤ 5. In the 2D model, the eigenfunctions are dominantly toroidal and extend to the bottom of the convection zone, with their highest kinetic energy density near the base of the convection zone (Fig. 2a). This is unlike the kinetic energy density of the p modes, which always peaks near the surface. The correct tilt of the spiral structure is only obtained in the 2D model (Fig. 2a). In this model, the highlatitude modes become baroclinically unstable (Fig. C.14) due to the latitudinal entropy gradient resulting from the thermal wind balance (Knobloch & Spruit 1982; Bekki et al., in prep.)^{1}. In the 1D model, only the m = 1 highlatitude modes are selfexcited (Fig. C.14) as a result of a shear instability at high latitudes; however, the tilt of the spiral is not consistent with the observations.
Criticallatitude inertial modes are found for both the 1D and 2D models. Their amplitudes are maximum near their critical latitudes; they are known as “center modes” in 1D hydrodynamics (Gizon et al. 2020a). The kinetic energy density of the m = 2 mode of the 2D model at −92 nHz (−73 nHz observed) is concentrated near the base of the convection zone near 45° latitude (Fig. 2b). This is a very important place in the Sun, as it is where the toroidal magnetic field generation should be strongest (Spruit 2011).
Equatorial Rossby modes are the easiest to identify: their frequencies are close to the classical dispersion relation for uniform rotation, σ = −2Ω/(m + 1). The 2D model supports modes with a different number of nodes in the radial direction, n. The frequencies of the observed equatorial Rossby modes span the range between the n = 0 and the n = 1 branches of the dispersion relation (Fig. 3). For example, the m = 3 mode is identified as a fundamental mode (n = 0).
4. Conclusion
We observed and identified three families of globalscale inertial modes in the solar convection zone, within the search range ω/2π ≤ 400 nHz and 1 ≤ m ≤ 10. Some of these modes are selfexcited in the models. We also found extended regions in frequency space where closely packed modes exist. The modes we have identified are sensitive to the physical conditions deep in the convection zone (see plots of the kinetic energy density in Fig. 2). The eigenfrequencies and surface eigenfunctions of the high and criticallatitude inertial modes have diagnostic potential for the latitudinal entropy gradient, the superadiabaticity (Fig. C.17), and the turbulent viscosity (Fig. C.18), which are largely unconstrained by traditional pmode helioseismology. We find that the observed inertial modes are compatible with δ < 2 × 10^{−7} and ν_{t} ≤ 100 km^{2} s^{−1} at the bottom of the convection zone. These observational upper limits are substantially below the expectation from mixing length theory – by approximately one order of magnitude each (ChristensenDalsgaard et al. 1996; MuñozJaramillo et al. 2011) – and they imply that the convective motions in the lower half of the convection zone are weak. This might correspond to the slightly subadiabatic conditions seen below 0.8 R_{⊙} in recent numerical simulations of solar convection (see e.g. Hotta 2017; Käpylä et al. 2017; Bekki et al. 2017). While our upper limit on the turbulent velocities ( m s^{−1} for a correlation time τ = 1 month) is well below the mixing length value, it is just above the lower limit required to drive solar differential rotation (8 m s^{−1} according to Miesch et al. 2012). A lower convection zone that is only marginally unstable (or even stable) would allow a flux transport dynamo to wind up and transport the magnetic field in this region. We expect that the characteristics of the observed inertial modes, including amplitudes and lifetimes, will allow us to infer δ(r) and ν_{t}(r) and understand in which regime of rotating convection the Sun operates (Hindman et al. 2020).
The formation of a spiral at high latitudes by baroclinic instability has also been discussed in the context of Venus’ atmosphere (Kashimura et al. 2019).
Acknowledgments
Author contributions: this project was initiated and supervised by LG and RHC. The helioseismic ring parameter fits were provided by RSB for HMI and KJ for GONG. BP and ZCL measured the mode parameters. YB solved the 2D eigenvalue problem. DF and LH solved the 1D eigenvalue problem. LG, RHC and ACB wrote the draft paper. All authors contributed to the final manuscript. We are very grateful to John Leibacher for useful comments. LG thanks the Max Planck Institute for Astrophysics for the opportunity to present a preliminary account of these results at the 2020 Biermann Lectures. The HMI data are courtesy of NASA/SDO and the HMI Science Team. This work utilizes GONG data from the National Solar Observatory (NSO), which is operated by AURA under a cooperative agreement with NSF and with additional financial support from NOAA, NASA, and USAF. This work was supported in part by NASA contract NAS502139 to Stanford University. YB is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen, and acknowledges partial support from the Japan Student Services Organization (JASSO). We acknowledge support from ERC Synergy Grant WHOLE SUN 810218. LG acknowledges support from NYUAD Institute Grant G1502. LG, DF and BP acknowledge funding by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through SFB 1456/432680300 Mathematics of Experiment, project C04. The computational resources were provided by the German Data Center for SDO through German Aerospace Center (DLR) grant 50OL1701. LG, ACB, and CD acknowledge support from DLR under PLATO Data Center grant 50OO1501.
References
 Aerts, C. 2021, Rev. Mod. Phys., 93, 015001 [Google Scholar]
 Alvan, L., Strugarek, A., Brun, A. S., Mathis, S., & Garcia, R. A. 2015, A&A, 581, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baruteau, C., & Rieutord, M. 2013, J. Fluid. Mech., 719, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & RabelloSoares, M. C. 2011a, J. Phys. Conf. Ser., 271, 012008 [Google Scholar]
 Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & RabelloSoares, M. C. 2011b, J. Phys. Conf. Ser., 271, 012009 [Google Scholar]
 Bogart, R. S., Baldner, C. S., & Basu, S. 2015, ApJ, 807, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R. H., & Schüssler, M. 2016, A&A, 591, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charbonneau, P., Dikpati, M., & Gilman, P. A. 1999, ApJ, 526, 523 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2002, Rev. Mod. Phys., 74, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Däppen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [CrossRef] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [Google Scholar]
 Gizon, L., Duvall, T. L., Jr., & Larsen, R. M. 2001, in Recent Insights into the Physics of the Sun and Heliosphere: Highlights from SOHO and Other Space Missions, eds. P. Brekke, B. Fleck, & J. B. Gurman, 203, 189 [Google Scholar]
 Gizon, L., Fournier, D., & Albekioni, M. 2020a, A&A, 642, A178 [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020b, Science, 368, 1469 [Google Scholar]
 Hathaway, D. H., & Upton, L. A. 2021, ApJ, 908, 160 [CrossRef] [Google Scholar]
 Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Hindman, B. W., Featherstone, N. A., & Julien, K. 2020, ApJ, 898, 120 [CrossRef] [Google Scholar]
 Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
 Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Kashimura, H., Sugimoto, N., Takagi, M., et al. 2019, Nat. Commun., 10, 23 [CrossRef] [Google Scholar]
 Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
 Kutsenko, A. S. 2021, MNRAS, 500, 5159 [CrossRef] [Google Scholar]
 Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29 [Google Scholar]
 Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ, 135, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Liang, Z.C., Gizon, L., Birch, A. C., & Duvall, T. L., Jr. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
 Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618 [Google Scholar]
 Miesch, M. S., Featherstone, N. A., Rempel, M., & Trampedach, R. 2012, ApJ, 757, 128 [NASA ADS] [CrossRef] [Google Scholar]
 MuñozJaramillo, A., Nandy, D., & Martens, P. C. H. 2011, ApJ, 727, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 182, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Proxauf, B., Gizon, L., Löptien, B., et al. 2020, A&A, 634, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spruit, H. C. 2011, in The Sun, the Solar Wind, and the Heliosphere, eds. M. P. Miralles, & J. Sánchez Almeida, 4, 39 [Google Scholar]
 Watson, M. 1981, Geophys. Astrophys. Fluid Dyn., 16, 285 [CrossRef] [Google Scholar]
Appendix A: Parameters of observed modes
Solar inertial modes detected in HMI ringdiagram flow maps for 2010–2020.
Appendix B: Normal modes of the differentially rotating Sun
B.1. 2D eigenvalue solver
In the Carrington frame, the linearized equations for the conservation of momentum, mass, and energy, together with the equation of state, are as follows:
where
is the material derivative and Ω(r, θ) is a differential rotation model close to the helioseismic measurements averaged over 2010−2020 (Larson & Schou 2018). Linear perturbations are denoted with primes. The background model is based on a standard solar model (ChristensenDalsgaard et al. 1996), except for the superadiabaticity δ = ∇ − ∇_{ad}, which is a constant parameter in the convection zone (the radiative zone is very stable below with δ ≈ −0.1). In the above equation, H_{p} is the pressure scale height, and c_{v} and c_{p} are the heat capacities per unit mass at constant volume and constant pressure. The viscous stress tensor accounts for wave attenuation, where δ_{ij} is the Kronecker delta. The energy equation includes advection and thermal diffusion. In our model, the viscous and thermal diffusivities are those resulting from the turbulence and, therefore, were considered to be equal.
The latitudinal entropy gradient is obtained by assuming that the differential rotation is the result of a thermal wind balance (Miesch et al. 2006):
where z = r cos θ is the coordinate along the rotation axis.
Boundary conditions need to be applied at θ = 0 and π. Since we only considered modes with m ≠ 0, we imposed u′ = 0 and ρ′=p′=s′ = 0. The numerical domain is bounded above by the photosphere and below by the radiative interior, both of which are strongly stably stratified and where radial flows are difficult to drive because of the strong, restoring buoyancy force. Therefore, we used an impenetrable and stressfree boundary condition at both radial boundaries.
We looked for solutions to the above problem where each physical quantity is proportional to exp(imϕ − iσt), where σ is the complex mode angular frequency and m is the integer longitudinal wavenumber. We discretized the spatial derivatives with secondorder central differences with 16 radial and 72 latitudinal grid points. The above equations were combined in matrix form into a complex eigenvalue problem, which was solved using the LAPACK routine. We focused on the lowfrequency solutions. We refer to the modes of oscillation obtained in this problem as the “modes of the 2D model”.
B.2. 1D eigenvalue solver
We also considered fluid motions that are purely toroidal (horizontal), restricted to a spherical surface of radius r. Keeping the density constant, the linearized momentum equations in a frame rotating at Ω_{carr} are as follows:
where the material derivative D_{t} is given above by Eq. (B.5), and Δ is the horizontal part of the Laplacian. For purely toroidal modes, we introduced the stream function Ψ(θ, ϕ, t) such that
The two above equations can be combined to obtain
We looked for solutions of the form
where m is the longitudinal wavenumber and σ is the (complex) angular frequency. The equation for ψ is of fourthorder and requires four boundary conditions. The condition that the flow vanishes at the poles implies
In order to discretize the problem, we projected ψ onto a basis of associated Legendre polynomials. For the numerical value of the eddy viscosity at the surface, we used the value ν_{t} = 500 km^{2} s^{−1} (Gizon et al. 2020a), unless otherwise specified. The resulting eigenvalue problem was solved for each m using the eigenvalue solver scipy.linalg.eig. We refer to the modes of oscillation obtained in this problem as the “modes of the 1D model”.
Appendix C: Supplementary figures
Fig. C.1.
Rotational frequencies of solar active regions versus latitude, measured in the Carrington frame (from May 2010 to December 2016, Kutsenko 2021). The data (black dots) have been symmetrized in latitude. The cyan ellipses contain 90% of the active regions. The ellipses are extended to higher latitudes by 10° and down to the equator to include flows around active regions. The resulting region in frequency–latitude space is given by the purple contour, which we denote via the equation ω = Ω_{AR}(θ)−Ω_{Carr}. We note that the HMI data used in the main text cover a longer observation period (from May 2010 to September 2020); however, the purple contour is not significantly affected by the very few active regions from the nearly quiet period 2017–2020 (about 5% of all cycle 24 active regions). 
Fig. C.2.
Power spectra for m = 1. Top row: power for the four components , , , and . The purple contours delineate the regions where inflows into active regions produce excess power (see Fig. C.1). The two blue curves show m(Ω − Ω_{Carr})/2π at the surface and at r = 0.95 R_{⊙}. Second row: power spectral density averaged over 0° −30°. The gray curves show the power spectra at full resolution (3.06 nHz), and the black curves show them at a quarter of the resolution. The 95% confidence levels are shown by the red horizontal lines. The cyan curves are for the quietSun period only (2.6 years from 4 February 2018 to 6 September 2020). Third row: power spectral density averaged over 15° −45°. Fourth row: power spectral density averaged over 37.5° −67.5°. In the three lower rows, the dots and the shaded areas (see legend) indicate the significant peaks and the excess power ranges not related to magnetic activity, given in Table A.1. 
Fig. C.10.
Same as Fig. C.2, but for m = 9. 
Fig. C.11.
Same as Fig. C.2, but for m = 10. 
Fig. C.12.
Same as Fig. 1, but using GONG data. The red dots mark the HMI frequencies for comparison. 
Fig. C.13.
Same as Fig. 1, but for the m = 8 equatorial Rossby mode. In the top plot, it is important to notice the range of excess power at low latitudes, below the critical latitudes at the surface. 
Fig. C.14.
Eigenfrequencies in the complex plane for m = 1 from the 2D solver (δ = 0, ν_{t} = 250 km^{2} s^{−1}) and the 1D solver (ν_{t} = 250 km^{2} s^{−1}). Modes with positive imaginary frequencies are selfexcited (unstable). The red vertical line shows the observed frequency of the highlatitude symmetric mode at −86.3 nHz. The red symbols indicate the modes from the models, which have frequencies and surface eigenfunctions close to those observed. 
Fig. C.15.
Same as Fig. 2, but for the complementary velocity components. 
Fig. C.16.
Observed and model radial vorticity for the selected modes of Fig. 2. First, second, and rightmost columns: radial vorticity ζ_{r} = (∇ × u)_{r} for the observations, the 2D model, and the 1D model, respectively. The remaining columns in the middle show meridional cuts of ζ_{r}, radial velocity u_{r}, and the kinetic helicity h_{k} = ⟨u ⋅ ζ⟩ for the 2D model. 
Fig. C.17.
Parameter study (2D model) for different values of the superadiabaticity δ, at fixed ν_{t} = 100 km^{2} s^{−1}. The modes are those shown in Fig. 1. The spiral patterns in u_{ϕ} of the m = 1 highlatitude and m = 2 criticallatitude modes are sensitive to a small change in δ. To obtain a pattern consistent with the observations, δ < 2 × 10^{−7} is implied. The case δ = 10^{−6} is excluded by both the eigenfunctions and the eigenfrequencies. The m = 3 equatorial Rossby mode is almost independent of δ because it is nearly purely horizontal (quasitoroidal). 
Fig. C.18.
Parameter study (2D model) for different values of the turbulent viscosity ν_{t}, for a convection zone that is adiabatically stratified (δ = 0). The modes are those shown in Fig. 1. The frequencies of the m = 1 highlatitude mode and the m = 2 criticallatitude modes are sensitive to the choice of ν_{t}. The smaller values of ν_{t} (≤100 km^{2} s^{−1}) give a better agreement with the observed frequencies (respectively −86.3 nHz and −73.4 nHz). The m = 3 equatorial Rossby modes is essentially insensitive to ν_{t}. 
Movie
Movie 1 associated with Fig. 2 (Fig2_movie_online) (Access here)
All Tables
All Figures
Fig. 1.
Power spectra showing selected modes of oscillation in the Carrington frame. Each column corresponds to a particular m and velocity component, as indicated at the top. Each row shows a different representation of the power spectrum. Top row: the power spectral density is plotted as a function of frequency and latitude. The two blue curves show m(Ω − Ω_{Carr})/2π at the surface and at r = 0.95 R_{⊙}, where Ω(r, θ) is the solar angular velocity in the inertial frame. The purple contour delineates the region in frequency–latitude space affected by inflows into active regions, m(Ω_{AR} − Ω_{Carr})/2π (see Fig. C.1). Second row: the power at each latitude is normalized by its average value over the frequency range between the red bars; this shows that each mode has excess power over a large range of latitudes. The red arrows point to the critical latitudes of ±38° at the surface for the mode with frequency −73 nHz. Third row: the power is averaged over the selected latitude bands specified on the plots, and the frequency resolution is reduced to 12.24 nHz. The red dots point to modes that are not activityrelated (see text) and they are listed in Table A.1. 

In the text 
Fig. 2.
Observed and model eigenfunctions for the modes shown in Fig. 1. Left column: observed velocity ( for the m = 1 and m = 2 modes, for the m = 3 mode). Middle columns: corresponding eigenfunctions of the 2D model for ν_{t} = 100 km^{2} s^{−1} and δ = 0, at the surface and through the central meridian, together with the kinetic energy density. The thick black curves show the critical latitudes. Rightmost column: eigenfunctions of the 1D model at the surface. The retrograde propagation of these modes in the Carrington frame is illustrated as an online movie. The other velocity components are shown in Fig. C.15 and the radial vorticity is shown in Fig. C.16. 

In the text 
Fig. 3.
Mode frequencies in the Carrington frame for the observations and the 2D model (Re σ). The symbols show the observed modes (diamonds for symmetric modes and squares for antisymmetric modes). The red symbols show the highlatitude modes, the orange symbols the criticallatitude modes, and the black symbols the equatorial Rossby modes. The rose and grayshaded areas show the observed frequency ranges of excess power (last column of Table A.1). For reference, the blueshaded area gives the range of rotation rates at the equator between the surface and 0.95 R_{⊙}. The curves give the dispersion relations for the modes of the 2D model with ν_{t} = 100 km^{2} s^{−1} and δ = 0. The red curve is the dispersion relation for the highlatitude modes. The solid and dashedblack curves are for the fundamental (n = 0) and first overtone (n = 1) equatorial Rossby modes. 

In the text 
Fig. C.1.
Rotational frequencies of solar active regions versus latitude, measured in the Carrington frame (from May 2010 to December 2016, Kutsenko 2021). The data (black dots) have been symmetrized in latitude. The cyan ellipses contain 90% of the active regions. The ellipses are extended to higher latitudes by 10° and down to the equator to include flows around active regions. The resulting region in frequency–latitude space is given by the purple contour, which we denote via the equation ω = Ω_{AR}(θ)−Ω_{Carr}. We note that the HMI data used in the main text cover a longer observation period (from May 2010 to September 2020); however, the purple contour is not significantly affected by the very few active regions from the nearly quiet period 2017–2020 (about 5% of all cycle 24 active regions). 

In the text 
Fig. C.2.
Power spectra for m = 1. Top row: power for the four components , , , and . The purple contours delineate the regions where inflows into active regions produce excess power (see Fig. C.1). The two blue curves show m(Ω − Ω_{Carr})/2π at the surface and at r = 0.95 R_{⊙}. Second row: power spectral density averaged over 0° −30°. The gray curves show the power spectra at full resolution (3.06 nHz), and the black curves show them at a quarter of the resolution. The 95% confidence levels are shown by the red horizontal lines. The cyan curves are for the quietSun period only (2.6 years from 4 February 2018 to 6 September 2020). Third row: power spectral density averaged over 15° −45°. Fourth row: power spectral density averaged over 37.5° −67.5°. In the three lower rows, the dots and the shaded areas (see legend) indicate the significant peaks and the excess power ranges not related to magnetic activity, given in Table A.1. 

In the text 
Fig. C.3.
Same as Fig. C.2, but for m = 2. 

In the text 
Fig. C.4.
Same as Fig. C.2, but for m = 3. 

In the text 
Fig. C.5.
Same as Fig. C.2, but for m = 4. 

In the text 
Fig. C.6.
Same as Fig. C.2, but for m = 5. 

In the text 
Fig. C.7.
Same as Fig. C.2, but for m = 6. 

In the text 
Fig. C.8.
Same as Fig. C.2, but for m = 7. 

In the text 
Fig. C.9.
Same as Fig. C.2, but for m = 8. 

In the text 
Fig. C.10.
Same as Fig. C.2, but for m = 9. 

In the text 
Fig. C.11.
Same as Fig. C.2, but for m = 10. 

In the text 
Fig. C.12.
Same as Fig. 1, but using GONG data. The red dots mark the HMI frequencies for comparison. 

In the text 
Fig. C.13.
Same as Fig. 1, but for the m = 8 equatorial Rossby mode. In the top plot, it is important to notice the range of excess power at low latitudes, below the critical latitudes at the surface. 

In the text 
Fig. C.14.
Eigenfrequencies in the complex plane for m = 1 from the 2D solver (δ = 0, ν_{t} = 250 km^{2} s^{−1}) and the 1D solver (ν_{t} = 250 km^{2} s^{−1}). Modes with positive imaginary frequencies are selfexcited (unstable). The red vertical line shows the observed frequency of the highlatitude symmetric mode at −86.3 nHz. The red symbols indicate the modes from the models, which have frequencies and surface eigenfunctions close to those observed. 

In the text 
Fig. C.15.
Same as Fig. 2, but for the complementary velocity components. 

In the text 
Fig. C.16.
Observed and model radial vorticity for the selected modes of Fig. 2. First, second, and rightmost columns: radial vorticity ζ_{r} = (∇ × u)_{r} for the observations, the 2D model, and the 1D model, respectively. The remaining columns in the middle show meridional cuts of ζ_{r}, radial velocity u_{r}, and the kinetic helicity h_{k} = ⟨u ⋅ ζ⟩ for the 2D model. 

In the text 
Fig. C.17.
Parameter study (2D model) for different values of the superadiabaticity δ, at fixed ν_{t} = 100 km^{2} s^{−1}. The modes are those shown in Fig. 1. The spiral patterns in u_{ϕ} of the m = 1 highlatitude and m = 2 criticallatitude modes are sensitive to a small change in δ. To obtain a pattern consistent with the observations, δ < 2 × 10^{−7} is implied. The case δ = 10^{−6} is excluded by both the eigenfunctions and the eigenfrequencies. The m = 3 equatorial Rossby mode is almost independent of δ because it is nearly purely horizontal (quasitoroidal). 

In the text 
Fig. C.18.
Parameter study (2D model) for different values of the turbulent viscosity ν_{t}, for a convection zone that is adiabatically stratified (δ = 0). The modes are those shown in Fig. 1. The frequencies of the m = 1 highlatitude mode and the m = 2 criticallatitude modes are sensitive to the choice of ν_{t}. The smaller values of ν_{t} (≤100 km^{2} s^{−1}) give a better agreement with the observed frequencies (respectively −86.3 nHz and −73.4 nHz). The m = 3 equatorial Rossby modes is essentially insensitive to ν_{t}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.