Open Access
Issue
A&A
Volume 682, February 2024
Article Number A39
Number of page(s) 11
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/202348380
Published online 31 January 2024

© The Authors 2024

Licence Creative CommonsOpen 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.

This article is published in open access under the Subscribe to Open model.

Open access funding provided by Max Planck Society.

1. Introduction

Inertial modes are global-scale low-frequency modes of oscillation in a rotating fluid whose restoring force is the Coriolis force (e.g., Greenspan et al. 1968). Recently, various kinds of inertial modes have been observed on the Sun. These include the equatorial Rossby modes (Löptien et al. 2018; Liang et al. 2019; Proxauf et al. 2020; Mandal & Hanasoge 2020; Hathaway & Upton 2021), the critical-latitude modes, and the high-latitude modes (Gizon et al. 2021). All of these modes propagate in a retrograde direction (opposite to solar rotation) when seen from the Carrington rotation frame. It is expected that these inertial modes can be used to probe the interior of the Sun (e.g., Goddard et al. 2020; Gizon et al. 2021; Bekki et al. 2022a). In particular, using the baroclinically unstable m = 1 high-latitude mode, Gizon et al. (2021) derived the observational constraint of the mean superadiabaticity δ, one of the most important unknown parameters in the Sun’s convection zone. They deduced δ < 2 × 10−7, implying that the Sun’s convection zone is closer to adiabatic than typically assumed.

Recently, Hanson et al. (2022, hereafter hanson2022) have reported to detect another family of modes near the surface of the Sun, namely, retrograde-propagating modes of equatorially antisymmetric radial vorticity, ζr. They can be most clearly seen in the l = m + 1 component of the radial vorticity power spectrum, where l is the spherical harmonic degree and m is the azimuthal order. hanson2022 reported that these modes propagate in a retrograde direction about three times faster than the equatorial Rossby modes with the same azimuthal order, m, which cannot be explained by the classical Rossby modes. In this paper, we refer to them as HHS22 modes.

A possible identification of the HHS22 modes was first provided by Triana et al. (2022, hereafter triana2022) using a linear model of rotating fluid in a spherical shell. The HHS22 modes are identified as a particular class of non-toroidal inertial modes which involve substantial radial motions in the middle convection zone. The computed linear dispersion relation of the HHS22 modes shows a good agreement with their observed frequencies. However, the linear model of triana2022 was highly simplified; for instance, the model assumes an incompressible (uniform density) fluid and uniformly rotating convection zone, which are not appropriate in the case of the Sun.

A similar mode identification was later reported by Bhattacharya & Hanasoge (2023, hereafter bhattacharya2022) in which the solar-like background stratification was taken into account using the anelastic approximation. However, the latitudinal differential rotation of the Sun was still omitted in their model for simplicity, which is known to have a substantial impact on the properties of inertial modes (e.g., Baruteau & Rieutord 2013; Gizon et al. 2021; Bekki et al. 2022a; Fournier et al. 2022; Philidet & Gizon 2023; Bhattacharya et al. 2023). bhattacharya2022 found that the computed frequencies of the HHS22 modes are lower than the observed ones by about 100 nHz, in contrast to the match found by triana2022. It is necessary to understand the origin of this discrepancy and to investigate the missing physics to properly reproduce the observed features of the HHS22 modes.

In this paper, we study the properties of the HHS22 modes using a more realistic linear model of the Sun’s convection zone (Bekki et al. 2022a), which takes into account both the solar density stratification and the solar differential rotation determined by global helioseismology (Larson & Schou 2018). We then show that their frequencies are strongly affected by the superadiabaticity δ in the bulk convection zone.

Recently, the solar inertial modes are also studied using fully nonlinear simulations of the rotating convection. Bekki et al. (2022b) have mostly focused on the equatorial Rossby modes and the columnar convective (also known as thermal Rossby) modes. Matilsky et al. (2022) have shown that the equatorial Rossby modes can contribute to the confinement of the solar tachocline via dynamo action in the radiative interior. However, the HHS22 modes have never been studied yet. In the appendix of this paper, we also show that the HHS22 modes can be found to exist in the fully-nonlinear simulations.

The organization of the paper is as follows. In Sect. 2, our linear eigenmode solver is explained. The effects of the background density stratification and the solar differential rotation on the HHS22 modes are examined in Sect. 3.1. The effects of turbulent viscosity is briefly discussed in Sect. 3.2. We then show how the solar observations can be used to infer the superadiabaticity in the bulk convection zone δcz in Sect. 3.3. In Appendix A, we further report that the HHS22 modes can be detected in the fully-nonlinear simulations of rotating convection. The conclusions are summarized in Sect. 4.

2. Method: Linear eigenmode analysis

We used the code developed by Bekki et al. (2022a), which numerically solves the linear eigenvalue problem for a rotating fluid in a spherical shell 0.71 R < r < 0.985 R. Here, R is the solar radius. The model is hydrodynamic; that is, the effects of magnetic field are neglected for simplicity. The eigenvalue equations are solved for azimuthal orders 1 ≤ m ≤ 19. To study the HSS22 modes, we seek for the retrograde-propagating modes (ℜ[ω]< 0) whose eigenfunctions satisfy the following criteria: (1) the eigenfunction of radial velocity, vr, is dominantly l = m in the middle convection zone, and has no radial node at the equator; (2) the eigenfunction of latitudinal velocity, vθ, is dominantly l = m + 1 at the surface, and has no radial node at low latitudes; and 3) the eigenfunction of longitudinal velocity, vϕ, is dominantly l = m or l = m + 2 at the surface. When the above criteria are satisfied by multiple eigenmodes, we select the least-damped mode (with largest growth rate ℑ[ω]) among them. For further details, see Bekki et al. (2022a, Sect. 2).

