Time-distance helioseismology of solar Rossby waves

Context. Solar equatorial Rossby waves (r modes) have recently been observed in the horizontal flow field near the solar surface. This discovery used the techniques of granulation-tracking and ring-diagram analysis applied to six years of SDO/HMI data. Aims. Here we apply time-distance helioseismology to the SOHO/MDI and SDO/HMI data, which cover 21 years of observations from May 1996 to April 2017. The goal of this study is to provide an independent confirmation of the solar r modes over two solar cycles and in deeper layers of the Sun. Methods. We measure south-north helioseismic travel times along the equator, which are sensitive to subsurface north-south flows. To reduce noise, the travel times are averaged over travel distances from 6$^\circ$ to 30$^\circ$; the mean distance corresponds to a p-mode lower turning point of 0.91 $R_\odot$. The 21-year time series of travel-time measurements is split into three seven-year subsets and transformed to obtain power spectra in a co-rotating frame. Results. The power spectra all show peaks near the frequencies of the classical Rossby waves for azimuthal wavenumbers in the range $3 \leq m \leq 15$. The mode frequencies and linewidths of the modes with $m \leq 9$ are consistent with a previous study. Modes with $m \geq 10$ are shifted toward less negative frequencies by 10--20 nHz. While most of these modes have e-folding lifetimes on the order of a few months, the longest lived mode, $m=3$, has an e-folding lifetime of more than one year. The rms north-south flow speed associated with r modes along the equator is estimated to be in the range of 1--3 m s$^{-1}$, with the largest values for $m\sim10$. No evidence for the $m = 2$ mode is found in the power spectrum, implying that the rms velocity of this mode is below $\sim$0.5 m s$^{-1}$.


Introduction
Spheroidal oscillations associated with solar f , p, and g modes have been studied intensively in the past few decades (see, e.g., Basu 2016, for a recent review). Toroidal oscillations, including r modes, have been discussed in the literature as well (e.g., Papaloizou & Pringle 1978;Unno et al. 1989). For these modes the Coriolis force is the dominant restoring force; they are similar to the Rossby waves observed in the Earth's atmosphere and oceans. Provost et al. (1981), Smeyers et al. (1981), and Saio (1982) derived linearized equations for the analysis of r modes in uniformly and slowly rotating stars. Wolff & Blizard (1986) applied the analysis of r modes to the solar case and discussed the possibility of different radial orders. Later, Wolff (1998) studied the effect of differential rotation on r modes and suggested that nonsectoral modes might be suppressed while sectoral modes are the least affected because their motions are concentrated in the equatorial regions.
In the corotating frame and to lowest order in the mean solar rotation rate Ω, the frequencies of r modes are given by The movie associated to Fig. 1 is available at https://www. aanda.org where is the harmonic degree, and m is the azimuthal order. A consequence of this is that the horizontal flow field arising from the r-mode oscillations characterized by ( , m) drifts at a phase velocity of ω/m = −2Ω/[ ( +1)] < 0. The negative phase velocity means the drift direction is retrograde. Löptien et al. (2018) used six years of observations from the Helioseismic and Magnetic Imager on board the Solar Dynamical Observatory (SDO/HMI: Scherrer et al. 2012;Schou et al. 2012) to provide a direct and unambiguous detection of Rossby waves at the surface and in the outer 20 Mm of the Sun. Using both granulation-tracking (e.g., Löptien et al. 2016Löptien et al. , 2017 and ring-diagram analysis (e.g., Bogart et al. 2015), Löptien et al. (2018) find radial vorticity patterns along the equator which propagate retrograde in the corotating frame with a dispersion relation that is consistent with Eq. (1) for the m = case (sectoral modes). Löptien et al. (2018) find no evidence for the existence of the nonsectoral modes.
The r-mode velocity for sectoral modes at the equator is solely in the north-south direction (see Saio 1982, for sketches of the flow field of r modes). Figure 1 shows the flow field and vorticity of classical r-mode oscillations for the case of m = = 5. For this paper we have used time-distance helioseismology to measure subsurface flows in the meridional direction along the equator. The time-distance helioseismology pipeline of Liang et al. (2018) applies not only to the recent SDO/HMI data but also to SOHO/MDI data (Scherrer et al. 1995). The use of time-distance Fig. 1. Classical sectoral r mode with m = 5, seen from the equatorial plane of a uniformly rotating solar model. The three panels show the southward flow v θ , the prograde flow v φ , and the radial vorticity ζ r (from left to right) in the corotating frame. The color scale is the same for v θ and v φ with red positive and blue negative while the colors for the ζ r indicate radially outward (red) or inward (blue) directions. The black dotted lines represent constant longitudes fixed in the corotating frame. The green rectangle marks the equatorial area (±15 • ) in which we measure v θ in this paper. A movie showing the patterns propagating in the retrograde direction in the corotating frame is available online.
technique and the SOHO/MDI data would provide an independent confirmation of the findings of Löptien et al. (2018).

