Numerical study of non-toroidal inertial modes with l = m + 1 radial vorticity in the Sun’s convection zone

Various types of inertial modes have been observed and identiﬁed on the Sun, including the equatorial Rossby modes, critical-latitude modes, and high-latitude modes. Recent observations have further reported the detection of equatorially antisymmetric radial vorticity modes that propagate in a retrograde direction about three times faster than those of the equatorial Rossby modes, when seen in the corotating frame with the Sun. Here, we study the properties of these equatorially antisymmetric vorticity modes using a realistic linear model of the Sun’s convection zone. We ﬁnd that they are essentially non-toroidal, involving a substantial radial ﬂow at the equator. Thus, the background density stratiﬁcation plays a critical role in determining their dispersion relation. The solar di ﬀ erential rotation is also found to have a signiﬁcant impact by introducing the viscous critical layers and conﬁning the modes near the base of the convection zone. Furthermore, we ﬁnd that their propagation frequencies are strikingly sensitive to the background superadiabaticity, δ , because the buoyancy force acts as an additional restoring force for these non-toroidal modes. The observed frequencies are compatible with the linear model only when the bulk of the convection zone is weakly subadiabatic ( − 5 × 10 − 7 (cid:46) δ (cid:46) − 2 . 5 × 10 − 7 ). Our result is consistent with but tighter than the constraint independently derived in a previous study ( δ < 2 × 10 − 7 ), employing the high-latitude inertial mode. It is implied that, below the strongly superadiabatic near-surface layer, the bulk of the Sun’s convection zone might be much closer to adiabatic than typically assumed or it may even be weakly subadiabatic.


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. 2022b).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 HHS22) have reported to detect another family of modes near the surface of the Sun, i.e., 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.HHS22 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 call them HHS22 modes hereafter.
A possible identification of the HHS22 modes was first provided by Triana et al. (2022, hereafter TGBR22) using a linear model of rotating fluid in a spherical shell.They 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 TGBR22 was highly simplified, e.g., the model assumes an incompressible (uniform density) fluid and uniformly-rotating convection zone, which are not appropriate in the Sun.
A similar mode identification was later reported by Bhattacharya & Hanasoge (2023, hereafter BH23) in which the solarlike 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. 2022b;Fournier et al. 2022;Philidet & Gizon 2023;Bhattacharya et al. 2023).BH23 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 TGBR22.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. 2022b) which takes into account both the solar den- Notes.Model 1 assumes the incompressible fluid and uniform rotation, corresponding to that of TGBR22.In model 2, the solar-like density stratification is taken into account but the uniform rotation is still assumed, which is similar to that of BH23.In model 3, we consider both the background density stratification and the solar differential rotation determined by global helioseismology (Larson & Schou 2018).
sity stratification and the solar differential rotation determined by global helioseismology (Larson & Schou 2018).We will 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. (2022a) have mostly focused on the equatorial Rossby modes and the columnar convective (also known as thermal Rossby) modes.Matilsky et al. (2022) have implied that the equatorial Rossby modes might 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 Appendix of this paper, we will 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 § 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 § 3.1.The effects of turbulent viscosity is briefly discussed in § 3.2.We then show how the solar observations can be used to infer the superadiabaticity in the bulk of the convection zone in § 3.3.In Appendix A, we further report that the HHS22 modes are found to exist in the fully-nonlinear simulations of rotating convection.The conclusions are summarized in § 4.