In this paper, we first carry out three sets of linear eigenmode calculations with different model setups. The first setup (model 1) consists of incompressible fluid and uniform rotation, which were assumed in the previous study of triana2022. In the second setup (model 2), we still assume the uniform rotation but include the solar background density stratification, which is similar to the numerical setup of bhattacharya2022. The third setup (model 3) finally takes into account both the realistic solar density stratification and the solar differential rotation determined by global helioseismology (Larson & Schou 2018). These are summarized in Table 1.

Table 1.

Summary of the linear analysis model setups.

3. Results

3.1. Comparison with previous studies

In this section, we compare the results from models 1–3. In all cases, we include the spatially constant viscous diffusivity of ν = 1012 cm2 s−1. For simplicity, we assume that the background is purely adiabatic (i.e., δ = 0). Hence, there is no entropy variation neither in the radial nor the latitudinal directions.

3.1.1. Effects of density stratification

Figure 1 shows the linear dispersion relations of the HHS22 modes computed from the models 1–3, measured in the Carrington frame rotating at Ω0/2π = 456.0 nHz. We find that the dispersion relation from the model 1 matches almost exactly to that of triana2022; the differences are found to be less than 1.5%. It is confirmed that the observed frequencies of the HHS22 modes can be nicely fitted by the dispersion relation of our model 1 where the over-simplifying assumptions are used. Interestingly, however, when the solar density stratification is taken into account, the dispersion relation of the HHS22 modes deviates from the observations: in our model 2, the frequencies become much lower than the observed ones (by about 100 nHz at m = 10). These frequencies are consistent with the results reported in bhattacharya2022 with the differences on the order of few percent. The small discrepancy can be attributed to the differences in the model setups such as the different radial profiles of the turbulent diffusivities and the inclusion of the radiative interior below the convection zone.

thumbnail Fig. 1.

Dispersion relation of the HHS22 modes obtained from the linear analysis for 1 ≤ m ≤ 19. Green points represent the results from our model 1 where we assume an incompressible (constant density) fluid and an uniform rotation. Blue points represent the results from our model 2 where the solar density stratification is included but the uniform rotation is still assumed. Orange points represent the results from our model 3 where both the density stratification and the differential rotation of the Sun are included. Red points also show the results from model 3 but with a weakly subadiabatic bulk convection zone (δcz = −5 × 10−7) and a strongly superadiabatic near-surface layer (δsf = 10−3). Lime and cyan solid curves show the results from triana2022 and from bhattacharya2022, respectively. For comparison, we also show the observed frequencies of the HHS22 modes reported in hanson2022 where the gray diamonds and squares denote the measurements with ring-diagram analysis (RDA) and mode-coupling analysis (MCA), respectively. All the frequencies are measured in the Carrington frame rotating at Ω0/2π = 456 nHz.

Figure 2 shows the meridional eigenfunctions of the HHS22 mode at m = 10. All the eigenfunctions are normalized such that the maximum of vθ is 1.0 m s−1 at the surface, as suggested by the observation hanson2022. Those of model 1 (Fig. 2a) agree with the results of triana2022 (see Fig. 4 top panels therein). As already reported, the HHS22 modes have a radial vorticity ζr which is north-south antisymmetric across the equator. At the equator, it is seen that the latitudinal diverging (converging) motion at the surface involves a substantial radial upflow (downflow) in the middle convection zone. This clearly manifests that the HHS22 modes are not toroidal at all. Therefore, their mode properties are different from the l = m + 1 classical Rossby modes which are quasi-toroidal (even though their surface eigenfunctions appear similar to those of l = m + 1 classical Rossby modes).

thumbnail Fig. 2.

Meridional cuts of the eigenfunctions of the m = 10 HHS22 mode obtained from the linear analysis. The velocity, pressure, and vorticity eigenfunctions are expressed as ν(r, θ)exp[i( − ωt)], p1(r, θ)exp[i( − ωt)], and ζ(r, θ)exp[i( − ωt)], and the solutions are shown in the meridional plane at t = 0 and ϕ = 0. The units of the colorbars are m s−1 for three velocity components, 104 dyn cm−2 for the pressure perturbation, and 10−8 s−1 for the vorticity components. The eigenfunctions are normalized such that the maximum of vθ is 1.0 m s−1. Panels (a–c) show the results from models 1–3, respectively. Note that the background entropy stratification is adiabatic in all cases. In panel (c), the black solid curves show the locations of the critical latitudes.

To further investigate the consequences of non-toroidalness, we show the eigenfunctions of z-vorticity ζz in the rightmost panels of Fig. 2, where z denotes a direction parallel to the rotational axis. It is seen that the HHS22 mode involves a north-south symmetric z-vortical motion near the equator, in addition to the north-south antisymmetric r-vortical motion. We note that the amplitude of ζz is comparable to that of ζr. This z-vortical motion invokes additional β-effects, namely, the topographic β-effect originating from the spherical curvature (e.g., Busse 2002) and the compressional β-effect originating from the background density stratification (e.g., Glatzmaier et al. 2009; Verhoeven & Stellmach 2014; Gastine et al. 2014; Bekki et al. 2022a). Outside the tangential cylinder of the Sun’s convection zone, it is known that the compressional β-effect dominates over the topographic β-effect (Bekki 2022, see Fig. 3.32 therein). At the equator, the z-vorticity equation can be expressed as:

ζ z t = β comp v r + [ ] , with β comp = 2 Ω 0 H ρ , $$ \begin{aligned} \frac{\partial \zeta _{z}}{\partial t}=\beta _{\mathrm{comp}}v_{r} + [ \ldots ], \quad \mathrm{with} \quad \beta _{\mathrm{comp}}=-\frac{2\Omega _{0}}{H_{\rho }}, \end{aligned} $$(1)

where Hρ denotes the density scale height of the background and we only retain the term related to the compressional β-effect. When the simplifying assumption of incompressible fluid is used in model 1, the compressional β-effect is ignored (βcomp = 0). However, when the density stratification is included in model 2, a negative value of βcomp promotes a prograde propagation which acts against the retrograde propagation of the modes (Glatzmaier & Gilman 1981). This decreases the retrograde propagation frequencies of the HHS22 modes compared to the model 1, and consequently, the dispersion relation deviates from the observations. Our study suggests that a great agreement reported by triana2022 is largely due to the unrealistic assumption of incompressible fluid in the Sun’s convection zone, which ignores the compressional β-effect.

