Probing Depth Variations of Solar Inertial Modes through Normal Mode Coupling

Recently discovered inertial waves, observed on the solar surface, likely extend to the deeper layers of the Sun. Utilizing helioseismic techniques, we explore these motions, allowing us to discern inertial mode eigenfunctions in both radial and latitudinal orientations. We analyze 8 yr of space-based observations (2010–2017) taken by the Helioseismic and Magnetic Imager on board the Solar Dynamic Observatory using normal mode coupling. Couplings between the same and different-degree acoustic modes and different frequency bins are measured in order to capture the various length scales of the inertial modes. We detect inertial modes at high latitude with azimuthal order t = 1 and frequency ∼ −80 nHz, measured in a corotating frame with a rotation frequency of 453.1 nHz. This mode is present in the entire convection zone. The presence of Rossby modes may be seen down to a depth of ∼0.83R ⊙, and the Rossby signal is indistinguishable from noise below that depth for high azimuthal order. We find that the amplitudes of these modes increase with depth down to around 0.92R ⊙ and decrease below that depth. We find that the latitudinal eigenfunctions of Rossby modes deviate from sectoral spherical harmonics if we use a similar approach as adopted in earlier studies. We find that spatial leakage and even pure noise in the measurements of nonsectoral components can also explain the abovementioned characteristics of the latitudinal eigenfunctions. This realization underscores the necessity for careful interpretation when considering the latitudinal eigenfunctions of Rossby modes. Exploring the depth-dependent characteristics of these modes will enable us to capture interior dynamics distinctly, separate from p-mode seismology.


INTRODUCTION
The Sun supports a large number of inertial modes, e.g., equatorial Rossby modes, high-latitude inertial modes and critical latitudinal modes (Gizon et al. 2021).Previously, Löptien et al. (2018) detected equatorial Rossby modes (first observed in earth by Rossby 1939) in the Sun in which Coriolis force is the restoring mechanism (Papaloizou & Pringle 1978;Saio 1982;Provost et al. 1981).Since Coriolis force is the primary driver of these modes, their frequency is also of the same order as the rotation frequency of the Sun.

Equatorial Rossby modes
Equatorial Rossby modes have a well-defined dispersion relation based on their harmonic degree, s and azimuthal order, t, in a co-rotating frame.In this study, we opt for tracking at the surface equatorial rotation rate of Ω/(2π) = 453.1 nHz.We utilize harmonic degree and azimuthal order, denoted as (ℓ, m), to characterize solar acoustic modes, and employ (s, t) to represent inertial waves in subsequent sections.Observed Rossby modes closely follow this dispersion relation.Various methods of helioseismology have been used to confirm this finding, i.e., mode coupling (Hanasoge & Mandal 2019;Mandal & Hanasoge 2020;Mandal et al. 2021, hereafter HM19, MH20, MHG21 repectively), ringdiagram analysis (Hill 1988;Proxauf et al. 2020;Hanson et al. 2020), and time-distance helioseismology (Duvall et al. 1993;Liang et al. 2019).Although these prior studies were able to characterize the modes at the surface, the depth dependence has not been determined.Ring-diagram analysis and local correlation tracking are only sensitive to the near-surface layers.Proxauf et al. (2020) attempted to estimate the depth variations of Rossby modes from the surface down to 8 Mm using ring-diagram analysis, finding that amplitudes decrease with depth.Hathaway & Upton (2021), using supergranulation tracking, found mode amplitudes to be invariant with depth.Though Provost et al. (1981) derived the equations for the depth dependence of these modes for a uniformly rotating medium, solving that equation is not straightforward.Damiani et al. (2020) solved the equation derived by Provost et al. (1981) numerically for different polytropes in the inviscid limit.They found that only sectoral modes with no radial nodes may exist and that the radial dependencies of the modes scale as r t .Since the Sun is not a uniformly rotating inviscid polytrope, their results cannot be directly applied to the Sun.Moreover, careful consideration is required at the base of the convection zone, a sharp layer that separates the convectively stable radiative interior from the unstable zone above.Additionally, the Sun rotates differentially both in latitude and radius, whereas Equation 11 of Provost et al. (1981) was derived for a uniformly rotating medium.Since Rossby modes extend in both latitude and radius, differential rotation arguably (Dziembowski & Kosovichev 1987;Gizon et al. 2020) significantly influences these modes.Recently, Bekki et al. (2022); Triana et al. (2022); Bhattacharya et al. (2023) developed numerical methods which take into account solar stratification, turbulent diffusivities, differential rotation and a latitudinal entropy gradient to study inertial waves in the Sun.In particular, Bekki et al. (2022) assumed spatially uniform values for viscosity, thermal diffusivity and superadiabaticity to simplify the problem, which are not true for the Sun.Additionally, their computational domain excluded the tachocline and near-surface shear layer.How important those are in determining the radial eigenfunction still needs to be determined.Bekki et al. (2022); Bhattacharya et al. (2023) found that there exist different radial orders of Rossby modes, e.g., n = 0 and n = 1, which are associated with different frequencies at low azimuthal orders, but almost identical frequencies at high azimuthal order.In their study, Dikpati et al. (2022), based on 3D hydrodynamic shallow-water simulations of a rotating thin spherical shell, suggested that Rossby waves within supergranular layers may be excited through the inverse cascade of kinetic energy.Exploring the depth-dependent characteristics of Rossby waves and other inertial waves becomes significant in advancing our understanding of solar interior dynamics, with farreaching implications for predicting space-weather events.Since our current understanding of the radial dependence of Rossby waves is still nascent, we primarily focus on inferring their behavior from SDO/HMI data and do not attempt to explain it analytical or numerical means.We also characterize the latitudinal eigenfunctions of these modes in this work.

