© Z.-C. Liang et al. 2019

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

$$\begin{array}{c}\hfill \omega \approx -\frac{2m\mathrm{\Omega }}{\ell (\ell +1)}\text{for}m>0,\end{array}$$(1)

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. 2016, 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 technique and the SOHO/MDI data would provide an independent confirmation of the findings of Löptien et al. (2018).

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.

2. Data analysis and results

2.1. 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 band-pass 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⁻¹ 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. 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.

Fig. 2. 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.

To enhance the signal-to-noise ratio, we average the travel-time 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)¹. 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₀ + 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

$$\begin{array}{c}\hfill {\psi }_{\mathit{ij}}={\psi }_{00}+i\mathrm{\Delta }\psi -({\mathrm{\Omega }}_{\mathrm{cr}}-{\mathrm{\Omega }}_{\oplus })(t_j-t_0),\end{array}$$(2)

where i ∈ { − 5, −4, …, 5} (11 measurements along the equator each day) and Δψ = 10° is the spatial sampling rate. The longitude ψ₀₀ is the Carrington longitude of the central meridian as seen by the observer at time t₀, which is the time when we start using MDI observations:

$$\begin{array}{c}\hfill t_0=1\text{May}1996\text{12:00:00}\_\text{TAI}.\end{array}$$(3)

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

$$\begin{array}{c}\hfill {\varphi }_{\mathit{ij}}={\psi }_{\mathit{ij}}+({\mathrm{\Omega }}_{\mathrm{cr}}-{\mathrm{\Omega }}_{\mathrm{eq}})(t_j-t_0),\end{array}$$(4)

where ϕ_ij is the longitude measured in the new frame. Figure 3 is a schematic plot of the corresponding window function W(ϕ_ij, 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.

Fig. 3. Section of the window function W(ϕij, tj) in the frame that rotates at Ωeq/2π = 453.1 nHz. At time tj the longitudes of the travel times are denoted by ϕij (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.

We denote by δτ(ϕ_ij, t_j) the travel times interpolated at the longitudes ϕ_ij 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₀)], where m is the azimuthal wavenumber.

2.2. Power spectra

The full data set δτ(ϕ_ij, 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

$$\begin{array}{cc}\hfill t^{(1)}& =t_0,\hfill \\ \hfill t^{(2)}& =1\text{May}2003\text{12:00:00}\_\text{TAI},\hfill \\ \hfill t^{(3)}& =1\text{May}2010\text{12:00:00}\_\text{TAI}.\hfill \end{array}$$

The window functions for each 7-yr period is

$$\begin{array}{c}\hfill W_{\mathrm{T}}^{(k)}({\varphi }_{\mathit{ij}},t_j)=C_{\mathrm{T}}(t_j-t^{(k)})W({\varphi }_{\mathit{ij}},t_j),k=\{1,2,3\},\end{array}$$(5)

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).

Fig. 4. Window function CT(t) (thick solid line) that selects an observation period of T = 2556Δt ≈ 7 yr, tapered with one-year raised cosines at both ends.

The power spectrum of each 7-yr period is

$$\begin{array}{c}\hfill P_{\mathrm{T}}^{(k)}(m,\omega )=\frac{1}{N}{{\displaystyle \left|\sum _{i,j}{W}_{\mathrm{T}}^{(k)}({\varphi }_{\mathit{ij}},t_j)\delta \tau ({\varphi }_{\mathit{ij}},t_j)e^{-\mathrm{i}m{\varphi }_{\mathit{ij}}+\mathrm{i}\omega t_j}\right|}}^2,\end{array}$$(6)

where the azimuthal order m is in the range |m| ≤ 18 (spatial Nyquist frequency). For each $P_{\mathrm{T}}^{(k)}(m,\omega )$, the temporal Nyquist frequency is 5787 nHz and the frequency resolution is 4.5 nHz.

