Solar inertial modes: Observations, identification, and diagnostic promise

The oscillations of a slowly rotating star have long been classified into spheroidal and toroidal modes. The spheroidal modes include the well-known 5-min 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 high-latitude inertial modes, the critical-latitude inertial modes, and the equatorial Rossby modes. In the model, the high-latitude and critical-latitude modes have maximum kinetic energy density at the base of the convection zone, and the high-latitude modes are baroclinically unstable due to the latitudinal entropy gradient. As a first application of inertial-mode helioseismology, we constrain the superadiabaticity and the turbulent viscosity in the deep convection zone.


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, 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 Send offprint requests to: L. Gizon, e-mail: gizon@mps.mpg.de 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. 2020b, 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. 2020a) 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 indepen- In the 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. In the 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. dent 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.

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 hr 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, 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 PS D = |û s j (θ, m, ω)| 2 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. 2020b).
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 lat-itude 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 u + θ (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.

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 −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 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).
The high-latitude inertial modes are analogous to the "wall modes" in plane Poiseuille flows (Gizon et al. 2020b). 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.). 1 In the 1D model, only the m = 1 high-latitude 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.
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. 2020b). 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 frequen- 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.95R . 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 high-latitude modes. The solid and dashed-black curves are for the fundamental (n = 0) and first overtone (n = 1) equatorial Rossby modes.
cies 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).

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 ( 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 (≈ √ 3ν t /τ ≤ 11 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).
Acknowledgements. 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.
A&A proofs: manuscript no. main Appendix A: Parameters of observed modes Notes. In the Carrington frame, the linearized equations for the conservation of momentum, mass, and energy, together with the equation of state, are as follows: 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 D i j = ρν t [∂ i u j + ∂ j u i − 2 3 (∂ k u k )δ i j ] accounts for wave attenuation, where δ i j 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 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 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 low-frequency solutions. We refer to the modes of oscillation obtained in this problem as the "modes of the 2D model".

Appendix 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 (Ω sin 2 θ) ∂Ψ ∂φ = ν t ∆ 2 Ψ. (B.10) We looked for solutions of the form 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 ψ = dψ dθ = 0 at θ = 0 and π. (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 2 s −1 (Gizon et al. 2020b), 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".  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).   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 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 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 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). 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 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 .