High latitude inertial modes
Analyzing 10 years of SDO/HMI data, Gizon et al. (2021) detected a swirling pattern with its maximum velocity occurring at high latitudes, concluding that it is a normal mode of the Sun.Similar features were detected earlier by Hathaway et al. (2013); Bogart et al. (2015); Hathaway & Upton (2021), although they were reported as giant-cell convection at the time.Though several analyses with different techniques have been used to confirm the existence of equatorial Rossby modes, been few studies have focused on characterizing inertial modes at high latitude, mainly because current instruments dedicated to helioseismology have limited coverage near the pole and systematical biases become significant near the limb (e.g., center-to-limb systematics in time-distance helioseismology Zhao et al. 2012).We apply the mode-coupling method in this work to detect this mode in the convection zone.
Characterizing the eigenfunctions of these inertial modes is important as they may provide information about the solar interior that is otherwise not possible to obtain by classical helioseismology.These modes are sensitive to turbulent viscosity and superadiabaticity in the convection zone (Bekki et al. 2022).Characterizing these modes will therefore enable us to map latitudinal entropy gradient, turbulent viscosity and superadiabaticity in the convection zone.Although MH20 performed several tests that involve accurately recovering the depth profiles of Rossby modes for synthetic measurements with added noise (also see in the appendix), we additionally verify our technique by recovering the s = 3 spatial scale of differential rotation from the surface down to the radiative interior and compare it with inferences obtained from global helioseismology.
2. DATA ANALYSIS High-resolution line-of-sight Doppler observations, Φ, taken by SDO/HMI, are analyzed here.In normal-mode coupling, we transform the data fully to the spectral domain, first to the spherical-harmonic domain to obtain Φ ℓm , where ℓ and m are the harmonic degree and azimuthal order of the p-mode respectively (for details, see Larson & Schou 2015).We then temporally Fourier transform these data to obtain Φ ω ℓm , where ω represents the temporal frequency.We obtain time series of these data products from the JSOC website1 .As p-modes are excited by turbulent convection, we expect the cross-correlation to contain power only for the quantity |Φ ω * ℓm Φ ω ℓm |.For simpler notation, we express Φ ω ℓm without explicitly showing its dependence on the radial order, n.Cross-spectral measurements such as ⟨Φ ω * ℓm Φ ω+σ ℓ+∆ℓm+t ⟩ in general will not carry any information because of the stochastic nature of the p-mode excitation, unless the difference in frequency, σ, harmonic degree, ∆ℓ and azimuthal order, and t match with the temporal and spatial scales of perturbations in the medium.This work examines the coupling, ⟨Φ ω * ℓm Φ ω+σ ℓ+∆ℓm+t ⟩ between modes of the same radial order, n.A general flow field, which we consider as a perturbation to our background medium, may be expressed in terms of poloidal and toroidal components where θ and ϕ are co-latitude and longitude respectively, and ∇ 1 = θ∂ θ + φ sin θ ∂ ϕ .The first two terms with u σ st and v σ st determine the poloidal component, and the final term, w σ st , captures toroidal motions.Mass conservation ties u σ st and v σ st , and as a consequence, Equation 2 has only two independent components.Here, we choose u σ st and w σ st as the independent parameters.The cross-spectral measurement, ⟨Φ ω * ℓm Φ ω+σ ℓ+∆ℓm+t ⟩, in presence of a flow as given by Equation 2 may be expressed as (Hanasoge et al. 2017) where the expression for H σ ℓℓ+∆ℓmt is given by Equation 5 of MHG21.K u st and K w st are sensitivity kernels for poloidal and toroidal components respectively.These kernels are dependent on the modes of interest (n, ℓ, m) and (n, ℓ + ∆ℓ, m + t).For the sake of simplicity in notation, we have refrained from listing all the dependencies in K u st and K w st , which may be simplified using their asymptotic expressions (Vorontsov 2011;Hanasoge 2018).
where the expression for γ ℓ+∆ℓsℓ tm is given by Equation 4 of MHG21.The formula for K nℓ (r) can be found in Hanasoge (2018) in terms of the eigenfunctions of the p-modes.The factor f is given by for odd s + ∆ℓ and the expression for g is for even s + ∆ℓ.Since Rossby and inertial modes at high latitude are mostly toroidal in nature, our choice of ∆ℓ is such that s + ∆ℓ is always odd.For toroidal flow, Equation 2 reduces to  (n, ℓ) for the frequency bin σ/(2π) = −80.4nHz, which also corresponds to the peak frequency of inertial waves at high latitude.Since B-coefficients show alternating sign with ℓ, we multiply it by (−1) ℓ and sum over all harmonic degrees, ℓ.We subsequently average this expression over all radial orders, n.The corresponding errors in the measured values (±1σ around the mean) is indicated by the pink shaded area.The inverted profile, w σ st , obtained using the measured B-coefficients has been used to in forward modeling and shown using green solid line.
Because of the selection criterion described subsequent to Equation 6, we may only choose values ∆ℓ = 0, 2, . . . to determine odd s, and ∆ℓ = 1, 3, . . . to detect even harmonic degrees.HM19 used ∆ℓ = 0 to determine the odd harmonic degrees of Rossby modes.Subsequently, MHG21 used ∆ℓ = 1, 3 to measure even harmonic degrees of Rossby modes.We neglect the poloidal component from Equation 3, We define the term in the first bracket as the B-coefficient (Hanasoge 2018) We see from Equation 10 that B coefficients depend on the modes (n, ℓ) and (n, ℓ + ∆ℓ) but not on azimuthal order, m and frequency, ω.Therefore, we estimate the B coefficients from Equation 9(Woodard 2016) using The summation applies over all m and ω to improve the signal-to-noise ratio.We ensure that the frequency interval over which the sum is carried out satisfies the criterion We analyze a span of 8 years of SDO/HMI data, covering the period from 2010 to 2017.Our focus involves investigating p and f-modes over mode-frequency and harmonic-degree ranges of ranges [1000, 5600] µHz and ℓ ∈ [10, 180], respectively.To comprehensively examine the low-frequency evolution of toroidal flows, we focus on the frequency range σ/(2π) ∈ [−500, 0] nHz when we track at the surface equatorial rotation rate of Ω/(2π) = 453.1 nHz.We solve the inverse problem set out in Equation 10, and thereby determine w σ st through the application of a regularized least-square technique, as outlined in MH20.By recovering w σ st as a function of radius, we gain insights into the depth-dependent behaviour of these modes.