We denote the 2D Fourier transform of $W_{\mathrm{T}}^{(k)}({\varphi }_{\mathit{ij}},t_j)$ by ${\hat{W}}_{\mathrm{T}}^{(k)}(m,\omega )$. Because the resulting $P_{\mathrm{T}}^{(k)}(m,\omega )$ involves the convolution of ${\hat{W}}_{\mathrm{T}}^{(k)}(m,\omega )$, it is necessary to investigate the power spectrum of the window function, $|{\hat{W}}_{\mathrm{T}}^{(k)}{(m,\omega )|}^2$, from which the spectral leakage arises. We note that ${\hat{W}}_{\mathrm{T}}^{(k)}(-m,-\omega )={\hat{W}}_{\mathrm{T}}^{(k)\ast }(m,\omega )$ since the window function is real.

The mean power spectrum of the window functions, $\sum _{k=1}^3{|{\hat{W}}_{\mathrm{T}}^{(k)}(m,\omega )|}^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_{\mathrm{T}}^{(k)}$ involves a convolution in Fourier space of the solar data with ${\hat{W}}_{\mathrm{T}}^{(k)}$, 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.

Fig. 5. Panel a: mean power spectrum of the window functions $W_{\mathrm{T}}^{(k)}$. 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.

Figure 6 shows the $P_{\mathrm{T}}^{(k)}(m,\omega )$ for the three different periods, as well as the mean spectrum

$$\begin{array}{c}\hfill \overline{P}(m,\omega )=\frac{1}{3}(P_{\mathrm{T}}^{(1)}(m,\omega )+P_{\mathrm{T}}^{(2)}(m,\omega )+P_{\mathrm{T}}^{(3)}(m,\omega )).\end{array}$$(7)

Fig. 6. Power spectra of south-north travel-time shifts measured in the frame rotating at Ωeq/2π = 453.1 nHz. Panels a–c: $P_{\mathrm{T}}^{(k)}(m,\omega )$ from three time periods. Panel d: mean spectrum $\overline{P}(m,\omega )$ from an average of the three $P_{\mathrm{T}}^{(k)}(m,\omega )$ (seven years each). The blue lines highlight the dispersion relation of the classical Rossby waves described in Eq. (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.

The power distribution in all three $P_{\mathrm{T}}^{(k)}(m,\omega )$ 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_{\mathrm{T}}^{(k)}(m,\omega )$, suggesting a possible temporal evolution of Rossby waves over the solar cycles. After the averaging, the ridge structure in the $\overline{P}(m,\omega )$ 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₀-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 $\overline{P}(m,\omega )$ 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).

Fig. 7. Enlargement of Fig. 6d around the frequency range of interest (see Fig. 8 for the line profiles of individual modes in the range 3 ≤ m ≤ 14). The red circles show the mode frequencies estimated from Lorentzian fits (see Sect. 2.3). The errors in the mode frequencies are given roughly by the size of the red circles (see errors in Table 1). For comparison, the mode frequencies measured by Löptien et al. (2018) are also indicated by green circles.

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.

2.3. 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,

$$\begin{array}{c}\hfill F_{\mathrm{m}}(\omega )=\frac{A}{1+{[(\omega -{\omega }_{\mathrm{m}})/({\mathrm{\Gamma }}_{\mathrm{m}}/2)]}^2}+B,\end{array}$$(8)

is fit to the $\overline{P}(m,\omega )$ 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.

Fig. 8. Power spectra of south-north travel-time shifts $\overline{P}(m,\omega )$ for modes in the range 3 ≤ m ≤ 14 (black curves, with the frequency resolution 1/T = 4.5 nHz). The travel times are measured in the frame rotating at equatorial rotation rate Ωeq/2π = 453.1 nHz. The red lines are the fits Fm(ω) given by Eq. (8). The green vertical lines indicate the mode frequencies from Löptien et al. (2018). The orange arrows mark the excess low-frequency power at low m that might be caused by active regions or large-scale convection. The purple arrows mark the leaks from the m − 1 r-mode power and the low-frequency power at m − 1. Each power spectrum is accompanied by a plot of the ratio $\overline{P}(m,\omega )/F_{\mathrm{m}}(\omega )$ in gray in the lower panel, which is expected to have a mean of unity (dashed line) and a constant variance if the fit is not biased.

