Time-distance helioseismology of solar Rossby waves

TL;DR

This study uses 21 years of helioseismic data to confirm the existence of solar Rossby waves in deeper layers, analyzing their frequencies, lifetimes, and velocities, providing independent validation and new insights into their properties.

Contribution

It applies time-distance helioseismology to long-term data to confirm and extend previous surface observations of solar Rossby waves into deeper solar layers.

Findings

Rossby wave frequencies match previous surface studies for m ≤ 9

Modes with m ≥ 10 are shifted toward less negative frequencies

Longest-lived mode (m=3) persists over more than a year

Abstract

Context. Solar Rossby waves (r modes) have recently been discovered in the near-surface horizontal flow field using 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 combined SOHO/MDI and SDO/HMI data sets, which cover 21 years of observations from May 1996 to April 2017. The goal of this study is to provide an independent confirmation over two solar cycles and in deeper layers of the Sun. Methods. We have measured south-north helioseismic travel times along the equator, which are sensitive to subsurface north-south flows. To reduce noise, the travel times were 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 was split into three seven-year subsets and transformed to obtain power spectra in a corotating frame. Results. The power spectra all show peaks near the frequencies of the classical sectoral 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 whereas 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. For each mode, the rms velocity at the equator is in the range of 1--3 m s$^{-1}$ , with the largest values for $m\sim10$. No evidence for the $m=2$ sectoral mode is found in the power spectrum, implying that the rms velocity of this mode is below $\sim$0.5 m s$^{-1}$.

Tables

Table 1. Table 1: Measured characteristics of solar sectoral r modes.

	This work									Löptien et al. (2018)
	1996–2017 (MDI & HMI)									2010–2016 (HMI)
						period in	e-folding
$m$	${\nu }_{\mathrm{start}}$	${\omega }_m/2\pi$	${\Gamma }_m/2\pi$	$A$	$B$	${\Omega }_{\mathrm{eq}}$ frame	lifetime	$\delta {\tau }_{\mathrm{rms}}$	$v_{\mathrm{rms}}$	${\omega }_m/2\pi$
	(nHz)	(nHz)	(nHz)			(days)	(days)	(s)	(m s-1)	(nHz)
3	$-330$	$-253\pm 2$	$7_{-3}^{+4}$	${165}_{-77}^{+173}$	$18\pm 1$	$-46$	$505$	$0.07$	$1.0$	$-259$
4	$-350$	$-198\pm 5$	${34}_{-13}^{+15}$	${51}_{-14}^{+27}$	${15}_{-2}^{+1}$	$-58$	$109$	$0.08$	$1.3$	$-{194}_{-4}^{+5}$
5	$-330$	$-156\pm 2$	$8_{-4}^{+5}$	${123}_{-56}^{+111}$	$15\pm 1$	$-74$	$445$	$0.06$	$1.0$	$-157\pm 4$
6	$-310$	$-{135}_{-5}^{+4}$	${20}_{-11}^{+12}$	${50}_{-17}^{+48}$	${19}_{-2}^{+1}$	$-86$	$181$	$0.06$	$1.1$	$-129\pm 8$
7	$-280$	$-110\pm 4$	${40}_{-10}^{+12}$	${89}_{-22}^{+31}$	$11\pm 1$	$-105$	$92$	$0.11$	$2.2$	$-112\pm 4$
8	$-240$	$-91\pm 3$	${19}_{-6}^{+7}$	${131}_{-44}^{+72}$	${18}_{-2}^{+1}$	$-128$	$198$	$0.09$	$1.9$	$-90\pm 3$
9	$-236$	$-82\pm 5$	${33}_{-12}^{+13}$	${47}_{-13}^{+24}$	$12\pm 1$	$-142$	$112$	$0.08$	$1.7$	$-86\pm 6$
10	$-225$	$-{60}_{-6}^{+5}$	${54}_{-16}^{+18}$	${44}_{-10}^{+16}$	$9\pm 1$	$-193$	$68$	$0.09$	$2.2$	$-75\pm 5$
11	$-225$	$-48\pm 7$	${84}_{-24}^{+27}$	${24}_{-4}^{+7}$	$5\pm 1$	$-239$	$44$	$0.08$	$2.2$	$-75\pm 7$
12	$-209$	$-36\pm 8$	${75}_{-28}^{+31}$	${13}_{-2}^{+5}$	$5\pm 1$	$-323$	$49$	$0.06$	$1.7$	$-59\pm 6$
13	$-190$	$-26\pm 7$	${82}_{-25}^{+29}$	${13}_{-2}^{+4}$	$3\pm 1$	$-449$	$45$	$0.06$	$2.0$	$-40\pm 10$
14	$-206$	$-35\pm 5$	${27}_{-12}^{+13}$	${11}_{-4}^{+8}$	$4\pm 0.4$	$-326$	$135$	$0.03$	$1.2$	$-{56}_{-7}^{+6}$
15	$-197$	$-{22}_{-3}^{+2}$	${11}_{-5}^{+7}$	${19}_{-8}^{+19}$	$4\pm 0.3$	$-519$	$345$	$0.03$	$1.1$	$-{47}_{-6}^{+7}$