High-latitude mode with azimuthal order t = 1
The retrograde inertial mode with azimuthal order t = 1, occurring at high latitudes (above 60 • ), exhibits the most pronounced amplitude, measuring approximately 10 m/s.In order to investigate this prominent mode, we choose ∆ℓ = 0 and an azimuthal order of t = 1 to compute B-coefficients using Equation 11.This choice facilitates the determination of w σ s,t=1 for all odd harmonic degrees, encompassing s values ranging from 1 to 20.Subsequently, we employ Equation 8 to evaluate the hemispherically symmetric component u θ and the antisymmetric component u ϕ .
In line with the terminology outlined by Gizon et al. (2021), these are referred to as u + θ and u − ϕ , respectively.For the scope of this study, we direct our attention to the symmetric mode.To visualize the signal within the measured Bcoefficients, we select the peak frequency σ = −80.4nHz, where the signal is most prominent.Notably, these measured B-coefficients exhibit an alternation in signs relative to harmonic degrees.Consequently, we depict the values (−1) ℓ B σ st in Figure 1.Importantly, the measured B-coefficient deviates from zero, displaying an increasing amplitude as ℓ grows.This trend implies the presence of a signal associated with the high-latitude mode.The amplification in intensity with higher values of ℓ is attributed to the ℓ 3/2 factor present in the sensitivity kernel equation (refer to Equation 8 in Hanasoge 2018).As long as the influence of the remaining terms in the sensitivity kernel remains subordinate, this increase in amplitude as ℓ rises is to be expected.We perform an inversion of the measured B-coefficients, yielding the radial profile w σ st (r).Subsequently, we proceed to forward model the inverted w σ st (r) and draw a comparison with the measured B-coefficients, as depicted in Figure 1. Figure 2 showcases the representation of u − ϕ on a spherical surface, revealing a discernible spiraling pattern at higher latitudes.A similar spiral pattern was earlier reported by Hathaway et al. (2013); Bogart et al. (2015); Hathaway & Upton (2021); Gizon et al. (2021).Additionally, we plot both u + θ and u − ϕ in Figure 3, presenting them in a frequency versus latitude diagram.The power of this mode is indeed focused near the pole.Upon normalizing the power at each latitude, we observe a consistent presence of power around the frequency of −80 nHz across all latitudes, as depicted in Figure 3.This widespread distribution suggests that this mode is global.Notably, the data used for the construction of Figure 3 spans a duration of four years, from 2014 to 2017.In an alternate assessment, we partition the 8 years of SDO/HMI data, spanning from 2010 to 2017, into two separate 4-year segments.In both scenarios, we observe that both u + θ and u − ϕ exhibit significant power centered around 80 nHz.As anticipated, the signal from the lower depths appears noisier due to increased noise, yet we are able to measure it all the way down to the base of the convection zone.To determine the central frequency and line-width of this mode, we create a power spectrum depicted in the lower panel of Figure 3.We then apply a fitting procedure using a Lorentzian function augmented by a constant background term, defined as follows: In this equation, A represents the amplitude, σ 0 denotes the central frequency, Γ signifies the full width at half maximum, and B stands for the background power.To execute the fitting process, we employ the curve fit module available in scipy.optimize.By applying this method to the first four years of SDO/HMI data (spanning 2010 − 2013), we obtain a central frequency of σ 0 /(2π) = −80.1 ± 2.1 nHz and a line-width of Γ/(2π) = 29.5 ± 6.8 nHz.Upon analyzing the data from SDO/HMI spanning the years 2014 to 2017, we determine the central frequency to be σ 0 /(2π) = −79.7 ± 0.5 nHz and the line-width as Γ/(2π) = 7.4 ± 2.0 nHz.Subsequently, by averaging the power spectra from these two instances and performing fitting, we establish the central frequency as σ 0 /(2π) = −81.0± 1.4 nHz and the line-width as Γ/(2π) = 18.6 ± 4.4 nHz.We depict the spatial representation of u − ϕ on a spherical surface in Figure 2, showcasing its depth variation.To assess the accuracy of our depth profile reconstruction, we subject it to a validation test against the differential rotation profile obtained through global helioseismology, as detailed in Appendix A.