3.1.2. Effects of solar differential rotation

It is shown in Fig. 1 that the inclusion of solar differential rotation modifies the dispersion relation of the HHS22 modes (from model 2 to model 3). The significant modification occurs at higher m, where the retrograde propagation frequencies become larger than the observed values by about 50 − 80 nHz. This is a direct consequence of the Doppler frequency shift by the differential rotation.

The eigenfunctions of the m = 10 HHS22 mode from the model 3 are shown in Fig. 2c. It is known that the inclusion of latitudinal differential rotation gives rise to critical latitudes where the phase speed of the mode becomes equal to the local differential rotation speed (Baruteau & Rieutord 2013; Gizon et al. 2020; Bekki et al. 2022a; Fournier et al. 2022; Philidet & Gizon 2023). The locations of the critical latitudes are denoted by black solid curves. Compared to those from the model 2 (Fig. 2b), the mode eigenfunctions from the model 3 are distorted by the existence of critical layers.

To better see the impact of the critical layers, we show the radial profiles of the root-mean-square (RMS) velocity eigenfunctions of the HHS22 modes for different azimuthal orders m in Fig. 3a. For small azimuthal orders m (≲6), the radial eigenfunctions from models 2 and 3 are similar: In both cases, the RMS velocity increases with radius and peaks at the surface. However, as m increases, the velocity eigenfunctions from the model 3 are more and more confined into the lower convection zone, leading to a significant deviation from those of model 2. This is because the HHS22 modes tend to be more and more localized around the critical layers. For sufficiently large m (≳10), the critical layers exist near the base of the convection zone at the equator. Figure 3b further compares the volume-integrated kinetic energy of the HHS22 modes between the model 2 and model 3 where the maximum horizontal velocity amplitudes are fixed to 1.0 m s−1 at the surface in both cases. It is shown that, as a consequence of the mode confinement near the base, the HHS22 modes have much more kinetic energy in model 3 than in model 2. Therefore, our study suggests that the solar differential rotation needs to be properly taken into account in order to evaluate the dynamical importance of the HHS22 modes in the Sun’s convection zone. In the remaining part of the paper, we only consider the most realistic model 3, which takes into account both the solar density stratification and the solar differential rotation.

thumbnail Fig. 3.

(a) Root-mean-square (RMS) amplitudes of the velocity eigenfunctions of the HHS22 modes from the linear analysis, where the average is taken over the spherical surfaces. Dashed and solid curves show the results from model 2 (without differential rotation) and model 3 (with differential rotation). Different colors show different azimuthal orders m. The eigenfunctions are normalized such that the maximum horizontal velocity at r = 0.985 R is 1.0 m s−1. (b) Spectra of volume-integrated kinetic energy of the HHS22 modes. The blue and red points denote those from models 2 and 3, respectively.

3.2. Dependence on turbulent viscosity

In this section, the dependence of the HHS22 modes on the turbulent viscosity ν is briefly discussed. Here, we use the model 3 with the adiabatic background (i.e., δ = 0 throughout the convection zone). Figure 4 compares the eigenfequencies of the HHS22 modes obtained for two different values of turbulent viscosity, ν = 1011 cm2 s−1 (blue) and ν = 1012 cm2 s−1 (red). For simplicity, ν is assumed to be spatially constant throughout the convection zone. Although the mode linewidths are decreased by about 60 − 100 nHz with the decrease of ν, the dispersion relation of the HHS22 modes is only marginally affected (difference of less than 3%). Regardless of the values of turbulent viscosity used, the discrepancies between our linear model and the observations remain in the propagation frequencies of the HHS22 modes at 8 ≤ m ≤ 14.

thumbnail Fig. 4.

(a) Propagation frequencies and (b) linewidths of the HHS22 modes obtained from the linear calculations with two different values of turbulent viscosity, ν. Blue and red points show the results with spatially uniform viscosity of ν = 1011 cm2 s−1 and for 1012 cm2 s−1, respectively. The model 3 is used (with solar density stratification and solar differential rotation) and the background is assumed to be adiabatic (δ = 0). The observed values reported by hanson2022 are shown by gray diamonds.

3.3. Dependence on superadiabaticity δ

Next, we investigate the effects of the non-adiabatic stratification in the Sun’s convection zone on the HHS22 modes. A deviation from the adiabatic stratification is measured by the superadiabaticity δ = ∇ − ∇ad, where ∇ = dlnT/dlnp. The superadiabaticity in the Sun’s convection zone is typically estimated to be very small, δ ≈ 𝒪(10−6) (Ossendrijver 2003), except for a very thin near-surface layer where the solar granulation is vigorously driven by the strong surface radiative cooling (e.g., Nordlund et al. 2009). According to the solar standard model S (Christensen-Dalsgaard et al. 1996; Stix 2002), δ is expected to increase up to 𝒪(10−3) at r ≈ 0.99 R and even further up to 𝒪(10−1) in the photosphere.

In this section, we carried out a set of linear calculations with varying superadiabaticity δ. The most realistic model (model 3) was used with the turbulent viscosity of ν = 1012 cm2 s−1. We used the following radial profiles of δ(r), which mimics a sharp transition from the close-to-adiabatic bulk convection zone to the strongly superadiabatic surface,

δ ( r ) = δ cz + ( δ sf δ cz ) exp [ ( r r max d sf ) 2 ] , $$ \begin{aligned} \delta (r)=\delta _{\mathrm{cz}} + (\delta _{\mathrm{sf}}-\delta _{\mathrm{cz}}) \exp {\left[-\left(\frac{r-r_{\mathrm{max}}}{d_{\mathrm{sf}}} \right)^{2} \right]}, \end{aligned} $$(2)