Method: linear eigenmode analysis
We use the code developed by Bekki et al. (2022b), which numerically solves the linear eigenvalue problem for a rotating fluid in a spherical shell 0.71R ⊙ < r < 0.985R ⊙ .Here, R ⊙ is the solar radius.The model is hydrodynamic, i.e., the effects of magnetic field are ignored 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: -The eigenfunction of radial velocity v r is dominantly l = m in the middle convection zone, and has no radial node at the equator.-The eigenfunction of latitudinal velocity v θ is dominantly l = m + 1 at the surface, and has no radial node at low latitudes.-The eigenfunction of latitudinal 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. (2022b, § 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, 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 background stratification is included but the uniform rotation is still assumed.Orange points represent the results from our model 3 where both the solar stratification and the solar differential rotation 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 TGBR22 and from BH23, respectively.For comparison, we also show the observed frequencies of the HHS22 modes reported in HHS22 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.which were assumed in the previous study of TGBR22.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 BH23.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.

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 ν = 10 12 cm 2 s −1 .For simplicity, we assume that the background is purely adiabatic, i.e., there is no entropy variation neither in the radial nor the latitudinal directions.

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.0nHz.We find that the dispersion relation from the model 1 matches almost ex- actly to that of TGBR22; 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 BH23 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.
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 (HHS22).Those of model 1 (Fig. 2a) agree with the results of TGBR22 (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 shows 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 look similar to those of l = m + 1 classical Rossby modes).
To further investigate the consequences of the nontoroidalness, 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, i.e., 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. 2022b).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 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 β 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 TGBR22 is largely due to the unrealistic assumption of incompressible fluid in the Sun's convection zone, which ignores the compressional β-effect.

Effects of solar differential rotation
It is shown in Fig. 1 that the inclusion of solar differential rotation changes 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. 2022b;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, 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.

Dependence on turbulent viscosity
In this section, the dependence of the HHS22 modes on the turbulent viscosity ν is examined.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, ν = 10 11 and 10 12 cm 2 s −1 .For simplicity, the viscosity ν 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.In this section, we carry out a set of linear calculations with varying superadiabaticity δ.The most realistic model 3 is used with the turbulent viscousity of ν = 10 12 cm 2 s −1 .We use the following radial profiles of δ(r) which mimics a sharp transition from the close-to-adiabatic bulk convection zone to the strongly superadiabatic surface, 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 d sf = 0.015R ⊙ and r max = 0.985R ⊙ .Therefore, the strong superadiabaticities are localized in a thin layer near the top boundary of our numerical domain.

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.95R ⊙ ) is fixed to be purely adiabatic, i.e., δ cz = 0. Within the parameter range investigated here, the frequencies of the HHS22 modes are shown to be almost independent of δ sf .
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 .Figs. 5b-d show the eigenfunctions of radial velocity v r , z-vorticity ζ r , and entropy perturbation s 1 in the equatorial plane, respectively.The radial velocity v r is almost unchanged being concentrated near the base of the convection zone.However, the strong entropy perturbation s 1 is generated in the nearsurface superadiabatic layer as δ sf increases (Fig. 5c).The buoyancy force associated with this entropy perturbation drives the zvortical motion in this near-surface layer (Fig. 5d).Nonetheless, it is shown that the overall structure of the mode eigenfunctions in the bulk of the convection zone (below 0.95R ⊙ ) is only little affected by the inclusion of the strongly superadiabatic nearsurface 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.985R ⊙ up to O(10 −1 ) near 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.

Impact of superadiabaticity in the bulk convection zone
Next, we vary the superadiabaticity in the bulk of the convection zone δ cz within a range of ±10 −6 while fixing a value for the near-surface superadiabaticity δ sf = 10 −3 .The radial profiles of δ(r) used in this study are shown in Fig. 6a.
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, i.e., 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 ra- dial 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. (2022b, see § 5.), in which their prograde propagation frequencies become slower (faster) when δ > 0 (< 0).It is shown in Fig. 6b that the observed frequencies of the HHS22 modes 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 is estimated to be N/2π ≈ 240 − 340 nHz (with g ≈ 520 m s −1 and H p ≈ 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 gravitoinertial 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.In the bulk of the convection zone, therefore, 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.

Summary and Discussion
In this study, we investigated the properties of the l = m + 1 radial vorticity modes recently observed by HHS22 (we call them HHS22 modes in this paper).We used a linear eigenmode solver of the Sun's convection zone developed by Bekki et al. (2022b).Our model can successfully reproduce the previous results of TGBR22 and BH23 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.95R ⊙ ) to a strongly superadiabatic value near the top surface (0.95R ⊙ ≲ r ≤ 0.985R ⊙ ).Surprisingly, the strong superadiabaticities near the top surface δ sf do not affect the properties of the HHS22 modes such as their propagation frequencies and the radial vorticity eigenfunction at the surface (Fig. 5).In contrast, 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 exists 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, i.e., −5 × 10 −7 ≲ δ cz ≲ −2.5 × 10 −7 .This constraint is consistent with but tighter than the previous constraint derived by 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 by models with more realistic effects included.For instance, we ignored 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.985R ⊙ ) with a very strong superadiabaticity on the order of δ ≈ O(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ä et al. 2017;Hotta 2017;Bekki et al. 2017;Karak et al. 2018;Nelson et al. 2018;Hotta et al. 2022;Käpylä 2023) (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).which best-fits the observation (we will call 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. 2022a) 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.show the results from our linear eigenmode analysis from the model 3-sub (δ cz = −5 × 10 −7 and δ sf = 10 −3 ).The frequencies are measured in the Carrington frame.The peak frequencies and linewidths (full widths at half maxima) are obtained by Lorentzian fits.The quantity max(v hori ) represents the maximum surface horizontal velocity amplitude of the mode eigenfunction extracted from the nonlinear simulation using singular-value-decomposition (SVD).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.
In order to extract the spatial eigenfunctions of the HHS22 modes from the simulation data, we apply 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. (2022a, § 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 slight 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.
As reported in Table A.1, 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. 2022a).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 δ 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ä et al. 2017;Bekki et al. 2017;Karak et al. 2018).In the nonlinear simulation of Bekki et al. (2022a), the formation of this weakly subadiabatic layer is insignificant, leading to the radially-averaged superadiabaticity value of δ mean = 1.2 × 10 −7 .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 likely because of the very efficient small-scale dynamo effect.