Radial dependencies of Rossby modes
Different values of ∆ℓ can be used in Equation 11 to measure B-coefficients in order to image Rossby modes in depth (MHG21).We primarily consider ∆ℓ = 0 case, which has the highest signal-to-noise.By inverting Equation 10 using B-coefficients where s = t, we generate the radial profile w σ st (r).We illustrate the power |w σ st | 2 at various depths in Figure 4. We expect significant power above the background close to the theoretical dispersion relation of Rossby waves, Equation 1.As we go deeper, the background starts to dominate, and below 0.83R ⊙ , distinguishing signals for higher azimuthal orders, denoted as t, becomes challenging.However, traces of low azimuthal orders are still discernible in the plot; refer to Figure 4. Our findings align with a similar conclusion when we investigate the weighted sum of measured B-coefficients using the phase-filter method described in Appendix B.
In order to determine how these mode amplitudes vary, we plot w σ st as a function of radius for the frequency bins at which we observe the maximum power for each mode in Figure 4. Rossby modes have significant power even in the deeper layers.We utilize a total of 8 years of SDO/HMI data spanning 2010 to 2017, which we split into two distinct 4-year segments.We consider average B-coefficients from these two data-sets and perform the inversion.We find that the amplitudes of all modes with odd harmonic degrees starting from s = 1 to s = 11 first increase with depth down to around 0.92R ⊙ , decreasing at deeper layers.We perform similar experiments with the measured B-coefficients for Rossby modes as we did in Section 3.1 in order to test whether the observed B-coefficients deviate from the predicted B-coefficients when w σ st assumes an r s variation.We first select frequency bins where the observed B-coefficient attains maximum power for each mode, and then add the signed B-coefficients so that they contribute positively in the bottom panel of Figure 4. We compute theoretical B-coefficients for the synthetic w σ st profile that varies as r s .Bottom panel of Figure 4 shows that the synthetic B-coefficient deviates notably from the measured B-coefficients.st are plotted as functions of radius for two distinct time periods, encompassing odd-harmonic degrees ranging from s = 1 -11.The chosen frequency bin σ corresponds to the maximum power for each mode, and its specific value is provided in each panel (measured in nHz).To construct these plots, a dataset spanning a total of eight years from 2010 to 2017, acquired from SDO/HMI, was divided into two four-year chunks.The inversion process made use of the average B-coefficients.Horizontal error bars in each panel represent the widths of the averaging kernels at the corresponding depths.Error bars are calculated based on two, four-year data series.Bottom panel: sum of measured B-coefficients for Rossby modes, encompassing all oddharmonic degrees from s = 3 to s = 9.Frequency bins at which we find maximum power for each Rossby mode have been used in order to perform a similar analysis as in Figure 1.Reconstructed B-coefficients with the inverted profile are shown as the green solid line.If the depth profiles of Rossby modes had r s dependencies and maximum amplitudes is the same as observed at depth 0.99R⊙, our estimated B-coefficients would correspond to the blue dashed line.
Equation 8 implies that latitudinal eigenfunctions are determined by derivatives of spherical harmonics, r × ∇ h Y st .In our prior studies (HM19; MHG21), we have considered Rossby modes to be purely sectoral and only analyzed Bcoefficients for s = t.Sectoral modes peak at the equator and decay towards higher latitudes but do not change sign.Whereas observational studies by Löptien et al. (2018); Proxauf et al. (2020) found that surface eigenfunctions of these modes, which they label only with azimuthal number, are close to sectoral spherical harmonics but slightly different since they change sign at latitudes of ∼ ±30 • .At fixed azimuthal order, the eigenfunctions contain contributions from both sectoral and non-sectoral components (s ̸ = t).The length scales of a perturbation are labeled by both harmonic degree, s and azimuthal order, t (see Equation 9) in normal-mode coupling and we analyze the measured B-coefficients for s ̸ = t.The latitudinal eigenfunctions of these modes deviate from sectoral spherical harmonics if contributions from the s ̸ = t channels are sufficiently large.We analyze the B-coefficients for all possible cases where s ̸ = t, with the s = |t| + 2 of particular interest.We find significant power as seen in Figure 5 and the observed power at the spatial scale (s, t = s − 2) follows the relation −2Ω/(t + 1), which is also the frequency of the sectoral mode with azimuthal order t.In order to determine whether it is signal or spatial leakage from spatial scale (s, s) to (s + 2, s), we perform similar synthetic tests as discussed in MH20.We consider power only for sectoral modes (s, s) and estimate B-coefficients that leak from sectoral to non-sectoral components, as shown in Figure 5.It is seen that significant power from odd sectoral modes (s, s) leak into spatial wave numbers (s + 2, s) at the same frequency bins.A similar leakage phenomenon was also earlier reported by Woodard (2021).In appendix C, we explain why there would be a leakage similar to that observed in Figure 5.We plot the latitudinal eigenfunction of the t = 5 Rossby mode in Figure 5 for two cases.In one case, we consider it to be purely sectoral and in another, we consider contributions from the non-sectoral component (s ̸ = t).It is seen that the sectoral-case eigenfunction does not have zero crossings at higher latitudes in the former scenario whereas it has zero crossings in the latter case.Contamination by leakage as discussed above will make the interpretation of latitudinal eigenfunctions difficult in the present analysis.Following the approach used by Proxauf et al. (2020), we estimate the temporal covariance of vorticity at each azimuthal order between the equator and all other latitudes according to where ⟨.⟩ t denotes temporal averaging, ξ t is the centered vorticity ξ ′ t = ξ t − ⟨ξ t ⟩ t and T denotes time.We show the estimated function C(R, θ) in Figure 6 for each t.We also find similar zero crossings at high latitudes as observed in previous studies and use a Monte-Carlo simulation to estimate errors in the inferences.
We consider another possibility, namely that w σ t+2 t is entirely noise.For this case, only contributions from sectoral components w σ tt obtained from the inversion have been considered.We use a Gaussian noise model with zero mean and fixed standard deviation for w σ t+2 t .We vary the noise amplitude and study its impact on the inference of the eigenfunction, C t .We show our findings in Figure 6.Zero crossings at high latitudes (∼ ±30) are also found in this case.We analytically explain this behavior in Appendix C.1.As the noise level increases, the absolute value of the minimum of C t also increases and the latitude at which it crosses zero approaches the equator.This study demonstrates that care must be taken to interpret the latitudinal eigenfunctions obtained using mode coupling, especially if we apply the definition given in Equation 15.Our result shows that the presence of signal, leakage, and noise in non-sectoral components can all contribute to zero crossing.Hence, it is imperative to meticulously isolate each component and assess its individual contribution.This task warrants separate investigation, which we acknowledge as a crucial aspect for future studies.We underscored the importance of recognizing that leakage might be applicable to other methodologies as well, but the extent of its impact necessitates method-specific analysis.