where δcz and δsf denote the values of superadiabaticity in the bulk of the convection zone and near the top surface, respectively. We use the values: dsf = 0.015 R and rmax = 0.985 R. Therefore, the strong superadiabaticities are localized in a thin layer near the top boundary of our numerical domain.

3.3.1. Impact of near-surface superadiabaticity

Figure 5a shows the linear dispersion relations of the HHS22 modes computed for different values of the near-surface superadiabaticity, δsf. In all cases, the bulk of the convection zone (below 0.95 R) is fixed to be purely adiabatic, δcz = 0. Within the parameter range investigated here, the frequencies of the HHS22 modes are shown to be almost independent of δsf.

thumbnail Fig. 5.

(a) Dispersion relation of the HHS22 modes computed for different values of superadiabaticity near the surface δsf. The bulk of the convection zone is assumed to be adiabatic, δcz = 0. (b–e) Eigenfunctions of the m = 10 HHS22 mode computed for δsf = 10−6 (left column) and for δsf = 10−3 (right column). Panels b–d show cuts of radial velocity vr (in m s−1), entropy perturbation s1 (in erg g−1 K−1), and z-vorticity ζz (in 10−8 s−1) in the equatorial plane (along the rotational axis) seen from the north pole. The black dashed curves denote the height r = 0.95 R, above which the strongly superadiabatic layer is located. (e) Mollweide projection of the radial vorticity eigenfunction ζr at the top boundary r = 0.985 R (in 10−8 s−1). All the eigenfunctions are normalized in the same way as in Fig. 2.

Figure 5b compares the eigenfunctions of the m = 10 HHS22 mode from the two representative cases with δsf = 10−6 and with δsf = 10−3. Panels (b–d) show the eigenfunctions of radial velocity vr, entropy perturbation s1, and z-vorticity ζr, in the equatorial plane, respectively. The radial velocity vr is almost unchanged given it is concentrated near the base of the convection zone. However, the strong entropy perturbation s1 is generated in the near-surface superadiabatic layer as δsf increases (Fig. 5c). The buoyancy force associated with this entropy perturbation drives the z-vortical motion in this near-surface layer (Fig. 5d). Nonetheless, it is demonstrated that the overall structure of the mode eigenfunctions in the bulk of the convection zone (below 0.95 R) is only little affected by the inclusion of the strongly superadiabatic near-surface layer. In fact, the surface eigenfunction of radial vorticity ζr remains unchanged by the increase of δsf up to 10−3 (Fig. 5e).

We must note that, in the real Sun, δ is expected to further increase above 0.985 R up to 𝒪(10−1) in the photosphere (e.g., Christensen-Dalsgaard et al. 1996; Stix 2002). An inclusion of this realistic photosphere may have a non-negligible impact on the HHS22 modes. However, resolving this very thin surface layer is numerically challenging and thus beyond the scope of this paper.

3.3.2. Impact of superadiabaticity in the bulk convection zone

Next, we vary the superadiabaticity in the bulk of the convection zone δcz within the range of ±10−6, while using a fixed value for the near-surface superadiabaticity δsf = 10−3. The radial profiles of δ(r) used in this study are shown in Fig. 6a.

thumbnail Fig. 6.

(a) Radial profiles of the superadiabaticity δ(r) defined by Eq. (2) for different values of the bulk superadiabaticity δcz. The superadiabaticity in the near-surface layer is fixed to δsf = 10−3. (b) Dispersion relations of the HHS22 modes computed for different superadiabaticity profiles. The observed frequencies reported by hanson2022 are denoted by gray diamonds. (c) Meridional eigenfunctions from the case with weakly subadiabatic bulk convection zone (δcz = −5 × 10−7). Black solid and gray dashed curves show the locations of the critical latitudes (by differential rotation) and the turning surfaces (by subadiabatic stratification), respectively.

Figure 6b manifests a striking sensitivity of the dispersion relation of the HHS22 modes to a tiny change in δcz. It is found that, as the stratification becomes more superadiabatic (subadiabatic), their dispersion relation shifts towards more positive (negative) direction in frequency. In other words, the HHS22 modes propagate in a retrograde direction with slower (faster) phase speed when δcz > 0 (< 0). This frequency shift can be understood by considering whether the buoyancy force associated with the radial motion acts as an additional restoring force or the opposite. The similar behavior is already reported and discussed in the case of columnar convective (thermal Rossby) modes by Bekki et al. (2022a as described in Sect. 5), in which their prograde propagation frequencies become slower (faster) when δ > 0 (< 0). It is shown in Fig. 6b that the observed frequencies can be nicely fitted by the linear dispersion relation with weakly subadiabatic bulk convection zone −5 × 10−7 ≲ δcz ≲ −2.5 × 10−7. This is within the range of the observational constraint of δ independently derived by Gizon et al. (2021) based on the m = 1 high-latitude inertial mode.

For the range of subadiabaticity −5 × 10−7 ≲ δcz ≲ −2.5 × 10−7, the Brunt-Väisälä frequency N = g | δ cz | / H p $ N=\sqrt{g|\delta_{\mathrm{cz}}|/H_{p}} $ is estimated to be N/2π ≈ 240 − 340 nHz (with g ≈ 520 m s−2 and Hp ≈ 57 Mm near the base of the convection zone), which is comparable to the mode frequencies at 8 ≤ m ≤ 14. Therefore, the HHS22 modes are expected to behave as gravito-inertial modes in which both Coriolis and buoyancy forces act as restoring forces. Some gravito-inertial modes are known to be trapped by turning surfaces (Friedlander & Siegmann 1982; Dintrans et al. 1999; Mirouh et al. 2016). The eigenfunctions of the m = 10 HHS22 mode are shown in Fig. 6c in the case of δcz = −5 × 10−7. We find that the turning surfaces (denoted by gray dashed curves) are located at higher latitudes than the critical latitudes of the differential rotation. Therefore, the turning surfaces only play a limited role in trapping the HHS22 modes which are already strongly confined into the equatorial region by the critical latitudes. Thus, in the bulk of the convection zone, no significant difference can be seen in the mode eigenfunctions from the adiabatic case (Fig. 2c). Near the surface, by contrast, strong z-vortical motions are apparent as a consequence of the strong superadiabaticity.