Time-distance analysis
We use medium-Dopplergrams taken by SOHO/MDI and SDO/HMI covering the period from May 1996 to April 2017. Background Doppler signals, such as solar rotation, are removed by subtracting a one-hour running mean for each pixel. A bandpass frequency filter is further applied to isolate the p modes within 2-5 mHz. Every 24-h time series of Dopplergrams are then projected into heliographic coordinates with a map scale of 0.6 • pixel −1 and tracked at the Carrington rotation rate. An improvement to the mapping procedure used in Liang et al. (2018) is made here by taking an error in the inclination angle of the solar rotation axis into account (see Appendix A.1).
We compute the cross-covariance function (CCF) between pairs of points arranged in the arc-to-arc geometry used by Liang et al. (2018) where two 30 • -wide concentric arcs, separated by an angular distance ∆, are aligned in the north-south direction. The CCFs between pairs of points on the opposing arcs are averaged and associated with the central point of the two arcs. The procedure is repeated for different central points located within ±15 • latitude at intervals of 0.6 • in longitude and latitude, and for the distance range ∆ = 6 • -30 • in steps of 0.6 • . Unlike the averaging scheme in Liang et al. (2018) where the CCFs were averaged over longitude and over days in each month, the daily CCFs are Gaussian smoothed with FWHM = 12 • in longitude and latitude, and subsampled along the equator at intervals of 10 • in longitude.
The south-north travel-time shift in this work is defined as δτ = τ s − τ n , where τ n and τ s are the northward and southward travel times of the first-skip wavelets in the CCF, respectively. We measure the south-north travel-time shifts from the spatially smoothed CCFs using the linearized one-parameter fitting method (Gizon & Birch 2002) as this algorithm is more robust to noise (Gizon & Birch 2004). More precisely, a 20-min interval around the first-skip wavelet in the CCF for positive time lag (i.e., the wavelet traveling in the northward direction) is selected as the reference function to compute the weight function derived by Gizon & Birch (2002). The weighted sum of the differences between the southward and northward wavelets in the CCF gives an estimate of the south-north travel-time shift. r/R Schematic plot of acoustic ray paths (red lines) that connect pairs of points across the equator in a meridional plane. The travel-time difference between the southward and northward propagating acoustic waves is sensitive to the north-south flow along the ray path. The larger the angular distance between the observation points at the surface, the deeper the lower turning point of the rays.
The measured south-north travel-time shifts for different travel distances are sensitive to subsurface north-south flows around the equator as depicted in Fig. 2. Data points within active regions are included in the averaging of CCFs in order to reduce the noise level. Also, the periods when the SOHO spacecraft was rotated by 180 • , which were not used in Liang et al. (2018), are used here for these alternate three-month gaps would result in strong leakage sidelobes in the Fourier domain.
To enhance the signal-to-noise ratio, we average the traveltime shifts over all travel distances with a weighting function that takes the noise correlation between different distances into account (see Appendix A.2). The weighted mean distance is about 14.6 • , for which the corresponding depth of the lower turning point from the ray approximation is about 63 Mm. The weighted average of the travel-time shifts suffers from a strong annual variation, mostly due to the center-to-limb effects (Zhao et al. 2012) 1 . To remove this time-varying background, we fit and subtract a periodic function from the measured travel-time shifts (see Appendix A.3). The Carrington coordinate system rotates at Ω cr /2π = 456.03 nHz in a sidereal frame. As viewed from the Earth, the Carrington coordinate system rotates at Ω cr −Ω ⊕ , where Ω ⊕ /2π = 31.69 nHz is the Earth's mean orbital frequency. The central times of the daily measurements are denoted by t j = t 0 + j∆t, where ∆t = 24 h and j ∈ {0, 1, . . . , 7669}, covering 21 years of data. At each time t j the Carrington longitudes of the travel-time measurements are given by  3. Section of the window function W(φ i j , t j ) in the frame that rotates at Ω eq /2π = 453.1 nHz. At time t j the longitudes of the travel times are denoted by φ i j (see main text). The window function is equal to one for i ∈ {−4, . . . 4}, one half at the boundaries i = ±5, and zero elsewhere. The window function is also zero for missing data. The temporal periodicity of the window is 2π/(Ω eq − Ω ⊕ ) = 27.46 days.
where i ∈ {−5, −4, . . . , 5} (11 measurements along the equator each day) and ∆ψ = 10 • is the spatial sampling rate. The longitude ψ 00 is the Carrington longitude of the central meridian as seen by the observer at time t 0 , which is the time when we start using MDI observations: t 0 = 1 May 1996 12:00:00_TAI. ( To compare with the results of Löptien et al. (2018), we transform to a frame that rotates at Ω eq /2π = 453.1 nHz, that is the solar surface equatorial rotation rate in the sidereal frame measured by global-mode helioseismology. The transformation between the two coordinate systems is given by where φ i j is the longitude measured in the new frame. Figure 3 is a schematic plot of the corresponding window function W(φ i j , t j ), which is equal to 1 for i ∈ {−4, . . . , 4}, tapered to 1/2 at the boundaries i = ±5 where noise is slightly higher, and zero elsewhere. The window function is also zero for missing data. We denote by δτ(φ i j , t j ) the travel times interpolated at the longitudes φ i j in the new frame. In practice, we implement the interpolation in the spatial Fourier domain by multiplying the spatial Fourier transform of the travel times by the phase factor exp[−im(Ω cr − Ω eq )(t j − t 0 )], where m is the azimuthal wavenumber.

Power spectra
The full data set δτ(φ i j , t j ) is divided into three consecutive periods of seven years. The starting times of these three periods, denoted by t (k) with k = {1, 2, 3}, are given by  The window functions for each 7-yr period is where C T (t j − t (k) ) is a taper time window that selects an observation period of T = N∆t ≈ 7 yr where N = 2556, starting from time t (k) (see Fig. 4).
The power spectrum of each 7-yr period is where the azimuthal order m is in the range |m| ≤ 18 (spatial Nyquist frequency). For each P (k) T (m, ω), the temporal Nyquist frequency is 5787 nHz and the frequency resolution is 4.5 nHz.
We denote the 2D Fourier transform of W (k) T (φ i j , t j ) by W (k) T (m, ω). Because the resulting P (k) T (m, ω) involves the convolution of W (k) T (m, ω), it is necessary to investigate the power spectrum of the window function, | W (k) T (m, ω)| 2 , from which the spectral leakage arises. We note that W (k) T (−m, −ω) = W (k) * T (m, ω) since the window function is real. The mean power spectrum of the window functions, 3 k=1 | W (k) T (m, ω)| 2 /3, consists of well-defined peaks, as can be seen in Fig. 5a (plot restricted to m ≥ 0). The maximum values of the peaks at m = 1 and m = 2 are 76% and 30% of the peak at the origin (m, ω) = (0, 0). The maximum power of the other peaks drops rapidly, 3% for m = 3 and 1% for m = 4. The decrease in the peak power with increasing m (as shown in Fig. 5c) corresponds to that of a square window covering 100 • longitudes. In Fig. 5a, adjacent peaks are shifted in frequency by (Ω eq − Ω ⊕ )/2π = 421.41 nHz due to the choice of coordinate frame. Since the power spectrum P (k) T involves a convolution in Fourier space of the solar data with W (k) T , the power spectrum of the data at (m, ω) leaks to multiple sidelobes at m + δm, ω − δm (Ω eq − Ω ⊕ ) , where |δm| = 1, 2, . . . Figure 5b shows a cut at m = 0 through the power spectrum of the window function. The width of the main lobe along the ω-axis is just one frequency bin and the frequency leaks to the neighboring frequency bins are nearly zero (two order of magnitude smaller than the central peak). Since we observe only a fraction of the Sun (see Fig. 3), there is spatial leakage to the neighboring m at frequency separations that are integer multiples of 421.41 nHz. Figure 5c shows these leaks, which implies a resolution in m of about 4. Figure 6 shows the P (k) T (m, ω) for the three different periods, as well as the mean spectrum T . The power spectrum is normalized to unity at (m, ω) = (0, 0). The maximum value of each peak is written right next to the peak. The frequency offset between the peaks with adjacent m is (Ω eq − Ω ⊕ )/2π = 421.41 nHz. Panel b: cut at m = 0 through the power spectrum from panel a. Panel c: maximum power of the peaks as a function of m, depicting the spatial leakage.
The power distribution in all three P (k) T (m, ω) peaks around the eigenfrequencies of r modes as described in Eq. (1) for the case of = m ≥ 3. However, there is no clear sign of = m ≤ 2. It is interesting to note that the mode amplitudes vary among the three P (k) T (m, ω), suggesting a possible temporal evolution of Rossby waves over the solar cycles. After the averaging, the ridge structure in the P(m, ω) along ω = −2Ω eq /(m + 1) becomes more prominent. The other two ridges of power seen in Fig. 6, which are separated from the central ridge by (δm, δω) = ±1, ∓(Ω eq − Ω ⊕ ) with a reduced amplitude (60∼80%), are due to the spectral leakage discussed previously. Excess power at low frequencies is present at low m, which also leaks and modulates the spectrum at frequencies between −300 and −500 nHz. It might be caused by large-scale convection or local flows around active regions that corotate with the Sun, and should not be mistaken for oscillation power. While the annual variation (31.7 nHz) in the background has been removed in the analysis, harmonics at ±63.4 nHz remain for m = 0 as the background fitting in Appendix A.3 only accounts for the first order term of B 0 -angle variation. We do not expect any aliasing at high m since we applied a Gaussian smoothing with FWHM = 12 • in longitude as mentioned in Sect. 2.1. Figure 7 shows the P(m, ω) in close-up around the frequency range of interest. The distribution of r-mode power seems to shift toward less negative frequencies by 10-20 nHz for m ≥ 10 compared to the mode frequencies measured by Löptien et al. (2018).
We note that our analysis is independent of the choice of reference frame. To illustrate this point, we present the power spectrum of data computed in the frame of the observer (as seen from Earth) in Appendix A.4. Compared to the power spectrum computed in the corotating frame, this power spectrum is shifted by m × 421.41 nHz for each m value. Figure A.4 displays a set of horizontal segments of size set by the resolution in m (≈ 4, see Fig. 5c). While the resolution in m is not one, one can easily retrieve the mode frequencies since the sectoral Rossby modes are well separated in frequency due to their long lifetimes.

Mode frequencies and linewidths of r modes
To quantify the line profiles of r modes, a model consisting of a Lorentzian function plus a constant background, is fit to the P(m, ω) at individual m. Here A is the maximum height of the Lorentzian function, ω m the mode frequency, Γ m the full width at half maximum, and B the constant background power.
To ensure that the Lorentzian fits are not affected by the spatial leaks from neighboring m (see Fig. 8, at negative frequencies, on the left) nor by the low-frequency power from active regions or convection (near zero frequency, on the right), we chose a fitting range [ν start , ν start + 300 nHz], with ν start given in Table 1.
The best fit is obtained by minimizing the sum with respect to the model parameters A, B, ω m , Γ m (Duvall & Harvey 1986). Although Eq. (9) was derived assuming that the quantity 2P(m, ω)/F m (ω) has a chi-squared distribution with two degrees of freedom for a single realization, we can minimize the same function for the case of multiple realizations (Anderson et al. 1990) as the P(m, ω) is an average over three power spectra. The minimization algorithm used is the downhill simplex method (e.g., Press et al. 1992, Sect. 10.4). Table 1 lists the resulting best fit parameters for 3 ≤ m ≤ 15. Figure 8 shows the fitted F m (ω) on top of the P(m, ω), along with the results of Löptien et al. (2018) for comparison. We have also computed the spectrum without splitting the 21-yr-long time series, the results of which appear to be similar to Fig. 8 but with better frequency resolution and higher noise level. The fitted mode frequencies and linewidths are remarkably consistent with that of Löptien et al. (2018) for 3 ≤ m ≤ 9. However, excess power on the right of the main peak is present for m ≥ 8, the extent of which for m ≥ 10 becomes so large that the overall profiles shift toward less negative frequencies.

Estimates of r-mode velocity
To obtain a rough estimate of the flow speed associated with a mode, first we compute in the ray approximation the traveltime shifts due to a prescribed toroidal flow field (see, e.g., Saio 1982). For the forward calculation, we choose a maximum horizontal flow speed at the surface of 2 m s −1 . We considered two toroidal flows: the first one is independent of depth and the second one decreases linearly with depth to vanish at 0.9 R . The forward-modeled travel-time shifts in the northsouth direction are Gaussian smoothed in longitude and latitude and averaged over travel distances in the same way as the measurements. The resulting travel-time shifts decrease with increasing m due to the smoothing in longitude. (1) with = m and Ω = Ω eq . The orange ellipse marks the excess low-frequency power at low m and the purple ellipse marks the spectral leakage from the low-frequency power. The gray scale is the same for the four panels and is shown in the color bar on the right. For clarity, the spectra are rebinned in frequency by a factor of three, such that the frequency resolution is 3/T = 13.6 nHz.  Table 1). For comparison, the mode frequencies measured by Löptien et al. (2018) are also indicated by green circles. to 0.05 s depending on m, from which we derive conversion constants of 15.6-41.7 m s −2 to convert from travel-time shifts to the surface flow speed. The conversion constants from the second model are in general larger by a factor of ∼1.6 for each mode.
Next, we obtain the rms travel-time shifts δτ rms from the fitted Lorentzian profile F m (ω) − B for each mode using Parseval's theorem. The effect of the incomplete data coverage and the spectral leakage is estimated by applying the same window functions and analysis to synthetic data, and is taken into account when computing the δτ rms .
Last, we use the above conversion constants to convert from δτ rms to the surface rms flow speed v rms . The conversion constants adopted are from the first model in which the flow is independent of depth, and thus the estimated v rms are conservative. The results for modes with 3 ≤ m ≤ 15 are listed in Table 1. For these modes, the surface velocity is on the order of 1-3 m s −1 , with larger values for 7 ≤ m ≤ 13. The modes with the lowest m values have velocities of order 1-2 m s −1 , as assumed in the RV study by Lanza et al. (2019).