DISCUSSION AND CONCLUSIONS
Through our mode-coupling analysis, we identify inertial modes at high latitudes characterized by an azimuthal order of t = 1 and a frequency of ∼ −80 nHz.This substantiates the discovery of high-latitude characteristics initially reported by Hathaway et al. (2013); Hathaway & Upton (2021), later confirmed as a normal mode by Gizon et al. (2021).Our findings indicate the penetration of this mode throughout the entire convection zone.To explore dynamics in the polar region effectively, local helioseismic techniques necessitate accurate measurements at high latitudes.Unfortunately, current instruments do not have good coverage at the high latitudes.Mode coupling makes use of global p modes, which are generally very reliable measurements, and can be used in principle to detect critical and other inertial modes reported by Gizon et al. (2021), which is part of set of future directions this work will take.We observe that our measured high latitude mode amplitude is approximately 4 m/s, whereas Gizon et al. (2021), Hathaway et al. (2013), andBogart et al. (2015) reported amplitudes on the order of ∼ 10 m/s, ∼ 20 m/s, and ∼ 0.5 m/s, respectively.The task of measuring at high latitudes poses challenges owing to the restricted coverage offered by current instruments, a limitation initially noted in global helioseismology.Additionally this discrepancy may stem from the varying sensitivities of different methods to high latitudes, resulting in different amplitude measurements.Notably, this specific mode's amplitude peaks exclusively near very high latitudes.While this study presents an initial effort towards understanding high-latitude inertial modes with azimuthal order m = 1, future work will address the challenges associated with measurements at high latitudes using normal-mode coupling.
Recently, Waidele & Zhao (2023) established a correlation between the power and frequency of Rossby modes and the solar cycle, employing time-distance helioseismology and ring-diagram analysis.This correlation can also be investigated through mode coupling analysis, as demonstrated by Hanasoge & Mandal (2019), who showed that even a two-year dataset from SDO/HMI can provide excellent signal-to-noise ratios for Rossby wave detection.Consequently, the mode-coupling technique emerges as a valuable tool for studying temporal variations in these waves and exploring potential correlations with solar cycle-related properties.In earlier work by MH20; MHG21, mode coupling was used to determine the eigenfrequencies and line-widths of Rossby modes.This work is aimed at investigating the depth structure of Rossby modes.We find that the amplitudes of these modes first increase with depth down to around 0.92R ⊙ .To ensure the reliability of these inferences, we validate the technique by performing an inversion for the s = 3 coefficient of differential rotation, as detailed in Appendix A, which is seen to match with global helioseismology results (Larson & Schou 2018).We do not intend to invert for spatial scales s > 3 since they are significantly affected by leakage.A more thorough analysis to remove this effect needs to be designed for this purpose (e.g., Kashyap et al. 2021).We also add another validation tests in Appendix D. We choose a profile characterized by a radial node, with modes exhibiting amplitudes of ∼ 1 m/s at the surface-a value in close proximity to the observed amplitude.In an effort to enhance the realism of the test, we perturbed the forward-modeled B-coefficients in alignment with the observed noise.We demonstrate that we are able to successfully recover the radial node, signifying that our analysis possesses the capability to reconstruct the input profile to a considerable extent.In Figure 4 we find that Rossby modes can be observed for r > 0.83R ⊙ .With increasing inference depth, the background becomes stronger, and at 0.83R ⊙ , it becomes impossible to discern these modes.We have used only 8-years of SDO/HMI data to determine the depth dependence of Rossby waves.Analyzing extensive time series data is essential for enhancing the signal-tonoise ratio in deeper layers.Our initial endeavor involves inferring inertial modes in these regions, with plans for a comprehensive analysis utilizing all available data series in the future.The error at the base of the convection zone is ∼ 15 cm/s, calculated for a single frequency bin.Considering error propagation across other frequency bins, that contribute to the background, is crucial for detecting these modes in deeper layers.We observe that background noise becomes notably prominent as we approach 0.83R ⊙ .If background noise remains low across these frequencies, it will allow us to probe deeper into the solar interior and identify signatures of modes in those deeper layers.Accounting for the potential temporal variability of Rossby waves, as demonstrated by Waidele & Zhao (2023), is crucial.Temporal averaging might therefore only provide an average profile for the specific period analyzed.
We also study the latitudinal eigenfunctions of these modes in Section 4. If modes are labeled only using azimuthal order, we expect contributions from both sectoral and non-sectoral components.We investigate the presence of signal in the non-sectoral components of these modes.We indeed find significant power in the t = s − 2 channel, one source for which is power leakage from sectoral components (s, s) to the non-sectoral components (s + 2, s).In a more comprehensive analysis, we will need to model and remove these leakage contributions in order to determine the latitudinal eigenfunction.Without mitigating these systematics, we estimate the latitudinal eigenfunctions (Left panel of Figure 6), inferring that they have a zero crossing at high latitudes around ±30 • .We also investigate how the conclusion will be affected if there were no leakage from sectoral to non-sectoral components but the non-sectoral component were purely noise.We consider Gaussian noise for non-sectoral components and we set the sectoral component to the inversion result.In this case also, we find similar zero crossings as shown in the right panel of Figure 6.This study demonstrates that care must be taken to interpret the latitudinal eigenfunctions obtained using mode coupling, especially if we apply the definition given in Equation 15.The amplitude of the Rossby mode is approximately 1 m/s, rendering their imaging at depth a challenging endeavor.Nevertheless, the high latitude inertial waves with azimuthal order t = 1 possess a greater strength compared to Rossby modes, facilitating a more manageable analysis.These high-latitude inertial modes extend throughout the entire convection zone.Looking ahead, we hold an optimistic outlook that by broadening the spectrum of couplings and refining the measurement methodology, we can achieve an enhanced signal-to-noise ratio.This, in turn, will enable us to image the deeper layers with increased accuracy.from mode coupling (solid lines) with the same spatial component from global helioseismology (dashed lines) for three different latitude ranges.The latitude range over which averaging is performed is stated in the panel.