4. Summary and discussion

In this study, we investigated the properties of the l = m + 1 radial vorticity modes recently observed by hanson2022 (referred to as HHS22 modes in this paper). We used the linear eigenmode solver of the Sun’s convection zone developed by Bekki et al. (2022a). Our model can successfully reproduce the previous results of triana2022 and bhattacharya2022 when the simplifying assumptions are used. We found that, in contrast to the classical l = m + 1 Rossby modes, the HHS22 modes are very sensitive to the background density stratification. This is because the HHS22 modes are essentially non-toroidal. They involve substantial z-vortical motion near the equator and the compressional β-effect plays a significant role. The Sun’s differential rotation is also found to affect the HHS22 modes by introducing the viscous critical layers, leading to a confinement of the HHS22 modes near the base of the convection zone at large azimuthal orders m ≳ 10 (Fig. 3).

We further examined possible effects of the background superadiabaticity, δ, on the HHS22 modes. In this study, we used a highly-simplified radial function of δ(r) which changes from a close-to-adiabatic value in the bulk of the convection zone (r ≲ 0.95 R) to a strongly superadiabatic value near the top surface (0.95 R ≲ r ≤ 0.985 R). Surprisingly, the strong superadiabaticity near the top surface δsf does not affect the properties of the HHS22 modes such as their propagation frequencies and the radial vorticity eigenfunction at the surface (Fig. 5). On the other hand, we found that their dispersion relation is quite sensitive to a small variation in the bulk superadiabaticity δcz. This difference can be attributed to the fact that the radial motions of the HHS22 modes dominantly exist in the lower convection zone as a consequence of the mode confinement towards the base by the differential rotation (Figs. 2c and 6c). We showed that in order to explain the observed frequencies of the HHS22 modes, the bulk of the convection zone needs to be weakly subadiabatic, namely, −5 × 10−7 ≲ δcz ≲ −2.5 × 10−7. This constraint is consistent with but tighter than the previous constraint derived by Gizon et al. (2021) using the m = 1 high-latitude inertial mode (δ < 2 × 10−7).

Our result suggests that the Sun’s bulk convection zone is likely much less superadiabatic than typically thought and possibly be even subadiabatic. This unconventional conclusion needs to be tested using an improved linear model with more realistic effects included. For instance, we neglected the effects of magnetic field in this study, which are known to affect the properties of the Rossby modes (e.g., Zaqarashvili et al. 2021). Moreover, we still cannot rule out the possibility that the above conclusion could be affected by an inclusion of the realistic solar photosphere (above 0.985 R) with a very strong superadiabaticity on the order of δ ≈ 𝒪(10−1). We also note that the spatially uniform δ below the strongly superadiabatic near-surface layer might be over-simplifying. A further parameter survey on various radial profiles of δ(r) will be required (Dey et al., in prep.).

Recent direct numerical simulations of solar convection have reported a formation of the weakly subadiabatic layer near the base of the convection zone as a consequence of the non-local convective heat transport (Käpylä & Rheinhardt 2017; Hotta 2017; Bekki et al. 2017; Karak et al. 2018; Nelson et al. 2018; Hotta et al. 2022; Käpylä 2024; see also Appendix A). Our study implies that the HHS22 modes are likely gravito-inertial modes originating from this weakly subadiabatic lower convection zone. A further investigation on the HHS22 modes will be desired both from observational and theoretical perspectives to better understand this weakly subadiabatic layer, as it might help us to explain the Sun’s anomalously weak convective power at large spatial scales (e.g., Hanasoge et al. 2012) and to resolve the Sun’s convective conundrum (e.g., O’Mara et al. 2016; Hotta et al. 2023).

Acknowledgments

The author thanks an anonymous referee for constructive comments. Y.B. is also grateful to P. Dey, R. Cameron and L. Gizon for helpful discussions and useful comments on the initial manuscript. Y.B. acknowledges a support from ERC Synergy Grant WHOLE SUN 810218. All the numerical computations were performed at GWDG and the Max-Planck supercomputer RZG in Garching.

