Detection of Rossby modes with even azimuthal orders using helioseismic normal-mode coupling

Retrograde Rossby waves, measured to have significant amplitudes in the Sun, likely have notable implications for various solar phenomena. Rossby waves create small-amplitude, very-low frequency motions (on the order of the rotation rate and lower), which in turn shift the resonant frequencies and eigenfunctions of the acoustic modes of the Sun. The detection of even azimuthal orders Rossby modes using mode coupling presents additional challenges and prior work therefore only focused on odd orders. Here, we successfully extend the methodology to measure even azimuthal orders as well. We analyze 4 and 8 years of Helioseismic and Magnetic Imager (HMI) data and consider coupling between different-degree acoustic modes (of separations 1 and 3 in harmonic degree). The technique uses couplings between different frequency bins to capture the temporal variability of the Rossby modes. We observe significant power close to the theoretical dispersion relation for sectoral Rossby modes (where the azimuthal order is same as harmonic degree, s = |t|). Our results are consistent with prior measurements of Rossby modes with azimuthal orders over the range t = 4 to 16 with maximum power occurring at mode t = 8. The amplitudes of these modes vary from 1 to 2 m/s. We place an upper bound of 0.2 m/s on the sectoral t = 2 mode, which we do not detect in our measurements. This effort adds credence to the mode-coupling methodology in helioseismology


Introduction
Rossby waves are named after their discoverer, Carl-Gustaf Rossby, who first explained the largest scale oscillatory motions on Earth atmosphere (Rossby 1939) as arising from the conservation of potential vorticity. Chelton & Schlax (1996) observed these oscillations in the ocean by analyzing variations in sea-surface height from satellite data. These large-scale motions have implications on terrestrial weather and can influence convection and differential rotation in stars if they have sufficient amplitudes (Plaskett 1966). Rossby-like waves (also known as r-modes in the astrophysical context) can be sustained by any rotating spherical fluid body, such as the Sun, in which the Coriolis force acts as a restoring force (Papaloizou & Pringle 1978;Provost et al. 1981;Saio 1982). These waves have frequencies comparable to the rotation rate. For a uniformly rotating star, Rossby modes follow the dispersion relation (Provost et al. 1981) in a co-rotating frame, where s and t are the harmonic degree and azimuthal order of these modes respectively 1 . The availability of long-term high-resolution observational data of the Sun from Mount Wilson observatory, the Michelson Doppler Imager (MDI) on-board the Solar and Heliospheric Observatory (SOHO) and Global Oscillation Network Group (GONG) have inspired several authors (e.g. Kuhn et al. 2000;Ulrich 2001;Sturrock et al. 2015) to search for these oscillatory motions in the Sun. Sturrock et al. (2015) observed oscillations in solar radius measurements and attributed them due to r modes with t = 1. All these earlier studies lacked the measurement of the dispersion relation of these modes, critical to characterizing and understanding them. Löptien et al. (2018) were the first to measure the dispersion relation of these waves in the Sun using surface-granulation tracking methods and ringdiagram analysis with 6 years of SDO/HMI data and unambiguously detected modes with azimuthal order starting from t = 3 to t = 15. The surface eigenfunctions are close to A&A proofs: manuscript no. aanda sectoral spherical harmonics, although they are more peaked about the equator as a result of the effect of differential rotation . The detection of Rossby modes was later confirmed by several authors using different methods, e.g., by Liang et al. (2019) who used time-distance helioseismology, by Hanson et al. (2020) and Proxauf et al. (2020) who used ring-diagram analysis, and by Hanasoge & Mandal (2019) who applied normal-mode coupling (for a description of the method, see Woodard 1989;Lavely & Ritzwoller 1992;Roth & Stix 2003;Woodard et al. 2013). Interested readers are referred to Zaqarashvili et al. (2021) and references therein for a detailed review on Rossby waves in the Sun and other astrophysical problems. Although very powerful in its scope and inferential quality, the method of mode coupling is only slowly gaining traction in helioseismology (see Schad & Roth 2020;Hanson et al. 2021;Kashyap et al. 2021, for some other recent developments with this method); for instance, there is limited work on characterizing its limitations and developing mitigation strategies. The detection of Rossby modes by Hanasoge & Mandal (2019) using this method added an important milestone to this methodology. Mandal & Hanasoge (2020) compared properties of these modes, e.g., frequency and line width, with prior work by Löptien et al. (2018); Liang et al. (2019) and developed an approach to mitigate systematical errors that arise in this method, e.g., leakage in the spatial domain because of our limited vantage of the Sun. Hanasoge & Mandal (2019); Mandal & Hanasoge (2020) considered coupling between acoustic modes with same harmonic degrees. They detected Rossby modes only in the sectoral power spectra. Their measurement was not sensitive to even harmonic degrees which limited their detection to Rossby modes with odd harmonic degrees only. Here, we consider couplings between acoustic modes with different harmonic degrees and report the first detection of sectoral Rossby modes with even harmonic degrees using this method.