2018
).We show our results in Figure 7.Despite using only one set of couplings and only one year's worth of data, it is possible to estimate w σ st for s = 3 from the surface down to the base of the convection zone.While we observe minor discrepancies, particularly in the vicinity of the surface beyond 0.97R ⊙ , the underlying reasons for these deviations remain unclear, necessitating further investigation in future studies.We do not attempt to estimate w σ st for s = 5 as we expect leakage from the more dominant s = 3 component.Therefore we need to take into account leakage if we want to analyze the s = 5 component, which we reserve for future work.The data analysis procedure used in previous studies with mode coupling are different from the data reduction technique applied here.This validation exercise gives us confidence in our attempt to infer radial eigenfunctions of inertial modes in the Sun.

B. ROSSBY MODES IN THE INTERIOR LAYERS VIA B-COEFFICIENTS
We estimate the sum of the signed B-coefficients, (−1) ℓ B σ st (n, ℓ, ℓ + ∆ℓ) over all ℓ used in the analysis and compute the sum of the squares of their absolute values over all radial orders, n, in order to estimate the quantity , where N σ st (n, ℓ, ℓ + ∆ℓ) corresponds to measurement noise.This quantity is only a function of two variables, σ and t.We apply a phase-speed filter, as these describe the regions through which acoustic modes traverse the interior.We have, ω/ℓ = c(r t )/r t , where r t corresponds to the inner turning points of the p-modes and c is the sound speed at that depth.We thus indirectly infer how deep the Rossby waves penetrate into the interior.We show our results in Figure 8, which characterizes different limits on ω/ℓ.We find that, when ω/ℓ > 45, it becomes difficult to discern signal from background (ω/ℓ = 45 corresponds to depth r t = 0.83R ⊙ .)