Discussion
Using an independent helioseismic method and a different data set, we have confirmed the existence of the equatorial global Rossby waves reported by Löptien et al. (2018). We have extended the observations to deeper layers (down to ∼63 Mm) and to a total period of 21 years by combining SOHO and SDO observations. The power spectra obtained from three seven-year periods covering cycles 23 and 24 all show signatures of r modes for 3 ≤ m = ≤ 15. The measured mode frequencies and linewidths are generally consistent with the granulation-tracking results of Löptien et al. (2018). However, in our data, excess power is observed on the low-frequency side of the line profiles for m ≥ 10 which leads to a systematic shift of the fitted mode frequencies with respect to Löptien  In order to check if the spectral leakage comes into play, we applied the same window functions and Fourier analysis to synthetic data. The resulting spectra do not show apparent frequency shifts, implying that the excess power is not due to the leakage of neighboring modes. Because r-mode frequencies may vary over the solar cycle, we also examined separately the spectrum shown in Fig. 6c computed over a similar period as that used by Löptien et al. (2018). The above-mentioned excess power remains present, which means that this difference with Löptien et al. (2018) is not due to the use of MDI data in earlier periods.
Systematic effects associated with surface magnetic field cause solar cycle variations in the travel-time measurements (Liang & Chou 2015a). We tried excluding the data points inside the active regions from the averaging of CCFs; however, the Notes. Parameters ω m (mode frequency), Γ m (full width at half maximum), A (maximum height of the Lorentzian function), and B (background power) that give the best fit to the power spectrum P(m, ω) for each m in the range 3 ≤ m ≤ 15. The frequency ν start defines the frequency range [ν start , ν start + 300 nHz] for the Lorentzian fit. The 68% confidence interval on each parameter is indicated by the upper and lower bounds estimated from Monte Carlo simulations. The modes with m < 10 have a quality factor ω m /Γ m > 1 in the corotating frame (the frame for Rossby waves). However, the modes with m ≥ 10 have a quality factor less than or comparable to one. The e-folding lifetime, 2/Γ m , is given for each mode, as well as the oscillation period measured in the Sun's corotating frame (rotating at Ω eq /2π = 453.1 nHz). Negative oscillation periods indicates retrograde propagation. The rms south-north travel-time shift δτ rms for each m is obtained from the fitted Lorentzian F m (ω) − B using Parseval's theorem with the window function taken into account. The rms horizontal flow speed of the r modes along the equator in the near-surface layers, v rms , is converted from the δτ rms using conversion constants estimated from forward modeling (see text). A typical error in the v rms is about 0.5 m s −1 based on the error propagation. The maximum flow speed along the equator for each mode is related to the rms flow speed by v max = √ 2 v rms . Also listed are the mode frequencies from Löptien et al. (2018) for comparison. resulting power spectrum was dominated by noise as the masking procedure removed a considerable amount of pixels. Even if we could resolve the systematic effects of the surface magnetic field, the local flows surrounding the active regions may still enter the travel-time measurements (e.g., Gizon et al. 2001;Gizon 2004). We suspect that the excess low-frequency power at low m is caused by active regions and associated local flows since the rotation rate of the active regions is rather close to zero in our chosen rotation frame. Also the low-frequency power is stronger during the first (May 1996-April 2003 and third (May 2010-April 2017) periods when the solar activity is relatively higher. However, we cannot strictly exclude the possibility that the low-frequency power results from large-scale convection. Examination of the power spectra of travel-time measurements at different latitudes may help clarify this issue.
The errors in mode frequency measurements listed in Table 1 were estimated from Monte Carlo simulations. These errors can also be estimated using Eq. (2) from Libbrecht (1992) which relates the errors to the fitted parameters A, B, Γ m and the total observation time. The estimation of the errors on mode frequencies from the two methods are consistent. Our signal-to-noise ratio A/B is smaller than the granulation-tracking observations of Löptien et al. (2018). Twenty-one years of data analyzed using time-distance helioseismology give similar frequency error estimates as six years analyzed using granulation tracking.
It is interesting to compare the mode amplitudes measured here with those reported by Löptien et al. (2018). To this end, we first computed the radial vorticity of a prescribed r-mode toroidal flow field as used in the forward modeling, from which we obtained the ratio of rms velocity to rms vorticity in the latitude range between ±20 • for each mode. We then used the ratio to convert the rms radial vorticity reported by Löptien et al. (2018) to the rms velocity. The rms velocity estimated from their rms vorticity is on the order of 1-2 m s −1 and is consistent with the v rms in Table 1.
We do not see any evidence for the m = 2 mode (see Fig. 9). The travel-time measurements suffer from both timeindependent and annual variations for m ≤ 2 due to centerto-limb effects. The m = 2 mode is, however, expected to have a period of about 3/2 of the rotation period and should be cleanly separated from the center-to-limb systematics in the Fourier domain. To place an upper limit on the velocity of a possible m = 2 mode, we generated synthetic m = 2 sinusoidal travel-time shifts with frequencies in the range between −350 and −150 nHz. Three different amplitudes that correspond to rms velocity of 0.4, 0.5, and 0.6 m s −1 were implemented. These synthetic data were added into the measurements and Fourier analysis was applied with the window functions being taken into account. Figure 9 shows the resulting power levels. We find that an m = 2 sectoral mode with v rms < 0.5 m s −1 would be difficult to identify in the power spectrum.
The e-folding lifetimes of r modes, 2/Γ m , are on the order of a few months for most of the modes observed here. The longest lived modes (m = 3 and m = 5) have lifetimes of more than one year. Figure 10 shows a schematic excitation event (an exponentially decaying cosine) for the m = 5 r mode as seen in the and −150 nHz. The red dashed line is the threshold for 95% confidence level; that is, the noise in the background only has a 5% chance of being higher than this threshold for at least one frequency bin. We note that the spike around −295 nHz (on the right side of the blue line) is above the background but much lower than the 95% confidence level. The three green lines indicate the power that would correspond to a m = 2 sectoral mode with rms velocity of 0.4, 0.5, or 0.6 m s −1 .
corotating frame. The physical mechanisms responsible for the excitation, or damping, of the r modes are not yet known. While the r-mode oscillations in the near-surface layers have been remarkably determined using the granulation-tracking and ring-diagram analysis, the depth dependence of solar r modes throughout the convection zone is largely unknown. Timedistance helioseismology (Duvall et al. 1993) as used in this work has shown to be an important tool to study the r modes deep inside the Sun. More measurements of deep r-mode oscillations such as the one presented here will provide some insight into the nature of r modes and, hopefully, further constrain theories of solar r modes.  Table 1). We note that if the oscillation were seen in the Earth's frame, the observed oscillation period would be about 6 days for m = 5.   The zigzag might be due to the discreteness of the windows that isolate the single-skip wavelet in the cross-covariance function when fitting the travel-time shifts, particularly for short-distance cases in which the slope of the single-skip ridge in the time-distance diagram is steepest. Bottom: weighting function, w(∆), as a function of distance. The corresponding radii of the lower turning points from the ray approximation are indicated at the top. The weighted mean distance (vertical dashed line) is about 14.6 • , which corresponds to a lower turning point of ∼0.91 R .
and w(∆) as a function of distance. The trend of σ is consistent with Fig. 2 of Beck & Giles (2005) except that the noise level in their result is much lower probably owing to the use of phase-speed filter and different longitude and latitude ranges in the averaging. For ∆ < 10 • , the σ increases rapidly because of the low spatial resolution of medium-Dopplergrams and thus lack of information about the acoustic waves residing in the near-surface layers. Consequently, the weights for short distances are small. For ∆ > 10 • , the σ increases with increasing distance because of geometrical spreading and damping. Also, the noise correlation between different distances (off-diagonal elements of Λ) becomes greater when both ∆ and ∆ are large (Gizon & Birch 2004), which further reduces the weights for large-distance cases.
A.3. Removal of the background variation Following the notation in Sect. 2.1, at time t j , the Carrington longitudes of the travel-time measurements are denoted by ψ i j , and the central meridian as seen by the observer is denoted by ψ 0, j . Because the discussion here involves the center-to-limb effects, Jan '05 Jan '10 Jan '15 it is convenient to define the separation from the central meridian by ψ i j ≡ ψ i j − ψ 0, j . The weighted average of the travel-time shifts resulting from Appendix A.2 is expressed in this coordinate system and denoted by δτ (1) (ψ i j , t j ).
At each longitude ψ i j , a function representing the temporal variation of the background is fitted to and subtracted from the δτ (1) (ψ i j , t j ) for MDI and HMI data separately. Here t ref is taken from Appendix A.1. The a 0 (ψ i j ) represents the time-independent background, and a 1 (ψ i j ) represents the linear component, if any, of the slowly varying background. The a 2 (ψ i j ) accounts for an annual variation due to systematic errors such as the center-to-limb effects that vary with the B 0 angle (CRLT_OBS); in addition, the a 3 (ψ i j ) is added since the magnitude of center-to-limb effects was found to be time dependent (Liang & Chou 2015b;Liang et al. 2018). The a 4 (ψ i j ) accounts for another annual variation caused by the error in the inclination of the rotation axis, if not completely removed in Appendix A.1; the resulting error in the roll angle (i.e., Eq. (A.2)) may introduce a leakage of the solar rotation signal into the south-north travel-time shifts (Beck & Giles 2005). Figure A.2 shows the δτ (1) (ψ i j , t j ), the fitted background S (ψ i j , t j ), and the residuals δτ (2) (ψ i j , t j ) = δτ (1) (ψ i j , t j ) − S (ψ i j , t j ). (A.6) The annual variation is clearly seen in the δτ (1) (ψ i j , t j ), suggesting the necessity of dealing with this systematic effect. Also, the Jan '00 Jan '05 Jan '10 Jan '15 magnitude of δτ (1) (ψ i j , t j ) from HMI data at larger |ψ i j | is systematically greater than the rest. The fitted S (ψ i j , t j ) accounts for a large amount of the periodic variation of the background. After removing the fitted background, the resulting δτ (2) (ψ i j , t j ) from the two data sets are more consistent with each other, and are ready for the Fourier analysis. Figure A.3 shows the S (ψ i j , t j ) at selected longitudes as well as the fitted values of the parameters as a function of longitude for the two data sets. It is apparent that a 0 (ψ i j ) and a 2 (ψ i j ) (i.e., the center-to-limb effects) dominate the background variation while a 4 (ψ i j ) is nearly zero (0.008 ± 0.021 s for MDI and −0.015 ± 0.038 s for HMI within the longitude range |ψ i j | ≤ 30 • ) since we have applied a correction δI = −0.08 • to the inclination angle of solar rotation axis as discussed in Appendix A.1. Without this correction, the value of a 4 (ψ i j ) is around 0.1-0.2 s. We also tried δI = −0.1 • used by Larson & Schou (2015) who rounded the value −0.095 • ± 0.002 • determined by Beck & Giles (2005) up to −0.1 • ; however, the a 4 (ψ i j ) becomes negative in this case (−0.035 ± 0.019 s for MDI and −0.058 ± 0.038 s for HMI).
The background variation in the measurements from HMI is larger than that from MDI since the magnitude of center-tolimb effects for HMI is greater (Liang et al. 2017). The value of a 2 (ψ i j ) is expected to be negative because when the solar north pole is tilted toward the observer (CRLT_OBS > 0) the solar equator is on the southward side of the images where the center-tolimb variation in the south-north travel-time shifts is negative.
Furthermore, the background variation shows an east-west asymmetry which changes over time, especially in the measurements from MDI data. The east-west asymmetry in the traveltime measurements from MDI was reported by Giles (2000, Sect. 4.7.2). He attributed the asymmetry to the nonuniform A3, page 10 of 11 focus of the MDI camera, which changed at times during the mission. We note that the granulation-tracking measurements also suffer from systematic errors similar to the center-to-limb effects, termed the shrinking-Sun effect (Lisle et al. 2004). The shrinking-Sun effect not only depends on the distance to the disk center, but also exhibits an asymmetry in the east-west direction (Löptien et al. 2016(Löptien et al. , 2017. Löptien et al. (2016) suggested that the east-west asymmetry in the shrinking-Sun effect could be caused by the influence of the solar rotation on the observing height. If so, it would have a similar effect on the travel-time measurements since the magnitude of center-to-limb effects is strongly dependent upon the line-formation height of the observables (Zhao et al. 2012;Baldner & Schou 2012).
A.4. Power spectrum as seen from Earth where ω m is from Table 1 m(Ω eq − Ω ⊕ ) −2Ω eq /(m + 1) + m(Ω eq − Ω ⊕ ) 0 20 40 60 80 100 Fig. A.4. Power spectrum of south-north travel-time shifts computed in the Earth's frame. The blue line highlights the dispersion relation of sectoral modes of classical Rossby waves in this frame. The mode frequencies (red circles) are shifted by m × 421.41 nHz when measured in the Earth's frame. The low-frequency power from active regions or convection (near zero frequency at low m in Figs. 6 and 7) is also shifted to m × 421.41 nHz (green crosses) when measured in the Earth's frame.