Data analysis
Line-of-sight Doppler velocity, Φ, observed by space-based, e.g., SOHO/MDI (Scherrer et al. 1995), SDO/HMI (Schou et al. 2012) and ground-based observatories, e.g., GONG are the main inputs to helioseismology. The data are transformed into the spherical-harmonic domain to obtain Φ ℓm , where (ℓ, m) are harmonic degree and azimuthal order of the p mode. As described earlier, we use s and t to indicate harmonic degree and azimuthal order of the Rossby waves to avoid confusion. A detailed discussion about how this data product is obtained from line-ofsight Doppler velocity data may be found in Larson & Schou (2015). We obtain these time series from the JSOC website 2 . We perform a temporal Fourier transform on the data to obtain Φ ω ℓm , where ω is the temporal frequency. The next step is to perform cross correlations across wave numbers, i.e., Φ ω * ℓm Φ ω+σ ℓ+∆ℓ m+t , where σ is varied to capture the time dependence of the perturbation and the difference between the two azimuthal order, t, captures the length scale of the perturbation. ∆ℓ is the difference between harmonic degrees of the two acoustic modes of interest. For the problem of Rossby waves, we know the dispersion relation analytically for a uniformly rotating fluid body with angular rotation rate Ω, which in the co-rotating frame (at the same rotational frequency), is given by Equation 1. The observed latitudinal eigenfunctions of these modes, labeled by only azimuthal order t, which may have contributions from sectoral (s = |t|) and non-sectoral modes (s |t|) (Löptien et al. 2018;Proxauf et al. 2020) peak at the equator and switch sign near 30 • latitude at the surface. Latitudinal eigenfunctions of these modes can not have any zero-crossings if these modes are purely sectoral. This can either be explained by considering presence of non-sectoral modes (Proxauf et al. 2020) or due to the effect of differential rotation on the Rossby mode eigenfunction . We in this work do not attempt to answer the above aspect of Rossby modes. Though normal mode coupling can easily distinguish between sectoral and non-sectoral Rossby modes, we in this work only focus on sectoral ones. In that case, the dispersion relation ( Equation 1 ) simplifies to . (2) We choose Ω/2π = 453.1 nHz, corresponding to the equatorial rotation rate of the Sun. We analyze the first 4 years of SDO/HMI data, from 2010 to 2014. We also analyze 8 years of SDO/HMI data from 2010 to 2018. The reason for choosing two different data sets is to demonstrate robustness of the method, i.e., that it can capture signals from a shorter time series. Though we have access to a few decades of high-resolution data for the Sun, we cannot say the same in the context of asteroseismology. Therefore, if a method works with a shorter time series of data, it can also be easily extended to asteroseismology. We analyze all the p-modes in the harmonic degree range ℓ ∈ [50, 180] and for all identified radial orders, n. We slightly modify our measurements for this particular work from that of Hanasoge & Mandal (2019); Mandal & Hanasoge (2020) as they considered correlations between acoustic modes with same harmonic degree which are only sensitive to odd harmonic degree of Rossby waves (see Equation 6 and 7 of Hanasoge 2018). For notational convenience, we use the same convention as in Hanasoge (2018) unless otherwise mentioned. Analyzing these correlation measurements at all temporal frequencies, ω, σ and for all pairs of acoustic modes with quantum numbers (n, ℓ, m) is very cumbersome. To simplify the analysis, we measure the B-coefficient, B σ st (Woodard 2016) which captures signal due to Rossby mode with harmonic degree, s and azimuthal order, t from couplings between acoustic modes with identical radial orders and ∆ℓ = 1, 3, which is defined as, where and where N ℓ is the mode normalization constant, L ℓ ′ m ′ ℓ m denotes leakage from mode (n, ℓ, m) to another mode (n, ℓ ′ , m ′ ) due to our limited vantage of the Sun (Schou & Brown 1994) and the expression of R ω ℓ m is (see Equation 11 of Hanasoge et al. 2017), where ω nℓm and Γ nℓ are the eigenfrequency and full width at half maximum of the mode, (n, ℓ, m). We analyze p-modes in the frequency range [1000,4600] µHz. In order to capture the temporal evolution of Rossby modes, we vary σ/2π ∈ [−350, 0.0] nHz, sufficient to capture the temporal evolution of all Rossby modes because its maximum frequency is σ/2π = −302 nHz for (s, t) = (2, 2) (we consider only even harmonic degrees in this work). For consistent analysis, we need to estimate noise in the measurement properly. We follow the analysis of Hanasoge (2018) to estimate the systematic noise in these measurements, and use Equation 3 to estimate variance. Taking a sum over all frequency bins ω in Equation 3 for all pairs of acoustic modes considered in the mode-coupling measurement is computationally very expensive. Therefore we restrict the sum to coupling between (n, ℓ, m) and (n, ℓ + ∆ℓ, m + t) modes in the frequency intervals The reason to subtract σ + tΩ in Equation 8 is to restrict the frequency interval to a full width around the resonance associated with (n, ℓ + ∆ℓ, m + t). We then connect these measured Bcoefficients with the perturbation in the medium, which in turn may be used for inference by performing an inversion as discussed in the next section.