C. EFFECT OF LEAKAGE IN DETERMINING LATITUDINAL EIGENFUNCTIONS OF ROSSBY MODES
The observed line-of-sight Doppler velocity cannot be exactly decomposed into spherical harmonics due to our inability to observe the far side of the Sun.This results in a smearing in spectral space termed "leakage", resulting in cross-mode talk, i.e., (ℓ ′ , m ′ ) leaks into the mode of interest, (ℓ, m), given by where L ℓ ′ m ′ ℓm quantifies the geometric leakage from mode (ℓ, m) to ℓ ′ , m ′ .Because of this limitation, our measured B-coefficients, B σ st , with harmonic degree, s and azimuthal order, t, are contaminated by contributions from other spatial scales s ′ , t ′ , mathematically expressed as We select a radial profile for w σ st of Rossby modes and employ Equation D17 to estimate b σ st .In the absence of any leakage (i.e., when L ℓ ′ m ′ ℓm = δ ℓℓ ′ δ mm ′ ), it becomes evident from Equation D16 that B σ st simplifies to b σ st .To assess the impact of leakage on the modification of B σ st , we adopt a profile for w σ st given by Ar(r − 0.85)f (σ, σ s ), where f (σ, σ s ) represents a Lorentzian profile with a central frequency σ s derived from the dispersion relation of Rossby waves.The parameter A is selected to ensure that the amplitude of the modes closely approximates 1 m/s at the surface, aligning with observed values.Following forward modeling with this w σ st profile to derive b σ st , we substitute it into Equation D16 to determine B σ st .Our tests reveal that when s ′ ̸ = s, there is no contribution since the modes have distinct frequencies, σ.Instead, contribution arises from neighboring modes, (ℓ ′ , m ′ ).The left panel of Figure 9 illustrates the disparity between B σ st and b σ st .The figure highlights that the impact of leakage becomes more noticeable for higher ℓ, yet the changes in values are not overly significant, which bodes well for the analysis.In order to check whether we can recover the original profile, we invert the following Equation to B σ st .This substitution serves as an approximation that we employ in our analysis.To make it more realistic we perturb B σ st according to noise from observation.We present our inverted profile in Figure 9 and compare it with the original profile.We successfully recover the original profile with a moderate degree of accuracy, including the nodal point radius.To quantify the uncertainties in the inverted profile, we perturb B σ st according to observed noise, generating 100 realizations.For each realization, we perform inversion to obtain w σ st and calculate the variance in the radial profile.This process allows us to determine error bars in the inverted profile of w σ st .Based on our analysis, we conclude that the approximation B σ st ≈ b σ st is reasonably accurate.In the future, we will explore avenues for refining this approximation to achieve a more precise profile.