Equations

$$\omega \approx -\frac{2m\Omega }{\ell (\ell +1)}\text{for}m>0,$$

$${\psi }_{ij}={\psi }_{00}+i\Delta \psi -({\Omega }_{\mathrm{cr}}-{\Omega }_{\oplus })(t_j-t_0),$$

$$t_0=\text{1 May 1996 12:00:00\_TAI}.$$

$${\varphi }_{ij}={\psi }_{ij}+({\Omega }_{\mathrm{cr}}-{\Omega }_{\mathrm{eq}})(t_j-t_0),$$

$$t^{(1)}$$

$$t^{(2)}$$

$$t^{(3)}$$

$$W_T^{(k)}({\varphi }_{ij},t_j)=C_T(t_j-t^{(k)})W({\varphi }_{ij},t_j),k=\{1,2,3\},$$

$$P_T^{(k)}(m,\omega )=\frac{1}{N}{\left|\sum _{i,j}{W}_T^{(k)}({\varphi }_{ij},t_j)\delta \tau ({\varphi }_{ij},t_j)e^{-\mathrm{i}m{\varphi }_{ij}+\mathrm{i}\omega t_j}\right|}^2,$$

$$\overline{P}(m,\omega )=\frac{1}{3}\left(P_T^{(1)}(m,\omega )+P_T^{(2)}(m,\omega )+P_T^{(3)}(m,\omega )\right).$$

$$F_m(\omega )=\frac{A}{1+[(\omega -{\omega }_m)/({\Gamma }_m/2){]}^2}+B,$$

$$\sum _{\omega }\left(\ln F_m(\omega )+\frac{\overline{P}(m,\omega )}{F_m(\omega )}\right)$$

$$CRLT\_OB{\mathtt{S}}^{\prime }$$

$$\mathtt{CROTA}{2}^{\prime }$$

$$w(\Delta )=\frac{{\sum }_{{\Delta }^{\prime }}{\Lambda }^{-1}(\Delta ,{\Delta }^{\prime })}{{\sum }_{\Delta ,{\Delta }^{\prime }}{\Lambda }^{-1}(\Delta ,{\Delta }^{\prime })},$$

$$\delta {\tau }^{(1)}(\psi ,t)=\sum _{\Delta }w(\Delta )\delta {\tau }^{(0)}(\psi ,t,\Delta ).$$

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

$$\delta {\tau }^{(2)}({\psi }_{ij}^{\prime },t_j)=\delta {\tau }^{(1)}({\psi }_{ij}^{\prime },t_j)-S({\psi }_{ij}^{\prime },t_j).$$

Full text

*Aims. * Here we apply time-distance helioseismology to the combined SOHO/MDI and SDO/HMI data sets, which cover 21 years of observations from May 1996 to April 2017. The goal of this study is to provide an independent confirmation over two solar cycles and in deeper layers of the Sun.

*Methods. * We have measured south-north helioseismic travel times along the equator, which are sensitive to subsurface north-south flows. To reduce noise, the travel times were averaged over travel distances from 6$\degr$ to 30$\degr$; 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 was split into three seven-year subsets and transformed to obtain power spectra in a corotating frame.

Results. * The power spectra all show peaks near the frequencies of the classical sectoral 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 whereas 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. For each mode, the rms velocity at the equator is in the range of 1 – 3 m s-1*, with the largest values for $m\sim 10$. No evidence for the $m=2$ sectoral mode is found in the power spectrum, implying that the rms velocity of this mode is below $\sim$0.5 m s*-1*.

*Conclusions. * This work confirms the existence of equatorial global Rossby waves in the solar interior over the past two solar cycles and shows that time-distance helioseismology is a promising technique to study them deep in the convection zone.

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 $\Omega$, the frequencies of r modes are given by

[TABLE]