Fig. 1 .
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 background stratification is included but the uniform rotation is still assumed.Orange points represent the results from our model 3 where both the solar stratification and the solar differential rotation 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 TGBR22 and from BH23, respectively.For comparison, we also show the observed frequencies of the HHS22 modes reported in HHS22 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.

Fig. 2 .
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 v(r, θ) exp [i(mϕ − ωt)], p 1 (r, θ) exp [i(mϕ − ωt)], and ζ(r, θ) exp [i(mϕ − ω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, 10 4 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 stratification is adiabatic in all cases.In panel (c), the black solid curves show the locations of the critical latitudes.

Fig. 3 .
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.985R ⊙ 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 model 2 and 3, respectively.

Fig. 4 .
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 ν = 10 11 cm 2 s −1 and for 10 12 cm 2 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 HHS22 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 ∇ = d ln T/d ln p.The superadiabaticity in the Sun's convection zone is typically estimated to be very small, δ ≈ O(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 O(10 −3 ) at r ≈ 0.99R ⊙ and even further up to O(10 −1 ) at the photosphere.

Fig. 5 .
Fig. 5. (a) Dispersion relations 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 v r (in m s −1 ), z-vorticity ζ z (in 10 −8 s −1 ), and entropy perturbation s 1 (in erg g −1 K −1 ) in the equatorial plane (along the rotational axis) seen from the north pole.The black dashed curves denote the height r = 0.95R ⊙ , above which the strongly superadiabatic layer is located.(e) Mollweide projection of the radial vorticity eigenfunction ζ r at the top boundary r = 0.985R ⊙ (in 10 −8 s −1 ).All the eigenfunctions are normalized in the same way as in Fig. 2.

Fig. 6 .
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 as shown in panel (a).The observed frequencies reported by HHS22 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.

Figures
Figures A.1a-c show the near-surface power spectra of the l = m component of the longitudinal velocity v ϕ , l = m + 1 component of the latitudinal velocity v θ , and l = m + 1 component of the 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 (HHS22) 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 will call 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. 2022a) 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.FigureA.2shows 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.2b.The measured peak frequencies and linewidths of the HHS22 modes in the nonlinear simula- Figures A.1a-c show the near-surface power spectra of the l = m component of the longitudinal velocity v ϕ , l = m + 1 component of the latitudinal velocity v θ , and l = m + 1 component of the 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 (HHS22) 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 will call 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. 2022a) 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.FigureA.2shows 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.2b.The measured peak frequencies and linewidths of the HHS22 modes in the nonlinear simula-

Fig. A. 1 .
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 (HHS22).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.

Fig. A. 2 .
Fig. A.2.The same power spectra as in Fig. A.1 but showing the slices at fixed azimuthal order m = 10.The red and blue solid lines deonte 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 green 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.

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

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 parenthesis