Inversion
Rossby waves are described as a toroidal flow and may be expressed as Y s t is the spherical harmonic degree with degree s and azimuthal order t and ∇ 1 =θ∂ θ +φ sin θ ∂ φ whereθ andφ are the unit vectors along increasing co-latitude, θ and longitude, φ direction. The term w σ s t (r) captures the depth dependence of mode. We only focus on sectoral Rossby modes in this work as discussed in Section 2 and consider azimuthal order, t = s. The successful measurement of Rossby modes requires this function to attain significant amplitudes beyond the background power in the vicinity of the theoretical dispersion relation, i.e., Equation 2. If we use the flow profile from Equation 9 as a perturbation to the background model, the mode eigenfunction computed using the background model will be coupled, which, after some algebra, may be estimated as (see for Hanasoge 2018) where f ℓ ′ −ℓ,s K nℓ corresponds to the sensitivity kernel for the coupling between mode (n, ℓ, m) and (n, ℓ ′ , m + t). f ℓ ′ −ℓ,s is obtained from the asymptotic expression of the kernel (Vorontsov 2011) as the following for odd s + ℓ ′ − ℓ. This form of the kernel is valid when s ≪ ℓ or s ≪ ℓ ′ . The usefulness of choosing the asymptotic over the exact form is due to convenience as it separates out the dependencies on s (Rossby wave degree), ℓ ′ − ℓ (difference between harmonic degrees of the acoustic modes of interest) and radius into a product of a function f ℓ ′ −ℓ,s and the kernel K nℓ (r). We demonstrate the validity of this assumption in Figure 1 (also see Hanasoge et al. 2017). We choose two different cases, ℓ ′ − ℓ = 1 and ℓ ′ − ℓ = 3 and compare kernels obtained using the asymptotic expression with the corresponding exact forms. It is seen that these asymptotic kernel is a valid choice for our problem as the maximum harmonic degree of Rossby modes s (< 20) that we are considering is less than the minimum harmonic degree of the acoustic modes (ℓ = 50) used here. In our previous work (Hanasoge & Mandal 2019;Mandal & Hanasoge 2020), we were unable to infer evendegree Rossby modes as we only considered coupling between same-degree p-modes (see condition in Equation 11). Couplings between p modes with ∆ℓ = ℓ ′ − ℓ = 1 and 3 enable us to detect Rossby waves of even harmonic degrees. In order to estimate w σ st from the measured B-coefficients, B σ st , we need to perform inversions. We use the Regularized Least Squares (RLS) method discussed in Mandal & Hanasoge (2020) to invert Equation 10.