Figure 2 .
Figure2.Upper left panel: eigenfunctions of solar inertial modes of azimuthal order t = 1.A spiraling pattern associated with this mode is seen at high latitude.The upper middle panel displays the radial profile of u − ϕ , while the corresponding error is presented in the right panel.Bottom panel: Normalized power spectra of u + θ over the entire latitude range (solid black line) for different time periods, stated in the title of each panel.The power has been normalized in the frequency versus latitude diagram for each latitude.Subsequently, we compute the mean spectrum across latitudes to generate the normalized power versus frequency diagram.We fit a Lorentzian function with a constant background to the power spectrum.The fitted values of frequency, σ/(2π) and line-width, Γ/(2π) are mentioned in each panel.Fits to the spectrum are shown as red solid lines.

Figure 3 .
Figure3.Upper panels: power spectra of u + θ on the left and u − ϕ on the right.The mode amplitude exhibits significant strength at high latitudes, reaching magnitudes on the order of 4 m/s.Bottom Panel: Depiction of normalized symmetric latitudinal velocity, u + θ (left panel), and antisymmetric azimuthal velocity, u − ϕ (right panel), corresponding to a symmetric mode with azimuthal order t = 1.The depths at which inferences are made are specified in the titles of the respective panels.This illustration is based on a four-year dataset derived from SDO/HMI spanning from 2014 to 2017.

Figure 4 .
Figure 4. Upper left panel: normalized power |w σ st | 2 for odd harmonic degrees s as obtained from inversion at various depths, each indicated in the panel header.Additionally, the black-dashed line in the graph represents the theoretical dispersion relation of Rossby waves in a uniformly rotating medium, with a rotation frequency of Ω/(2π) = 453.1 nHz.Upper right panel: absolute values of w σst are plotted as functions of radius for two distinct time periods, encompassing odd-harmonic degrees ranging from s = 1 -11.The chosen frequency bin σ corresponds to the maximum power for each mode, and its specific value is provided in each panel (measured in nHz).To construct these plots, a dataset spanning a total of eight years from 2010 to 2017, acquired from SDO/HMI, was divided into two four-year chunks.The inversion process made use of the average B-coefficients.Horizontal error bars in each panel represent the widths of the averaging kernels at the corresponding depths.Error bars are calculated based on two, four-year data series.Bottom panel: sum of measured B-coefficients for Rossby modes, encompassing all oddharmonic degrees from s = 3 to s = 9.Frequency bins at which we find maximum power for each Rossby mode have been used in order to perform a similar analysis as in Figure1.Reconstructed B-coefficients with the inverted profile are shown as the green solid line.If the depth profiles of Rossby modes had r s dependencies and maximum amplitudes is the same as observed at depth 0.99R⊙, our estimated B-coefficients would correspond to the blue dashed line.

Figure 5 .Figure 6 .
Figure 5. Upper left panel: Normalized power of |w σ st | 2 from inversions when t = s − 2. s varies from 1 to 15, assuming only odd values.Upper right panel: We observe leakage in our synthetic test, transitioning from mode (s, s) to mode (s + 2, s).The black dashed line shows the function −2Ω/(s + 1) for mode (s + 2, s) to highlight that leakage occurs at the same frequency as the mode, (s, s).Lower left panel: Observed eigenfunction for the t = 5 Rossby mode if we consider the mode to be sectoral.Lower right panel: Observed eigenfunction for Rossby mode if we consider contributions from different s at fixed t = 5.Bottom panel: Real part of Ct(R, θ) at all azimuthal orders, mentioned in each panel.Corresponding errors (±1σ around mean) in the analysis are denoted by the red shaded area.

Figure 9 .
Figure9.Left panel: we plot the signed B-coefficient, B σ st (depicted as a black solid line with error bars), for harmonic degree s = 3 and t = 3, and for radial order n = 4 (as specified in the panel title), corresponding to the frequency bin where the (s=3, t=3) Rossby modes peak.Corresponding b σ st which is not affected by leakage is shown by red solid line.We perform forward modeling using the inverted profile of w σ st based on Equation D18.The resulting reconstructed B-coefficient is represented by the blue solid line.Right panel: Inverted profile of w σ st with harmonic degree, s = 3 is shown by black solid line.Corresponding error in the inverted profile (±1σ) is denoted by gray shaded area.
Equation D17, with the only difference being the substitution of the left-hand side from b σ st Figure 7. Left panel: We show measured B-coefficient values from one year of SDO/HMI data with odd and even harmonic degrees, ℓ separately.It is seen that the signs for these two cases are opposite to those indicated by the expression for sensitivity kernel.Right panel: We compare inversions for ws=3;t=0