References

  1. Baruteau, C., & Rieutord, M. 2013, J. Fluid. Mech., 719, 47 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bekki, Y. 2022, Ph.D. Thesis, Georg-August-Universität Göttingen, Germany [Google Scholar]
  3. Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bekki, Y., Cameron, R. H., & Gizon, L. 2022a, A&A, 662, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bekki, Y., Cameron, R. H., & Gizon, L. 2022b, A&A, 666, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bhattacharya, J., & Hanasoge, S. M. 2023, ApJS, 264, 21 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bhattacharya, J., Hanson, C. S., Hanasoge, S. M., & Sreenivasan, K. R. 2023, ArXiv e-prints [arXiv:2308.12766] [Google Scholar]
  8. Busse, F. H. 2002, Phys. Fluids, 14, 1301 [NASA ADS] [CrossRef] [Google Scholar]
  9. Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [Google Scholar]
  10. Dintrans, B., Rieutord, M., & Valdettaro, L. 1999, J. Fluid. Mech., 398, 271 [NASA ADS] [CrossRef] [Google Scholar]
  11. Fournier, D., Gizon, L., & Hyest, L. 2022, A&A, 664, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Friedlander, S., & Siegmann, W. L. 1982, Geophys. Astrophys. Fluid Dyn., 19, 267 [NASA ADS] [CrossRef] [Google Scholar]
  13. Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156 [Google Scholar]
  14. Gastine, T., Heimpel, M., & Wicht, J. 2014, Phys. Earth Planet. Inter., 232, 36 [CrossRef] [Google Scholar]
  15. Gizon, L., Fournier, D., & Albekioni, M. 2020, A&A, 642, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Glatzmaier, G. A., & Gilman, P. A. 1981, ApJS, 45, 381 [NASA ADS] [CrossRef] [Google Scholar]
  18. Glatzmaier, G., Evonuk, M., & Rogers, T. 2009, Geophys. Astrophys. Fluid Dyn., 103, 31 [NASA ADS] [CrossRef] [Google Scholar]
  19. Goddard, C. R., Birch, A. C., Fournier, D., & Gizon, L. 2020, A&A, 640, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Greenspan, H., Batchelor, C., Ablowitz, M., et al. 1968, The Theory of Rotating Fluids (Cambridge: Cambridge University Press) [Google Scholar]
  21. Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928 [Google Scholar]
  22. Hanson, C. S., Hanasoge, S., & Sreenivasan, K. R. 2022, Nat. Astron., 6, 708 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hathaway, D. H., & Upton, L. A. 2021, ApJ, 908, 160 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
  25. Hotta, H., Kusano, K., & Shimada, R. 2022, ApJ, 933, 199 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hotta, H., Bekki, Y., Gizon, L., Noraz, Q., & Rast, M. 2023, Space Sci. Rev., 219, 77 [NASA ADS] [CrossRef] [Google Scholar]
  27. Käpylä, P. J. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202348325 [Google Scholar]
  28. Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23 [CrossRef] [Google Scholar]
  29. Karak, B. B., Miesch, M., & Bekki, Y. 2018, Phys. Fluids, 30, 046602 [NASA ADS] [CrossRef] [Google Scholar]
  30. Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29 [Google Scholar]
  31. Liang, Z.-C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
  33. Mandal, K., & Hanasoge, S. 2020, ApJ, 891, 125 [CrossRef] [Google Scholar]
  34. Matilsky, L. I., Hindman, B. W., Featherstone, N. A., Blume, C. C., & Toomre, J. 2022, ApJ, 940, L50 [NASA ADS] [CrossRef] [Google Scholar]
  35. Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mirouh, G. M., Baruteau, C., Rieutord, M., & Ballot, J. 2016, J. Fluid. Mech., 800, 213 [NASA ADS] [CrossRef] [Google Scholar]
  37. Nelson, N. J., Featherstone, N. A., Miesch, M. S., & Toomre, J. 2018, ApJ, 859, 117 [Google Scholar]
  38. Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
  39. O’Mara, B., Miesch, M. S., Featherstone, N. A., & Augustson, K. C. 2016, Adv. Space Res., 58, 1475 [CrossRef] [Google Scholar]
  40. Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
  41. Philidet, J., & Gizon, L. 2023, A&A, 673, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Proxauf, B., Gizon, L., Löptien, B., et al. 2020, A&A, 634, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Stix, M. 2002, The Sun: An Introduction (Heidelberg: Springer Berlin) [Google Scholar]
  44. Triana, S. A., Guerrero, G., Barik, A., & Rekier, J. 2022, ApJ, 934, L4 [NASA ADS] [CrossRef] [Google Scholar]
  45. Verhoeven, J., & Stellmach, S. 2014, Icarus, 237, 143 [NASA ADS] [CrossRef] [Google Scholar]
  46. Zaqarashvili, T. V., Albekioni, M., Ballester, J. L., et al. 2021, Space Sci. Rev., 217, 15 [Google Scholar]

Appendix A: Nonlinear simulation of solar rotating convection

In this appendix, we report an additional analysis of the fully-nonlinear numerical simulations of rotating convection by Bekki et al. (2022b). We must note that these direct numerical simulations tend to show several inconsistencies with the observations regarding the profiles of the large-scale mean flows (e.g., Nelson et al. 2018) and the surface velocity power spectra (e.g., Miesch et al. 2008; Hanasoge et al. 2012). This as-yet unsolved problem is known as the "convective conundrum" (e.g., O’Mara et al. 2016; Hotta et al. 2023). The purpose of this appendix is, therefore, not to reproduce the observations in a self-consistent manner but only to report that the HHS22 modes can be found in the nonlinear simulations of Bekki et al. (2022b).

A.1. Method

We use the data from the nonlinear numerical simulations of solar-like rotating convection carried out by Bekki et al. (2022b). In their simulations, a full set of compressible hydrodynamic equations are solved in a spherical shell (0.71R < r < 0.96R) rotating at the solar rotation rate. The luminosity is artificially reduced from the solar value by a factor of 20 to achieve the strongly rotationally-constrained regime (low Rossby number regime) and to obtain the solar-like differential rotation with faster equator and slower poles (e.g., Gastine et al. 2013). In Bekki et al. (2022b), total six simulations were carried out with the same numerical setup but with different initial random perturbations. Each simulation has been evolved for more than 25 solar years and we analyzed the 15-year-long data after the large-scale mean flows become statistically stationary. The results were averaged over these six realizations to increase the signal-to-noise ratio.

A.2. Results

Figures A.1(a–c) show the near-surface power spectra of the l = m component of longitudinal velocity vϕ, l = m + 1 component of latitudinal velocity vθ, and l = m + 1 component of radial vorticity ζr from the nonlinear simulations. The spectra are computed within a Carrington frame. Cyan diamonds denote the observed frequencies of the HHS22 modes hanson2022 and red points show the linear dispersion relation from the model 3 with the weakly subadiabatic bulk convection zone δcz = −5 × 10−7 which best-fits the observation (we refer to this model 3-sub hereafter). In the l = m spectra of vϕ and in the l = m + 1 spectra of ζr, the prograde-propagating columnar convective (thermal Rossby) modes are dominantly seen (Bekki et al. 2022b) but the HHS22 modes are not clearly visible. On the other hand, in the l = m + 1 spectra of vθ, we can see a concentration of the surface velocity power near the expected frequency range of the HHS22 modes.

thumbnail Fig. A.1.