Results
The Rossby signal in the measured B-coefficients may be calculated by performing a weighted sum over the signed B-coefficients (−1) ℓ B σ st (n, ℓ, ℓ + ∆ℓ), over all harmonic degrees, ℓ. We then take the sum of the squares of their absolute values over all radial orders, n, to estimate the quantity n | ℓ (−1) ℓ B σ st (n, ℓ, ℓ+∆ℓ)/N σ st (n, ℓ, ℓ+∆ℓ)| 2 , where N σ st (n, ℓ, ℓ+ ∆ℓ) corresponds to measurement noise. If the measured Bcoefficients capture the signal properly, we would expect to see significant power close to the theoretical dispersion relation (Equation 2), which is indeed the case, as seen in Figure 2. We show results for ∆ℓ = 1 and ∆ℓ = 3 separately. The sectoral mode s = 2 does not appear when considering of ∆ℓ = 3 since finite couplings are only possible when s ≥ ∆ℓ. Therefore it can only be captured from measurements with ∆ℓ = 1.
We obtain w σ st (Equation 9) after performing inversions described in Section 2 and plot |w σ st | 2 at the depth r = 0.98R ⊙ in Figure 3. The spectra show significant power close to the theoretical dispersion relation, similar to Figure 2, for all even azimuthal orders t = 4 to 16. We do not detect power in the t = 2 sectoral mode, as seen in Figure 5. Power close to the theoretical frequency of this mode does not stand out from the background. Prior works by Löptien et al. (2018), Liang et al. (2019), Proxauf et al. (2020) and Hanson et al. (2020) have also reached a similar conclusion. In order to place constraints on the amplitude of this mode, we highlight the anticipated t = 2 Rossbymode frequency in Figure 5, where no peaks rise beyond the background. Based on the strength of the background power, the non-detection of this mode implies an upper amplitude limit of 0.2 m/s. We also see significant power for sectoral mode t = 16, though it is slightly shifted from the theoretical frequency. Indeed, we do not expect all the modes to exactly follow the analytical dispersion relation (Equation 2) since it is derived for the case of a uniformly rotating sphere, whereas the Sun shows radial and latitudinal differential rotation. Other factors, e.g., convection and magnetic field, are not accounted for in Equation 2. The detection of the s = 16 mode, were it to be independently confirmed, would highlight the resolving power of this technique, since previous studies have not observed it.
After obtaining the profile of w σ st from inversions, we fit a Lorentzian function with a constant background, to |w σ s t | 2 of these sectoral modes for each harmonic degree, s. Here A, σ 0 and Γ are amplitude, frequency and line-width of the mode respectively. D is the constant background power. We use power spectra obtained from 8 years of HMI data and the ∆ℓ = 3 case for fitting. We follow the approach of Anderson et al. (1990) and minimize the following function with respect to the model parameters, A, σ 0 , Γ and D. We apply the fmin subroutine, which uses the downhill simplex algorithm, implemented in the scipy.optimize function to minimize Equation 13. Values of these mode parameters obtained from the fitting are tabulated in Table 1. Fit spectra are shown in Figure 4. We compare the mode frequencies from our work with prior results of Löptien et al. (2018); Liang et al. (2019) in Table 1. We find that the Rossby-mode power peaks around t = 8, which is consistent with the previous studies by Löptien et al. (2018) and Liang et al. (2019).

