© L. Gizon et al. 2021

1. Introduction

The free oscillations of a nonrotating spherical star have zero radial vorticity and are called spheroidal modes: they are the pressure (p), surface-gravity (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 (Christensen-Dalsgaard 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, quasi-toroidal 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 quasi-toroidal 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 27-day 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 11-year 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 quasi-toroidal normal modes of the Sun.

2. Observations

We use helioseismic maps of horizontal flows near the solar surface provided by the Stanford ring-diagram pipeline (Bogart et al. 2011a,b) applied to continuous high-resolution 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 low-resolution maps and 80.0° for the high-resolution 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 anti-symmetrizing (“−”) the data with respect to the equator, we computed the Fourier transform of the two velocity components in longitude and in time,

(1)

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 large-scale 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 high-latitude 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 activity-related (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 active-region inflows (Gizon et al. 2001) as modes of oscillation in the power spectra, we defined a region in frequency–latitude space based on the active-region 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 quiet-Sun 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 km2 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 p-mode 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⁻⁶, −2 × 10⁻⁷, 0, 2 × 10⁻⁷, 10⁻⁶ (Fig. C.17) and ν_t = 50, 100, 250, 500 km² s⁻¹ (Fig. C.18). We found that δ = 0 and ν_t = 100 km² s⁻¹ 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 high-latitude inertial modes (Fig. 2a), the critical-latitude 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 high-latitude modes, the orange symbols the critical-latitude modes, and the black symbols the equatorial Rossby modes. The rose- and gray-shaded areas show the observed frequency ranges of excess power (last column of Table A.1). For reference, the blue-shaded 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 km2 s−1 and δ = 0. The red curve is the dispersion relation for the high-latitude modes. The solid and dashed-black curves are for the fundamental (n = 0) and first overtone (n = 1) equatorial Rossby modes.

The high-latitude 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 high-latitude 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.)¹. In the 1D model, only the m = 1 high-latitude modes are self-excited (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.

Critical-latitude 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 global-scale inertial modes in the solar convection zone, within the search range |ω/2π| ≤ 400 nHz and 1 ≤ m ≤ 10. Some of these modes are self-excited 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 critical-latitude 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 p-mode helioseismology. We find that the observed inertial modes are compatible with δ < 2 × 10⁻⁷ and ν_t ≤ 100 km² s⁻¹ 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 (Christensen-Dalsgaard et al. 1996; Muñoz-Jaramillo 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⁻¹ 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⁻¹ 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).

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 Z-CL 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 NAS5-02139 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

Appendix A: Parameters of observed modes

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:

(B.1)

(B.2)

(B.3)

(B.4)

where

(B.5)

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 (Christensen-Dalsgaard 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):

(B.6)

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 stress-free 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 second-order 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 low-frequency 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:

(B.7)

(B.8)

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

(B.9)

The two above equations can be combined to obtain

(B.10)

We looked for solutions of the form

(B.11)

where m is the longitudinal wavenumber and σ is the (complex) angular frequency. The equation for ψ is of fourth-order and requires four boundary conditions. The condition that the flow vanishes at the poles implies

(B.12)

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² s⁻¹ (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 quiet-Sun 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.3. Same as Fig. C.2, but for m = 2.

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

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

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

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

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

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

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 km2 s−1) and the 1D solver (νt = 250 km2 s−1). Modes with positive imaginary frequencies are self-excited (unstable). The red vertical line shows the observed frequency of the high-latitude 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 ur, and the kinetic helicity hk = ⟨u ⋅ ζ⟩ for the 2D model.

Fig. C.17. Parameter study (2D model) for different values of the superadiabaticity δ, at fixed νt = 100 km2 s−1. The modes are those shown in Fig. 1. The spiral patterns in uϕ of the m = 1 high-latitude and m = 2 critical-latitude 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 (quasi-toroidal).

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 high-latitude mode and the m = 2 critical-latitude modes are sensitive to the choice of νt. The smaller values of νt (≤100 km2 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 activity-related (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 km2 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 high-latitude modes, the orange symbols the critical-latitude modes, and the black symbols the equatorial Rossby modes. The rose- and gray-shaded areas show the observed frequency ranges of excess power (last column of Table A.1). For reference, the blue-shaded 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 km2 s−1 and δ = 0. The red curve is the dispersion relation for the high-latitude modes. The solid and dashed-black 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 quiet-Sun 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 km2 s−1) and the 1D solver (νt = 250 km2 s−1). Modes with positive imaginary frequencies are self-excited (unstable). The red vertical line shows the observed frequency of the high-latitude 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 ur, and the kinetic helicity hk = ⟨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 km2 s−1. The modes are those shown in Fig. 1. The spiral patterns in uϕ of the m = 1 high-latitude and m = 2 critical-latitude 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 (quasi-toroidal).

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 high-latitude mode and the m = 2 critical-latitude modes are sensitive to the choice of νt. The smaller values of νt (≤100 km2 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