Power spectra from the nonlinear rotating convection simulation of (a) the l = m component of longitudinal velocity vϕ, (b) the l = m + 1 component of latitudinal velocity vθ, and (c) the l = m + 1 component of radial vorticity ζr near the surface r = 0.95R. The spectra are computed in the Carrington frame. The power is normalized at each m. Cyan diamonds denote the observed frequencies of the HHS22 modes hanson2022. Orange points show the dispersion relation of the HHS22 modes obtained from the linear calculation (model 3-sub with δcz = −5 × 10−7 and δsf = 10−3). The power ridge associated with the columnar convective modes is denoted by black arrows.

Figure A.2 shows the same power spectra as Fig. A.1 but at the fixed azimuthal order m = 10. It is shown that the columnar convective mode has a strong power peak at ω/2π ≈ 240 nHz with narrow linewidth of about 50 nHz, which can be seen in all the spectra (Figs. A.2a–c). In the l = m + 1 spectrum of vθ, another strong power concentration is seen at around ω/2π ≈ −200 nHz aside from that of the columnar convective modes, which we identified as a HHS22 mode. To further characterize the mode properties, we perform a Lorentzian fit to this power spectrum as shown by a cyan curve in Fig. A.2(b). The measured peak frequencies and linewidths of the HHS22 modes in the nonlinear simulations are reported in Table A.1 for 1 ≤ m ≤ 15. For comparison, we also report the values computed from the linear model 3-sub. It is found that, in our nonlinear simulations, the HHS22 modes have very broad linewidths of about 100 − 200 nHz, which are about twice larger than those from the linear analysis. This suggests that they are very short-lived in the nonlinear simulations.

thumbnail Fig. A.2.

Same power spectra as in Fig. A.1 but showing the slices at fixed azimuthal order m = 10. The red and blue solid lines denote the frequencies of the m = 10 HHS22 mode obtained from the linear calculation and measured by the observation, respectively. The cyan solid curve in panel b shows the Lorentzian fit for the spectra around the HHS22 mode power peak, and the cyan dashed line represents the fitted peak frequency. The same Lorentzian fitting is also performed for the columnar convective mode in panel (a) and shown by orange solid curve and dashed line.

Table A.1.

Properties of the HHS22 modes in our nonlinear rotating convection simulations for 1 ≤ m ≤ 15. The values in parentheses show the results from our linear model 3-sub (δcz = −5 × 10−7 and δsf = 10−3). The frequencies are measured in the Carrington frame. The peak frequencies and linewidths (FWHM) are obtained by Lorentzian fits. The quantity max(vhori) represents the maximum surface horizontal velocity amplitude of the mode eigenfunction extracted using the singular-value-decomposition (SVD).

In order to extract the spatial eigenfunctions of the HHS22 modes from the simulation data, we applied the singular-value decomposition (SVD) method to the l = m + 1 power spectrum of vθ for each m separately. For further details on the SVD eigenmode extraction, see Bekki et al. (2022b, § 3.2). Figure A.3 shows the extracted velocity eigenfunctions of the HHS22 mode at m = 10. Qualitatively, a good similarity can be confirmed in their overall spatial pattern of the velocity eigenfunctions with the linear model (Figs. 2c and 6c). However, there exist some differences such as the small-scale noise at high latitudes and the number of radial nodes at the equator. These can be attributed to the effects of stochastic turbulent convection and the difference in the latitudinal differential rotation profiles. In fact, the critical layers are formed closer to the equator (and even in the middle convection zone at high m) in the nonlinear simulations.

thumbnail Fig. A.3.

Velocity eigenfunctions of the m = 10 HHS22 mode extracted from the fully-nonlinear rotating convection simulation using SVD.

As reported in Table A.1, the typical velocity amplitudes of the HHS22 modes are vθ ≈ 0.5 − 3 m s−1 in the simulations, which are much weaker than those of the columnar convective modes but are comparable to those of the equatorial Rossby modes (Bekki et al. 2022b). It is speculated that the HHS22 modes are excited and damped by turbulent convection similarly to the equatorial Rossby modes (Philidet & Gizon 2023).

Lastly, we show in Fig. A.4 the radial profile of the superadiabaticity, δ(r), from the nonlinear rotating convection simulation. In the nonlinear simulation, the entropy stratification is subadiabatic near the base (r ≲ 0.75R) and superadiabatic in the upper bulk of the convection zone (0.75R ≲ r ≲ 0.95R). The formation of the weakly subadiabatic layer near the base is a direct consequence of the non-local convective heat transport (e.g., Käpylä & Rheinhardt 2017; Bekki et al. 2017; Karak et al. 2018). In the nonlinear simulation of Bekki et al. (2022b), the formation of this weakly subadiabatic layer is insignificant. However, recent magnetohydrodynamic simulations suggest that the subadiabatic layer tends to be enhanced and extended to upper convection zone as the numerical resolution is increased and the small-scale dynamo is better resolved (Hotta et al. 2022). In the real Sun, the mean entropy stratification is expected to be more subadiabatic than typically simulated due to the very efficient small-scale dynamo effect.

thumbnail Fig. A.4.

Radial profile of the superadiabaticity δ from the nonlinear rotating convection simulation. Black solid and dashed curves denote the area where the mean entropy stratification is superadiabatic (δ > 0) and subadiabatic (δ < 0), respectively.

All Tables

Table 1.

Summary of the linear analysis model setups.

Table A.1.

Properties of the HHS22 modes in our nonlinear rotating convection simulations for 1 ≤ m ≤ 15. The values in parentheses show the results from our linear model 3-sub (δcz = −5 × 10−7 and δsf = 10−3). The frequencies are measured in the Carrington frame. The peak frequencies and linewidths (FWHM) are obtained by Lorentzian fits. The quantity max(vhori) represents the maximum surface horizontal velocity amplitude of the mode eigenfunction extracted using the singular-value-decomposition (SVD).

All Figures

thumbnail Fig. 1.