Discussion and conclusions
Global-mode coupling, although a well known method, is still relatively unexplored in its applications to helioseismology. Given the challenges in understanding the data and given that it is a relatively novel measurement process, our prior work focused on same-degree couplings, which have the additional ben-efit of mitigating the leakage of background power into the resonances. However, this limited us to being able to only draw inferences of odd-harmonic-degree toroidal flows. We expand the set of inferences here by considering coupling between p modes with finite harmonic-degree separation, ∆ℓ = 1, 3. This allows us to detect sectoral Rossby modes with even azimuthal numbers in the Sun over the range t = 4 to 16. Similar to earlier findings (Löptien et al. 2018;Liang et al. 2019;Proxauf et al. 2020;Hanson et al. 2020), we also do not find evidence for the sectoral mode, t = 2. There is no current understanding why the t = 2 mode is undetectable. It is evident that power in this particular mode is very small or simply does not exist (see Figure 5), which is of interest in the context of understanding Rossby-mode excitation and dissipation. We also see evidence of sectoral-mode power beyond azimuthal number t = 15. Though we consider two cases ∆ℓ = 1 and 3 in this work, other cases, e.g., ∆ℓ > 3 may also be considered for this study, keeping in mind one caveat: modes with harmonic degree s < ∆ℓ will not be captured due to the selection rule discussed in Section 2. We find that signal to noise ratio is better for ∆ℓ = 3 than ∆ℓ = 1 (see Figure 3). We will investigate how signal to noise ratio varies for different ∆ℓ in future work and whether it can be improved by combining all the seismic data together in the analysis instead of separately, as done here.
We present our results only for near-surface layers in this work. However, because we analyze couplings between global modes, we are able to infer w σ st as a function of radius after inversion. The inference of the radial dependence of these modes We estimate the weighted sum n | ℓ (−1) ℓ B σ st (n, ℓ, ℓ + ∆ℓ)/N σ st (n, ℓ, ℓ + ∆ℓ)| 2 using the observed B-coefficients obtained after considering the coupling between acoustic modes with ∆ℓ = 1 (left panel) and 3 (right). Upper and lower panels show results from analyses of 4 and 8 years of SDO/HMI data, respectively. Power for each azimuthal order, t, has been normalized. The black dashed line is the theoretical dispersion relation for a uniformly rotating fluid with rotation frequency Ω/(2π) = 453 nHz. Modes with even azimuthal orders ranging from t = 2 to 16 are shown here.
at both odd and even azimuthal numbers is deferred to future work; this may shed light on the excitation mechanism of these modes. Although Mandal & Hanasoge (2020) have shown that spatial leakage does not affect inferences of the radial dependencies of these modes, care must still be taken as the observed data might contain other systematical errors, which may influence the inferences. This work, taken together with our previous work (Hanasoge & Mandal 2019), highlights mode-coupling as very useful in determining length and time scales of multi-scale dynamics in the Sun.  Notes. Analysis of first 8 years of SDO/HMI data. Fitted values of Mode frequency, σ 0 , amplitude √ A, full width at half maximum Γ to the observed B-coefficient spectra are listed in the co-rotating frame. MC stands for mode coupling here. The theoretical frequencies of modes in a co-rotating frame with tracking frequency 453 nHz are given in the second column in the table. Observed frequencies from two other studies, Löptien et al. (2018) and Liang et al. (2019) are listed in the fourth and fifth column in the table respectively. We plot the fitted spectrum in