The best fit is obtained by minimizing the sum

$$\begin{array}{c}\hfill {\displaystyle \sum _{\omega }}{\displaystyle (ln{F}_{\mathrm{m}}(\omega )+\frac{\overline{P}(m,\omega )}{F_{\mathrm{m}}(\omega )})},\end{array}$$(9)

with respect to the model parameters A, B, ω_m, Γ_m (Duvall & Harvey 1986). Although Eq. (9) was derived assuming that the quantity $2\overline{P}(m,\omega )/F_{\mathrm{m}}(\omega )$ 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 $\overline{P}(m,\omega )$ 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 $\overline{P}(m,\omega )$, 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.

2.4. 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 travel-time 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⁻¹. 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 north-south 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. The maximum travel-time shifts from the first flow model range from 0.13 s to 0.05 s depending on m, from which we derive conversion constants of 15.6–41.7 m s⁻² 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⁻¹, with larger values for 7 ≤ m ≤ 13. The modes with the lowest m values have velocities of order 1–2 m s⁻¹, as assumed in the RV study by Lanza et al. (2019).

3. 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 et al. (2018).

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 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⁻¹ 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 time-independent and annual variations for m ≤ 2 due to center-to-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⁻¹ 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⁻¹ would be difficult to identify in the power spectrum.

Fig. 9. Power spectrum of south-north travel-time shifts $\overline{P}(m,\omega )$ for m = 2 (black solid line). The blue vertical line indicates the frequency of the classical m = 2 sectoral r mode. The red solid line is the background B estimated by a fit to the power in the frequency range between −350 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.

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 corotating frame. The physical mechanisms responsible for the excitation, or damping, of the r modes are not yet known.