where $\ell$ 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 $(\ell,m)$ drifts at a phase velocity of $\omega/m=-2\Omega/[\ell(\ell+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=\ell$ 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=\ell=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).

2 Data analysis and results

2.1 Time-distance analysis

We use medium-$\ell$ 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$\degr$ 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$\degr$-wide concentric arcs, separated by an angular distance $\Delta$, 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 $\pm 15\degr$ latitude at intervals of 0.6$\degr$ in longitude and latitude, and for the distance range $\Delta=6\degr$–$30\degr$ in steps of 0.6$\degr$. 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 $\text{FWHM}=12\degr$ in longitude and latitude, and subsampled along the equator at intervals of 10$\degr$ in longitude.

The south-north travel-time shift in this work is defined as $\delta\tau=\tau_{\rm s}-\tau_{\rm n}$, where $\tau_{\rm n}$ and $\tau_{\rm 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$\degr$, 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 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$\degr$, 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).111 Because of the inclination of the solar rotation axis to the ecliptic, the solar equator does not pass through disk center when the solar tilt angle ($B_{0}$ angle) is not zero. Therefore the measurements along the equator are affected by center-to-limb effects, especially in spring and autumn.

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 $\Omega_{\rm cr}/2\pi=456.03$ nHz in a sidereal frame. As viewed from the Earth, the Carrington coordinate system rotates at $\Omega_{\rm cr}-\Omega_{\oplus}$, where $\Omega_{\oplus}/2\pi=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\Delta t$, where $\Delta t=24$ h and $j\in\{0,1,\dots,7669\}$, covering 21 years of data. At each time $t_{j}$ the Carrington longitudes of the travel-time measurements are given by

[TABLE]

where $i\in\{-5,-4,\dots,5\}$ (11 measurements along the equator each day) and $\Delta\psi=10^{\circ}$ is the spatial sampling rate. The longitude $\psi_{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:

[TABLE]

To compare with the results of Löptien et al. (2018), we transform to a frame that rotates at $\Omega_{\rm eq}/2\pi=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

[TABLE]

where $\phi_{ij}$ is the longitude measured in the new frame. Figure 3 is a schematic plot of the corresponding window function $W(\phi_{ij},t_{j})$, which is equal to 1 for $i\in\{-4,\dots,4\}$, tapered to $1/2$ at the boundaries $i=\pm 5$ where noise is slightly higher, and zero elsewhere. The window function is also zero for missing data.

We denote by $\delta\tau(\phi_{ij},t_{j})$ the travel times interpolated at the longitudes $\phi_{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[-\mathrm{i}m(\Omega_{\rm cr}-\Omega_{\rm eq})(t_{j}-t_{0})]$, where $m$ is the azimuthal wavenumber.

2.2 Power spectra

The full data set $\delta\tau(\phi_{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

[TABLE]

The window functions for each 7-yr period is

[TABLE]

where $C_{T}(t_{j}-t^{(k)})$ is a taper time window that selects an observation period of $T=N\Delta t\approx 7$ yr where $N=2556$, starting from time $t^{(k)}$ (see Fig. 4).

The power spectrum of each 7-yr period is

[TABLE]

where the azimuthal order $m$ is in the range $|m|\leq 18$ (spatial Nyquist frequency). For each $P_{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_{T}^{(k)}(\phi_{ij},t_{j})$ by $\widehat{W}_{T}^{(k)}(m,\omega)$. Because the resulting $P_{T}^{(k)}(m,\omega)$ involves the convolution of $\widehat{W}_{T}^{(k)}(m,\omega)$, it is necessary to investigate the power spectrum of the window function, $|\widehat{W}_{T}^{(k)}(m,\omega)|^{2}$, from which the spectral leakage arises. We note that $\widehat{W}_{T}^{(k)}(-m,-\omega)=\widehat{W}_{T}^{(k)*}(m,\omega)$ since the window function is real.

The mean power spectrum of the window functions, $\sum_{k=1}^{3}|\widehat{W}_{T}^{(k)}(m,\omega)|^{2}/3$, consists of well-defined peaks, as can be seen in Fig. 5a (plot restricted to $m\geq 0$). The maximum values of the peaks at $m=1$ and $m=2$ are 76% and 30% of the peak at the origin $(m,\omega)=(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$\degr$ longitudes. In Fig. 5a, adjacent peaks are shifted in frequency by $(\Omega_{\rm eq}-\Omega_{\oplus})/2\pi=421.41$ nHz due to the choice of coordinate frame. Since the power spectrum $P_{T}^{(k)}$ involves a convolution in Fourier space of the solar data with $\widehat{W}_{T}^{(k)}$, the power spectrum of the data at $(m,\omega)$ leaks to multiple sidelobes at \big{(}m+\delta m,\omega-\delta m\,(\Omega_{\rm eq}-\Omega_{\oplus})\big{)}, where $|\delta m|=1,2,\dots$.

Figure 5b shows a cut at $m=0$ through the power spectrum of the window function. The width of the main lobe along the $\omega$-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_{T}^{(k)}(m,\omega)$ for the three different periods, as well as the mean spectrum

[TABLE]

The power distribution in all three $P_{T}^{(k)}(m,\omega)$ peaks around the eigenfrequencies of r modes as described in Eq. 1 for the case of $\ell=m\geq 3$. However, there is no clear sign of $\ell=m\leq 2$. It is interesting to note that the mode amplitudes vary among the three $P_{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 $\omega=-2\Omega_{\mathrm{eq}}/(m+1)$ becomes more prominent. The other two ridges of power seen in Fig. 6, which are separated from the central ridge by (\delta m,\delta\omega)=\big{(}\pm 1,\mp(\Omega_{\rm eq}-\Omega_{\oplus})\big{)} with a reduced amplitude ($60\sim 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 $\pm 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 ${\rm FWHM=12^{\circ}}$ 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\geq 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\times 421.41$ nHz for each $m$ value. Figure 14 displays a set of horizontal segments of size set by the resolution in $m$ ($\approx 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,

[TABLE]

is fit to the $\overline{P}(m,\omega)$ at individual $m$. Here $A$ is the maximum height of the Lorentzian function, $\omega_{m}$ the mode frequency, $\Gamma_{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 [$\nu_{\rm start}$, $\nu_{\rm start}$+300 nHz], with $\nu_{\rm start}$ given in Table 2.

The best fit is obtained by minimizing the sum

[TABLE]

with respect to the model parameters $A,B,\omega_{m},\Gamma_{m}$ (Duvall & Harvey 1986). Although Eq. 9 was derived assuming that the quantity $2\overline{P}(m,\omega)/F_{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 2 lists the resulting best fit parameters for $3\leq m\leq 15$.

Figure 8 shows the fitted $F_{m}(\omega)$ 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\leq m\leq 9$. However, excess power on the right of the main peak is present for $m\geq 8$, the extent of which for $m\geq 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*-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_{\odot}$. 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*-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 $\sim$1.6 for each mode.

Next, we obtain the rms travel-time shifts $\delta\tau_{\rm rms}$ from the fitted Lorentzian profile $F_{m}(\omega)-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 $\delta\tau_{\rm rms}$.

Last, we use the above conversion constants to convert from $\delta\tau_{\rm rms}$ to the surface rms flow speed $v_{\rm rms}$. The conversion constants adopted are from the first model in which the flow is independent of depth, and thus the estimated $v_{\rm rms}$ are conservative. The results for modes with $3\leq m\leq 15$ are listed in Table 2. For these modes, the surface velocity is on the order of 1 – 3 m s*-1*, with larger values for $7\leq m\leq 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).

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 $\sim$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\leq m=\ell\leq 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\geq 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 2 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,\Gamma_{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 $\pm 20\degr$ 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_{\rm rms}$ in Table 2.

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\leq 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*-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_{\rm rms}<0.5$ m s*-1* would be difficult to identify in the power spectrum.

The e-folding lifetimes of r modes, $2/\Gamma_{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.

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.

Acknowledgements.

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.

Appendix A Supplementary material

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_{0}$ 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_{0}$ angle) and CROTA2 (instrument roll angle). Following Larson & Schou (2015 , Eqs.3 and 4), the first order corrections are implemented as

[TABLE]

where the primed keywords denote the updated values, $\delta I$ is the error in the inclination angle, $\delta P$ is an instrument-specific correction for CCD misalignment, $t$ and $t_{\mathrm{ref}}$ are the observation time and the time when CRLT_OBS is close to zero, both expressed in years (365.25 days). The $t_{\mathrm{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 $\delta P$ for MDI data is 0.2$\degr$ (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 $\delta I$, a value of $-0.08\degr$ 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 $\delta\tau^{(0)}(\psi,t,\Delta)$, where $\psi$ is the Carrington longitude, $t$ the observation time, and $\Delta$ the travel distance. We have dropped the indices $i$ and $j$ in the subscripts of $\psi$ 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 $\mathrm{Cov}[X,Y]=\langle XY\rangle-\langle X\rangle\langle Y\rangle$, 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 $\Lambda(\Delta,\Delta^{\prime})=\mathrm{Cov}[\delta\tau^{(0)}(\psi,t,\Delta),\delta\tau^{(0)}(\psi,t,\Delta^{\prime})]$. The weighting for each distance is then given by

[TABLE]

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

[TABLE]

The square root of $\Lambda(\Delta,\Delta)$ is the standard deviation of the travel-time shifts, denoted by $\sigma(\Delta)$, which gives an estimate of the noise level at each distance. Figure 11 shows the $\sigma(\Delta)$ and $w(\Delta)$ as a function of distance. The trend of $\sigma$ 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 $\Delta<10\degr$, the $\sigma$ increases rapidly because of the low spatial resolution of medium-$\ell$ 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 $\Delta>10\degr$, the $\sigma$ increases with increasing distance because of geometrical spreading and damping. Also, the noise correlation between different distances (off-diagonal elements of $\Lambda$) becomes greater when both $\Delta$ and $\Delta^{\prime}$ 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 $\psi_{ij}$, and the central meridian as seen by the observer is denoted by $\psi_{0,j}$. Because the discussion here involves the center-to-limb effects, it is convenient to define the separation from the central meridian by $\psi^{\prime}_{ij}\equiv\psi_{ij}-\psi_{0,j}$. The weighted average of the travel-time shifts resulting from Appendix A.2 is expressed in this coordinate system and denoted by $\delta\tau^{(1)}(\psi^{\prime}_{ij},t_{j})$.

At each longitude $\psi^{\prime}_{ij}$, a function representing the temporal variation of the background

[TABLE]

is fitted to and subtracted from the $\delta\tau^{(1)}(\psi^{\prime}_{ij},t_{j})$ for MDI and HMI data separately. Here $t_{\rm ref}$ is taken from Appendix A.1. The $a_{0}(\psi^{\prime}_{ij})$ represents the time-independent background, and $a_{1}(\psi^{\prime}_{ij})$ represents the linear component, if any, of the slowly varying background. The $a_{2}(\psi^{\prime}_{ij})$ 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}(\psi^{\prime}_{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_{4}(\psi^{\prime}_{ij})$ accounts for another annual variation caused by the error in the inclination of the rotation axis, if not completely removed in Sect. A.1; the resulting error in the roll angle (i.e., Eq. 11) may introduce a leakage of the solar rotation signal into the south-north travel-time shifts (Beck & Giles 2005).

Figure 12 shows the $\delta\tau^{(1)}(\psi^{\prime}_{ij},t_{j})$, the fitted background $S(\psi^{\prime}_{ij},t_{j})$, and the residuals

[TABLE]

The annual variation is clearly seen in the $\delta\tau^{(1)}(\psi^{\prime}_{ij},t_{j})$, suggesting the necessity of dealing with this systematic effect. Also, the magnitude of $\delta\tau^{(1)}(\psi^{\prime}_{ij},t_{j})$ from HMI data at larger $|\psi^{\prime}_{ij}|$ is systematically greater than the rest. The fitted $S(\psi^{\prime}_{ij},t_{j})$ accounts for a large amount of the periodic variation of the background. After removing the fitted background, the resulting $\delta\tau^{(2)}(\psi^{\prime}_{ij},t_{j})$ from the two data sets are more consistent with each other, and are ready for the Fourier analysis.

Figure 13 shows the $S(\psi^{\prime}_{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_{0}(\psi^{\prime}_{ij})$ and $a_{2}(\psi^{\prime}_{ij})$ (i.e., the center-to-limb effects) dominate the background variation while $a_{4}(\psi^{\prime}_{ij})$ is nearly zero ($0.008\pm 0.021$ s for MDI and $-0.015\pm 0.038$ s for HMI within the longitude range $|\psi^{\prime}_{ij}|\leq 30\degr$) since we have applied a correction $\delta I=-0.08\degr$ to the inclination angle of solar rotation axis as discussed in Appendix A.1. Without this correction, the value of $a_{4}(\psi^{\prime}_{ij})$ is around 0.1–0.2 s. We also tried $\delta I=-0.1\degr$ used by Larson & Schou (2015) who rounded the value $-0.095\degr\pm 0.002\degr$ determined by Beck & Giles (2005) up to $-0.1\degr$; however, the $a_{4}(\psi^{\prime}_{ij})$ becomes negative in this case ($-0.035\pm 0.019$ s for MDI and $-0.058\pm 0.038$ s for HMI).

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_{2}(\psi^{\prime}_{ij})$ is expected to be negative because when the solar north pole is tilted toward the observer ($\mathtt{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 & Toomre 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 14 shows the power spectrum computed in the frame of the observer (Earth’s frame).