Dispersion relation of the HHS22 modes obtained from the linear analysis for 1 ≤ m ≤ 19. Green points represent the results from our model 1 where we assume an incompressible (constant density) fluid and an uniform rotation. Blue points represent the results from our model 2 where the solar density stratification is included but the uniform rotation is still assumed. Orange points represent the results from our model 3 where both the density stratification and the differential rotation of the Sun are included. Red points also show the results from model 3 but with a weakly subadiabatic bulk convection zone (δcz = −5 × 10−7) and a strongly superadiabatic near-surface layer (δsf = 10−3). Lime and cyan solid curves show the results from triana2022 and from bhattacharya2022, respectively. For comparison, we also show the observed frequencies of the HHS22 modes reported in hanson2022 where the gray diamonds and squares denote the measurements with ring-diagram analysis (RDA) and mode-coupling analysis (MCA), respectively. All the frequencies are measured in the Carrington frame rotating at Ω0/2π = 456 nHz.

In the text
thumbnail Fig. 2.

Meridional cuts of the eigenfunctions of the m = 10 HHS22 mode obtained from the linear analysis. The velocity, pressure, and vorticity eigenfunctions are expressed as ν(r, θ)exp[i( − ωt)], p1(r, θ)exp[i( − ωt)], and ζ(r, θ)exp[i( − ωt)], and the solutions are shown in the meridional plane at t = 0 and ϕ = 0. The units of the colorbars are m s−1 for three velocity components, 104 dyn cm−2 for the pressure perturbation, and 10−8 s−1 for the vorticity components. The eigenfunctions are normalized such that the maximum of vθ is 1.0 m s−1. Panels (a–c) show the results from models 1–3, respectively. Note that the background entropy stratification is adiabatic in all cases. In panel (c), the black solid curves show the locations of the critical latitudes.

In the text
thumbnail Fig. 3.

(a) Root-mean-square (RMS) amplitudes of the velocity eigenfunctions of the HHS22 modes from the linear analysis, where the average is taken over the spherical surfaces. Dashed and solid curves show the results from model 2 (without differential rotation) and model 3 (with differential rotation). Different colors show different azimuthal orders m. The eigenfunctions are normalized such that the maximum horizontal velocity at r = 0.985 R is 1.0 m s−1. (b) Spectra of volume-integrated kinetic energy of the HHS22 modes. The blue and red points denote those from models 2 and 3, respectively.

In the text
thumbnail Fig. 4.

(a) Propagation frequencies and (b) linewidths of the HHS22 modes obtained from the linear calculations with two different values of turbulent viscosity, ν. Blue and red points show the results with spatially uniform viscosity of ν = 1011 cm2 s−1 and for 1012 cm2 s−1, respectively. The model 3 is used (with solar density stratification and solar differential rotation) and the background is assumed to be adiabatic (δ = 0). The observed values reported by hanson2022 are shown by gray diamonds.

In the text
thumbnail Fig. 5.

(a) Dispersion relation of the HHS22 modes computed for different values of superadiabaticity near the surface δsf. The bulk of the convection zone is assumed to be adiabatic, δcz = 0. (b–e) Eigenfunctions of the m = 10 HHS22 mode computed for δsf = 10−6 (left column) and for δsf = 10−3 (right column). Panels b–d show cuts of radial velocity vr (in m s−1), entropy perturbation s1 (in erg g−1 K−1), and z-vorticity ζz (in 10−8 s−1) in the equatorial plane (along the rotational axis) seen from the north pole. The black dashed curves denote the height r = 0.95 R, above which the strongly superadiabatic layer is located. (e) Mollweide projection of the radial vorticity eigenfunction ζr at the top boundary r = 0.985 R (in 10−8 s−1). All the eigenfunctions are normalized in the same way as in Fig. 2.

In the text
thumbnail Fig. 6.

(a) Radial profiles of the superadiabaticity δ(r) defined by Eq. (2) for different values of the bulk superadiabaticity δcz. The superadiabaticity in the near-surface layer is fixed to δsf = 10−3. (b) Dispersion relations of the HHS22 modes computed for different superadiabaticity profiles. The observed frequencies reported by hanson2022 are denoted by gray diamonds. (c) Meridional eigenfunctions from the case with weakly subadiabatic bulk convection zone (δcz = −5 × 10−7). Black solid and gray dashed curves show the locations of the critical latitudes (by differential rotation) and the turning surfaces (by subadiabatic stratification), respectively.

In the text
thumbnail Fig. A.1.

Power spectra from the nonlinear rotating convection simulation of (a) the l = m component of longitudinal velocity vϕ, (b) the l = m + 1 component of latitudinal velocity vθ, and (c) the l = m + 1 component of radial vorticity ζr near the surface r = 0.95R. The spectra are computed in the Carrington frame. The power is normalized at each m. Cyan diamonds denote the observed frequencies of the HHS22 modes hanson2022. Orange points show the dispersion relation of the HHS22 modes obtained from the linear calculation (model 3-sub with δcz = −5 × 10−7 and δsf = 10−3). The power ridge associated with the columnar convective modes is denoted by black arrows.

In the text
thumbnail Fig. A.2.

Same power spectra as in Fig. A.1 but showing the slices at fixed azimuthal order m = 10. The red and blue solid lines denote the frequencies of the m = 10 HHS22 mode obtained from the linear calculation and measured by the observation, respectively. The cyan solid curve in panel b shows the Lorentzian fit for the spectra around the HHS22 mode power peak, and the cyan dashed line represents the fitted peak frequency. The same Lorentzian fitting is also performed for the columnar convective mode in panel (a) and shown by orange solid curve and dashed line.

In the text
thumbnail Fig. A.3.

Velocity eigenfunctions of the m = 10 HHS22 mode extracted from the fully-nonlinear rotating convection simulation using SVD.

In the text
thumbnail Fig. A.4.

Radial profile of the superadiabaticity δ from the nonlinear rotating convection simulation. Black solid and dashed curves denote the area where the mean entropy stratification is superadiabatic (δ > 0) and subadiabatic (δ < 0), respectively.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.