Fig. 10. Schematic damped oscillation of a sectoral r mode with m = 5, seen at a fixed longitude in the Sun’s corotating frame. The value of the e-folding lifetime (vertical line) is as observed (see 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.

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. Time-distance 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.

Movie

Movie of Fig. 1 Access here

Acknowledgments

We thank B. Löptien, B. Proxauf and J. Schou for useful discussions. We also thank the referee for comments that helped improve the manuscript. The HMI data used are courtesy of NASA/SDO and the HMI science team. SOHO is a project of international cooperation between ESA and NASA. The data were processed at the German Data Center for SDO (GDC-SDO), funded by the German Aerospace Center (DLR). L.G. acknowledges partial support from the NYU Abu Dhabi Center for Space Science under grant G1502. We used the workflow management system Pegasus funded by The National Science Foundation under OCI SI2-SSI program grant #1148515 and the OCI SDCI program grant #0722019.

References

Appendix A: Supplementary m aterial

A.1. Inclination of the rotation axis

Because the inclination angle of the Sun’s rotation axis determined by Carrington (1863) has been found to be slightly in error (Beck & Giles 2005; Hathaway & Rightmire 2010), corrections to both the instrument roll angle and the solar tilt angle B₀ are needed, otherwise a strong annual variation appears in the measured travel-time shifts and modulate the resulting spectrum. Accordingly, an improvement to the mapping procedure used in Liang et al. (2018) is made by including corrections to the values of the two keywords CRLT_OBS (B₀ angle) and CROTA2 (instrument roll angle). Following Larson & Schou (2015, Eqs. (3) and (4)), the first order corrections are implemented as

$$\begin{array}{cc}& \mathtt{C}RLT\_OBS{}^{\prime }=\mathtt{C}RLT\_OBS+\delta Isin{\displaystyle \left(2\pi \frac{t-t_{\mathrm{ref}}}{1\mathrm{yr}}\right)},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& \mathtt{C}ROTA2{}^{\prime }=\mathtt{C}ROTA2+\delta P+\delta Icos{\displaystyle \left(2\pi \frac{t-t_{\mathrm{ref}}}{1\mathrm{yr}}\right)},\hfill \end{array}$$(A.2)

where the primed keywords denote the updated values, δI is the error in the inclination angle, δP is an instrument-specific correction for CCD misalignment, t and t_ref are the observation time and the time when CRLT_OBS is close to zero, both expressed in years (365.25 days). The t_ref is determined by fitting a sinusoid with a fixed period of 365.25 days to a time series of keyword CRLT_OBS, and is 7 June 1996 01:19:34_TAI for MDI data and 7 June 2010 14:17:20_TAI for HMI data. The adopted value of δP for MDI data is 0.2° (Liang et al. 2017, 2018), while that for HMI data is zero since the HMI team had calibrated and updated the value of CROTA2 (Couvidat et al. 2016; Hoeksema et al. 2018). As for δI, a value of −0.08° from Hathaway & Rightmire (2010) is adopted. The effect of this correction is quantified by fitting a periodic function to the south-north travel-time shifts (see Appendix A.3).

A.2. Weighting function

Here we use the notation introduced in Sect. 2.1 and denote the measured south-north travel-time shifts by δτ⁽⁰⁾(ψ, t, Δ), where ψ is the Carrington longitude, t the observation time, and Δ the travel distance. We have dropped the indices i and j in the subscripts of ψ and t for simplicity. To combine the travel-time shifts for different travel distances, we follow the procedure in Jackiewicz et al. (2008, Sect. 5.1). Given two random variables X and Y, the covariance is defined by Cov[X, Y]=⟨XY⟩−⟨X⟩⟨Y⟩, where the angle brackets denote the expectation value (ensemble average). The noise covariance matrix of the travel-time shifts between different distances is given by Λ(Δ, Δ′) = Cov[δτ⁽⁰⁾(ψ, t, Δ),δτ⁽⁰⁾(ψ, t, Δ′)]. The weighting for each distance is then given by

$$\begin{array}{c}\hfill w(\mathrm{\Delta })=\frac{{{\displaystyle \sum }}_{\mathrm{\Delta }\prime }{\mathrm{\Lambda }}^{-1}(\mathrm{\Delta },\mathrm{\Delta }\prime )}{{{\displaystyle \sum }}_{\mathrm{\Delta },\mathrm{\Delta }\prime }{\mathrm{\Lambda }}^{-1}(\mathrm{\Delta },{\mathrm{\Delta }}^{\prime })},\end{array}$$(A.3)

and the weighted average of the measured travel-time shifts is

$$\begin{array}{c}\hfill \delta {\tau }^{(1)}(\psi ,t)={\displaystyle \sum _{\mathrm{\Delta }}}w(\mathrm{\Delta })\delta {\tau }^{(0)}(\psi ,t,\mathrm{\Delta }).\end{array}$$(A.4)

The square root of Λ(Δ, Δ) is the standard deviation of the travel-time shifts, denoted by σ(Δ), which gives an estimate of the noise level at each distance. Figure A.1 shows the σ(Δ) 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.

Fig. A.1. Top: standard deviation of the measured travel-time shifts (square root of the diagonal elements of Λ) as a function of distance. 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⊙.

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 ψ_ij, and the central meridian as seen by the observer is denoted by ψ_0, j. Because the discussion here involves the center-to-limb effects, it is convenient to define the separation from the central meridian by ψ′_ij ≡ ψ_ij − ψ_0, j. The weighted average of the travel-time shifts resulting from Appendix A.2 is expressed in this coordinate system and denoted by δτ⁽¹⁾(ψ′_ij, t_j).

At each longitude ψ′_ij, a function representing the temporal variation of the background

$$\begin{array}{c}\hfill \begin{array}{cc}S({\psi }_{\mathit{ij}}^{\prime },t_j)\hfill & =a_0({\psi }_{\mathit{ij}}^{\prime })+a_1({\psi }_{\mathit{ij}}^{\prime })(t_j-t_{\mathrm{ref}})\hfill \\ \hfill & +[a_2({\psi }_{\mathit{ij}}^{\prime })+a_3({\psi }_{\mathit{ij}}^{\prime })(t_j-t_{\mathrm{ref}})]sin\left(2\pi \frac{t_j-t_{\mathrm{ref}}}{1\mathrm{yr}}\right)\hfill \\ \hfill & +a_4({\psi }_{\mathit{ij}}^{\prime })cos\left(2\pi \frac{t_j-t_{\mathrm{ref}}}{1\mathrm{yr}}\right),\hfill \end{array}\end{array}$$(A.5)

is fitted to and subtracted from the δτ⁽¹⁾(ψ′_ij, t_j) for MDI and HMI data separately. Here t_ref is taken from Appendix A.1. The a₀(ψ′_ij) represents the time-independent background, and a₁(ψ′_ij) represents the linear component, if any, of the slowly varying background. The a₂(ψ′_ij) accounts for an annual variation due to systematic errors such as the center-to-limb effects that vary with the B₀ angle (CRLT_OBS); in addition, the a₃(ψ′_ij) 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₄(ψ′_ij) 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 δτ⁽¹⁾(ψ′_ij, t_j), the fitted background S(ψ′_ij, t_j), and the residuals

Fig. A.2. Left: travel-time measurements δτ(1)(ψ′ij, tj), Gaussian smoothed with FWHM = 90 days in time for better visualization, where ψ′ij ≡ ψij − ψ0, j is the separation from the central meridian. Middle: fitted background S(ψ′ij, tj). Right: difference between the two panels on the left. The horizontal dashed lines in all panels indicate the division between the results from MDI and HMI data; the MDI data used in this work span from May 1996 to April 2010 while HMI data span from May 2010 to April 2017.

$$\begin{array}{c}\hfill \delta {\tau }^{(2)}(\psi {\prime }_{\mathit{ij}},t_j)=\delta {\tau }^{(1)}(\psi {\prime }_{\mathit{ij}},t_j)-S(\psi {\prime }_{\mathit{ij}},t_j).\end{array}$$(A.6)

The annual variation is clearly seen in the δτ⁽¹⁾(ψ′_ij, t_j), suggesting the necessity of dealing with this systematic effect. Also, the magnitude of δτ⁽¹⁾(ψ′_ij, t_j) from HMI data at larger |ψ′_ij| is systematically greater than the rest. The fitted S(ψ′_ij, t_j) accounts for a large amount of the periodic variation of the background. After removing the fitted background, the resulting δτ⁽²⁾(ψ′_ij, 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(ψ′_ij, 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₀(ψ′_ij) and a₂(ψ′_ij) (i.e., the center-to-limb effects) dominate the background variation while a₄(ψ′_ij) is nearly zero (0.008 ± 0.021 s for MDI and −0.015 ± 0.038 s for HMI within the longitude range |ψ′_ij|≤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₄(ψ′_ij) 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₄(ψ′_ij) becomes negative in this case (−0.035 ± 0.019 s for MDI and −0.058 ± 0.038 s for HMI).

Fig. A.3. Panel a: fitted S(ψ′ij, tj) at selected longitudes, where ψ′ij ≡ ψij − ψ0, j is the separation from the central meridian. The vertical dashed line indicates the division between the results from MDI and HMI data. Panels b–e: fitted values of parameters as a function of longitude for the MDI and HMI data sets, respectively. The error bars are estimated from the covariance matrix of the fitting parameters.

The background variation in the measurements from HMI is larger than that from MDI since the magnitude of center-to-limb effects for HMI is greater (Liang et al. 2017). The value of a₂(ψ′_ij) 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-to-limb 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 travel-time measurements from MDI was reported by Giles (2000, Sect. 4.7.2). He attributed the asymmetry to the nonuniform 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, 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

Figure A.4 shows the power spectrum computed in the frame of the observer (Earth’s frame).

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.

All Tables

All Figures

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.

In the text

Fig. 2. 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.

In the text

Fig. 3. Section of the window function W(ϕij, tj) in the frame that rotates at Ωeq/2π = 453.1 nHz. At time tj the longitudes of the travel times are denoted by ϕij (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.

In the text

	Fig. 4. Window function CT(t) (thick solid line) that selects an observation period of T = 2556Δt ≈ 7 yr, tapered with one-year raised cosines at both ends.
In the text

Fig. 5. Panel a: mean power spectrum of the window functions $W_{\mathrm{T}}^{(k)}$. 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.

In the text

Fig. 6. Power spectra of south-north travel-time shifts measured in the frame rotating at Ωeq/2π = 453.1 nHz. Panels a–c: $P_{\mathrm{T}}^{(k)}(m,\omega )$ from three time periods. Panel d: mean spectrum $\overline{P}(m,\omega )$ from an average of the three $P_{\mathrm{T}}^{(k)}(m,\omega )$ (seven years each). The blue lines highlight the dispersion relation of the classical Rossby waves described in Eq. (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.

In the text

Fig. 7. Enlargement of Fig. 6d around the frequency range of interest (see Fig. 8 for the line profiles of individual modes in the range 3 ≤ m ≤ 14). The red circles show the mode frequencies estimated from Lorentzian fits (see Sect. 2.3). The errors in the mode frequencies are given roughly by the size of the red circles (see errors in Table 1). For comparison, the mode frequencies measured by Löptien et al. (2018) are also indicated by green circles.

In the text

Fig. 8. Power spectra of south-north travel-time shifts $\overline{P}(m,\omega )$ for modes in the range 3 ≤ m ≤ 14 (black curves, with the frequency resolution 1/T = 4.5 nHz). The travel times are measured in the frame rotating at equatorial rotation rate Ωeq/2π = 453.1 nHz. The red lines are the fits Fm(ω) given by Eq. (8). The green vertical lines indicate the mode frequencies from Löptien et al. (2018). The orange arrows mark the excess low-frequency power at low m that might be caused by active regions or large-scale convection. The purple arrows mark the leaks from the m − 1 r-mode power and the low-frequency power at m − 1. Each power spectrum is accompanied by a plot of the ratio $\overline{P}(m,\omega )/F_{\mathrm{m}}(\omega )$ in gray in the lower panel, which is expected to have a mean of unity (dashed line) and a constant variance if the fit is not biased.

In the text

Fig. 9. Power spectrum of south-north travel-time shifts $\overline{P}(m,\omega )$ for m = 2 (black solid line). The blue vertical line indicates the frequency of the classical m = 2 sectoral r mode. The red solid line is the background B estimated by a fit to the power in the frequency range between −350 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.

In the text

	Fig. 10. Schematic damped oscillation of a sectoral r mode with m = 5, seen at a fixed longitude in the Sun’s corotating frame. The value of the e-folding lifetime (vertical line) is as observed (see 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.
In the text

Fig. A.1. Top: standard deviation of the measured travel-time shifts (square root of the diagonal elements of Λ) as a function of distance. 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⊙.

In the text

Fig. A.2. Left: travel-time measurements δτ(1)(ψ′ij, tj), Gaussian smoothed with FWHM = 90 days in time for better visualization, where ψ′ij ≡ ψij − ψ0, j is the separation from the central meridian. Middle: fitted background S(ψ′ij, tj). Right: difference between the two panels on the left. The horizontal dashed lines in all panels indicate the division between the results from MDI and HMI data; the MDI data used in this work span from May 1996 to April 2010 while HMI data span from May 2010 to April 2017.

In the text

Fig. A.3. Panel a: fitted S(ψ′ij, tj) at selected longitudes, where ψ′ij ≡ ψij − ψ0, j is the separation from the central meridian. The vertical dashed line indicates the division between the results from MDI and HMI data. Panels b–e: fitted values of parameters as a function of longitude for the MDI and HMI data sets, respectively. The error bars are estimated from the covariance matrix of the fitting parameters.

In the text

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.

In the text