Issue 
A&A
Volume 634, February 2020



Article Number  A44  
Number of page(s)  16  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201937007  
Published online  04 February 2020 
Exploring the latitude and depth dependence of solar Rossby waves using ringdiagram analysis
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: proxauf@mps.mpg.de
^{2}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3}
Center for Space Science, NYUAD Institute, New York University Abu Dhabi, Abu Dhabi, UAE
^{4}
W. W. Hansen Experimental Physics Laboratory, Stanford University, Stanford, CA 94305, USA
Received:
28
October
2019
Accepted:
7
December
2019
Context. Globalscale equatorial Rossby waves have recently been unambiguously identified on the Sun. Like solar acoustic modes, Rossby waves are probes of the solar interior.
Aims. We study the latitude and depth dependence of the Rossby wave eigenfunctions.
Methods. By applying helioseismic ringdiagram analysis and granulation tracking to observations by HMI aboard SDO, we computed maps of the radial vorticity of flows in the upper solar convection zone (down to depths of more than 16 Mm). The horizontal sampling of the ringdiagram maps is approximately 90 Mm (∼7.5°) and the temporal sampling is roughly 27 hr. We used a Fourier transform in longitude to separate the different azimuthal orders m in the range 3 ≤ m ≤ 15. At each m we obtained the phase and amplitude of the Rossby waves as functions of depth using the helioseismic data. At each m we also measured the latitude dependence of the eigenfunctions by calculating the covariance between the equator and other latitudes.
Results. We conducted a study of the horizontal and radial dependences of the radial vorticity eigenfunctions. The horizontal eigenfunctions are complex. As observed previously, the real part peaks at the equator and switches sign near ±30°, thus the eigenfunctions show significant nonsectoral contributions. The imaginary part is smaller than the real part. The phase of the radial eigenfunctions varies by only ±5° over the top 15 Mm. The amplitude of the radial eigenfunctions decreases by about 10% from the surface down to 8 Mm (the region in which ringdiagram analysis is most reliable, as seen by comparing with the rotation rate measured by globalmode seismology).
Conclusions. The radial dependence of the radial vorticity eigenfunctions deduced from ringdiagram analysis is consistent with a power law down to 8 Mm and is unreliable at larger depths. However, the observations provide only weak constraints on the powerlaw exponents. For the real part, the latitude dependence of the eigenfunctions is consistent with previous work (using granulation tracking). The imaginary part is smaller than the real part but significantly nonzero.
Key words: Sun: helioseismology / Sun: oscillations / Sun: interior / waves
© B. Proxauf et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
Recently, Löptien et al. (2018, hereafter LGBS18) discovered globalscale Rossby waves in maps of flows on the surface of the Sun. These waves are waves of radial vorticity that may exist in any rotating fluid body. Even though Rossby waves were predicted to exist in stars more than 40 years ago (Papaloizou & Pringle 1978; Saio 1982), solar Rossby waves were difficult to detect because of their small amplitudes (∼1 m s^{−1}) and long periods of several months. Solar Rossby waves contain almost as much vorticity as largescale solar convection. The dispersion relation of solar Rossby waves is close to the standard relation for sectoral modes, ω = −2Ω/(m + 1), where Ω is the rotation rate of a rigidly rotating star and m is the azimuthal order (Saio 1982). Rossby waves have a retrograde phase speed and a prograde group speed. In LGBS18, the authors also measured the horizontal eigenfunctions, which peak at the equator.
The detection of solar Rossby waves was confirmed by Liang et al. (2019, hereafter LGBD19) with timedistance helioseismology (Duvall et al. 1993) using data covering more than 20 years, obtained from the Solar and Heliospheric Observatory (SOHO) and from the Solar Dynamics Observatory (SDO; Pesnell et al. 2012). Alshehhi et al. (2019), in an effort to speed up ringdiagram analysis (RDA; Hill 1988) via machine learning, also saw globalscale Rossby waves. Hanasoge & Mandal (2019) and Mandal & Hanasoge (2019) provide another recent Rossby wave confirmation using a different technique of helioseismology known as normalmode coupling (Woodard 1989; Hanasoge et al. 2017).
Knowledge about the latitude dependence of Rossby wave eigenfunctions is incomplete, as LGBS18 studied only their real parts. In a differentially rotating star, the horizontal eigenfunctions are not necessarily spherical harmonics (and may not even separate in latitude and depth). Also, little is known observationally about the depth dependence of the Rossby waves. It would be well worth distinguishing between the few existing theoretical models of the depth dependence (Provost et al. 1981; Smeyers et al. 1981; Saio 1982; Wolff & Blizard 1986).
In this paper, we explore the latitude dependence of the eigenfunctions, as well as the phase and amplitude of solar Rossby waves as functions of depth from the surface down to more than 16 Mm using helioseismology. We use observations from the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) on board SDO, processed with RDA. From these we attempt to measure the eigenfunctions of the Rossby waves in the solar interior. For comparison near the surface, we also use data from local correlation tracking of granulation (LCT; November & Simon 1988).
2. Data and methods
We used maps of the horizontal velocity, derived from two different techniques applied to SDO/HMI observations. The first dataset consists of LCT (granulation tracking) flow maps at the surface (Löptien et al. 2017) and covers almost six years from May 20, 2010 to March 30, 2016. The second dataset comprises RDA flow maps from the HMI ringdiagram pipeline (Bogart et al. 2011a,b; see also Bogart et al. 2015). For comparisons with LCT, we took a period as close to the LCT period as possible, i.e., May 19, 2010 to March 31, 2016, while for all other results we used a longer period of more than seven years from May 19, 2010 to December 29, 2017; this corresponds to 102 Carrington rotations (CRs), i.e., CR 2097  2198.
2.1. Overview of LCT data
The LCT flow maps are obtained from and processed as described in Löptien et al. (2017). They are created by applying the Fourier LCT code (FLCT; Welsch et al. 2004; Fisher & Welsch 2008) to track the solar granulation in pairs of consecutive HMI intensity images. The image pairs are separated by 30 min. Several known systematic effects such as the shrinkingSun effect (Lisle et al. 2004; Löptien et al. 2016) and effects related to the SDO orbit are present in the LCT maps. Therefore the maps are decomposed into Zernike polynomials, a basis of 2D orthogonal functions on the unit disk, and the time series of the coefficient amplitudes for the lowest few Zernike polynomials are filtered to remove frequencies of one day and one year (associated with the SDO orbit) as well as all harmonics up to the Nyquist frequency. The zero frequency is also removed. The filtered maps are then tracked at the sidereal Carrington rate and remapped onto an equispaced longitudelatitude grid with a step size of 0.4° in both directions.
2.2. Overview of ringdiagram data
The ringdiagram pipeline (Bogart et al. 2011a,b) takes HMI Dopplergrams as input and remaps them onto tiles spanning 182 × 182 Mm (i.e., 15° each in latitude λ and longitude φ at the equator). The tiles overlap each other by roughly 50% in each direction such that the tile borders fall onto the centers of adjacent tiles. Both the latitude and longitude sampling are half the tile size. The latitude grid is linear and includes the equator, while the longitude grid is also linear, but is latitudedependent. Each tile is tracked for 1728 min (28.8 hr) at the sidereal Carrington rate. The temporal grid spacing is, on average, 1/24 of the synodic Carrington rotation period of 27.2753 days.
In the pipeline, for each tile a 3D local power spectrum is computed from the tracked Dopplergrams. The velocity fit parameters U_{x, nℓ} (prograde) and U_{y, nℓ} (northward) are extracted via a ringfit algorithm (Haber et al. 2000) for different solar oscillation modes, which are indexed by their radial order n and angular degree ℓ. The flow velocities u_{x} and u_{y} are inferred for various target depths via a 1D optimally localized averages (OLA) inversion. The inversion results for the sixparameter fits of the 15° tiles sample a range of target depths from 0.97 R_{⊙} to 1 R_{⊙} (step size 0.001 R_{⊙}), corresponding to a nonlinear grid of measurement depths (median of the ringdiagram averaging kernels) from 0.976 R_{⊙} to 1 R_{⊙}. In this paper, the term depth always refers to measurement depth and not to target depth.
The inversion results are stored in the Joint Science Operations Center (JSOC) data series hmi.V_rdvflows_fd15_frame. However, up to inversion module rdvinv v.0.91, the inversion results depended on the input tile processing order due to an array initialization bug. This caused significantly lower velocity uncertainties for tiles near latitude 7.5° and Stonyhurst longitude 37.5°, even when averaged over seven years, but also slightly affected the velocities. At the same disk locations the bug caused a correlation of u_{x} with the B_{0} angle. Since rdvinv v.0.92 is officially only applied since March 2018, we reinverted the entire dataset ourselves for the analysis shown in this work.
Apart from the array initialization bug, we found several other issues with the default HMI ringdiagram pipeline that have not yet been solved. Among these are underregularization in the inversion for some individual tiles, leading to relatively narrow averaging kernels and anomalously high noise. Finally, the number of ring fits used for the inversion depends strongly on disk position. This may lead to systematic effects and additional noise.
The ringdiagram velocities u_{x} reported at a certain measurement depth r at the equator for an angular rotation rate Ω(r) are equal to Ω(r)R_{⊙} instead of the local velocity Ω(r)r. Since we are interested in the latter, we multiplied u_{x} by r/R_{⊙}. By analogy, we also applied this factor to u_{y} and to all other latitudes. Additionally, the inversion does not account for the quantity β_{nℓ}, defined, for example, in Eq. (3.357) of Aerts et al. (2010). The quantity β_{nℓ} is related to the effect of the Coriolis force on the mode frequency splitting. For uniform rotation in particular, at fixed m, β_{nℓ} completely describes the effect of the rotation on the mode frequency splitting. Both issues are described in more detail in Appendix A.
2.3. Postprocessing of ringdiagram data
The ringdiagram data are organized in CRs, which undergo several processing steps, including the removal of systematic effects, an interpolation in longitude, an interpolation in time, and the removal of limb data.
Several systematic effects are present in the ringdiagram velocities, such as centertolimb effects that depend on the disk position of the tile (Baldner & Schou 2012; Zhao et al. 2012). There are timeindependent effects and systematics with a oneyear period, which are probably related to the B_{0} angle. To remove the systematics, we fit the time series at each position on the disk (in Stonyhurst coordinates) with sinusoids
and subtract the fits from the flow velocities. We used all available CRs to determine the fit parameters.
Because of the specific tile coordinate selection used by the ringdiagram pipeline (Bogart et al. 2011a), which seeks to optimally cover the visible disk, tile centers at different latitudes have Stonyhurst longitudes that are offset by multiples of 2.5° from each other. To obtain a latitudeindependent longitude grid, we interpolated the flow maps in Stonyhurst longitude using splines (Appendix B).
We also interpolated the ringdiagram flows in time similarly with splines to fill missing time steps due to instrumental issues (only 12 out of 2448 time steps), which cause a too low observational duty cycle (≤70%). We interpolated the data in the Carrington reference frame so as to use always roughly the same physical locations on the Sun. This mixes different systematics, which are primarily dependent on disk position, but we should already have removed the dominant contributions at this stage. We interpolated every missing time step from roughly the same number of data points (all available time steps within the corresponding disk passage) using splines (Appendix B).
The output uncertainties from the ringdiagram pipeline increase strongly toward the limb, in particular beyond an angular greatcircle distance of roughly 65° to the crossing of the central meridian with the equator (λ = 0° ,φ = 0°). We thus only used ringdiagram data within 65° of (λ = 0° ,φ = 0°).
2.4. From velocity maps to power spectra of radial vorticity
From this stage onward ringdiagram and LCT data are processed similarly. The processing steps include a shift to the equatorial rotation rate ν_{eq} = Ω_{eq}/2π = 453.1 nHz, the subtraction of the longitude mean, the calculation of the radial vorticity, a spherical harmonic transform (SHT), and a Fourier transform of the SHT coefficient time series.
The flow maps are shifted from the tracking rate (sidereal Carrington rate) to the surface sidereal equatorial rotation rate of 453.1 nHz, an average of zonal flows inferred from globalmode analysis of SDO/HMI observations (Larson & Schou 2018). We shifted the LCT data in Fourier space via a timedependent phase factor, applying the same convention for the Fourier transform as LGBS18. The ringdiagram data are first apodized by a raised cosine in angular greatcircle distance to the point (λ = 0° ,φ = 0°) to suppress nearlimb data and are shifted via splineinterpolation (Appendix B).
We next subtracted the longitude mean from the data to remove any remaining largescale flows. Differential rotation and meridional circulation should have already been subtracted in the RDA or LCT postprocessing, but any possible longitudeindependent flows still in the data are removed in this step.
Subsequently, we calculated the radial vorticity (via secondorder central finite differences) as follows:
where r is the measurement depth. We decomposed the resulting maps into spherical harmonics and performed a temporal Fourier transform of the spherical harmonic coefficients. Last, we calculated the power and the phase (where the phase range is the halfopen interval (−180° ,180°]). The sign convention is such that waves with positive m and negative frequency ν have a retrograde phase speed.
If not stated otherwise, the terms power spectrum or Fourier transform used in this paper always refer to the power spectrum or Fourier transform of the radial vorticity. Similarly, we discuss eigenfunctions of radial vorticity. These eigenfunctions are not spherical harmonics, however (LGBS18). In particular, while the modes can be meaningfully indexed by m, the angular degree ℓ is not observable. Throughout the paper ℓ thus only refers to the projection of the Rossby wave modes onto the corresponding spherical harmonic and not to the Rossby wave eigenfunction itself. We also use the terms latitudinal and radial eigenfunctions, which assumes separability in the r and λ coordinates. This assumption is addressed in more detail in Sect. 4.
3. Results
3.1. Radial vorticity maps
Figure 1 shows example vorticity maps from LCT surface flows and from RDA flows near the surface, at intermediate, and at large depths (0.7, 9.9, and 16.5 Mm), averaged over the first rotation in the dataset (May 20, 2010 to June 16, 2010). The LCT data have a much better horizontal resolution than the ringdiagram data and thus pick up smallscale convective contributions. To be able to compare LCT with RDA, we thus smooth the LCT vorticity with 1D Gaussian filters of width σ = 6° both in latitude and longitude.
Fig. 1. Radial vorticity maps from LCT at the surface and from RDA at depths 0.7, 9.9, and 16.5 Mm. The radial vorticity is averaged over one rotation (from May 20, 2010 to June 16, 2010). The LCT map is smoothed in latitude and longitude with a Gaussian filter (σ = 6°) to filter out smallscale convection. 

Open with DEXTER 
We do not expect perfect agreement of the two methods because of their different sensitivities to horizontal scales and to different depths. Nonetheless, the LCT map shows similar features as the nearsurface (0.7 Mm) ringdiagram map. While large absolute radial vorticities are visible at high latitudes (beyond ±50°) in the LCT but not in the ringdiagram data, the vorticities near the equator agree. As a test, we interpolate the LCT data to the RDA grid using a 2D bicubic spline. The correlation coefficient between the interpolated LCT and the ringdiagram maps decreases with the latitude width of the strip of pixels considered and there is a steep decrease beyond ±45°. The correlation is 0.92 when including only equatorial pixels, 0.79 for pixels within ±45°, and 0.59 for all pixels, i.e., within ±52.5°. The noise increases strongly with depth (see lower panels of Fig. 1), but the main vorticity features are still visible.
3.2. Power spectra of radial vorticity
Figure 2 shows the Rossby wave power of the ℓ = m component for the ringdiagram data near the surface (0.7 Mm) versus frequency and azimuthal order m (LGBS18 detected only the sectoral component of the Rossby waves). We divide the power, at each m, by the frequency average over [−300, 100] nHz near the surface (0.7 Mm). The visible power ridge corresponds to the Rossby wave signal. The mode frequency increases with m roughly according to the textbook dispersion relation for sectoral waves, ω = −2Ω_{eq}/(m + 1), as seen earlier by LGBS18.
Fig. 2. Sectoral power spectrum (ℓ = m) of the radial vorticity for RDA data at depth 0.7 Mm. The solid red line indicates the Rossby wave frequencies from LGBD19 for m = 3 and from LGBS18 for m ≥ 4. Frequency intervals for the Rossby wave region and the background region are indicated by the red and black dashed lines, respectively. The ridge of power at positive frequencies is due to the first side lobe of the window function. For better visibility of the Rossby waves at low m, the power is normalized at each m by the frequency average over [−300, 100] nHz. The color scale is truncated at 50% of the maximum value (black). 

Open with DEXTER 
Besides the Rossby wave signal there are other ridges, that, at fixed Δm = m − m′, are shifted from the Rossby waves by roughly Δm(ν_{eq}−1/(1 yr)), where ν_{eq} − 1/(1 yr)∼421.4 nHz. The main contribution to these side lobes comes from a temporal window function, which is introduced by solar rotation and not by time gaps; the time coverage is very good (see Sect. 2). This leads to side lobes of wave power from modes at m′ to modes at m. We only show the side lobes for Δm = +1, but we typically observe the side lobes above the noise between Δm = −2 and Δm = +3. In LGBD19, the authors use 21 years of timedistance data from a combined sample of observations from the Michelson Doppler Imager (MDI) on board SOHO and from SDO/HMI and they discuss the window function in detail.
Figure 3 shows the power versus frequency for different azimuthal orders m. We divide the power, at each m, by the frequency average of the ℓ = m = 8 mode over [−300, 100] nHz near the surface (0.7 Mm). The power decreases from 0.7 to 9.9 Mm, then increases toward 16.5 Mm, but the depth dependence is modest (≤20%). We also see that the wave power decreases with m faster for RDA than for LCT owing to the different sensitivity kernels, as found by LGBS18. The ℓ = m = 3 signal has a multipeak structure and is thus difficult to measure. We do not observe Rossby waves for ℓ = m ≤ 2; the dashdotted blue lines for ℓ = m = 1 and ℓ = m = 2 indicate the expected mode frequencies from the textbook dispersion relation.
Fig. 3. Sectoral power spectra (ℓ = m) showing the Rossby wave power in the LCT data (black line) and the RDA data at depths 0.7 and 16.5 Mm (cyan and red lines). The power is normalized by the average of the m = 8 power in the range [−300, 100] nHz. The dashdotted vertical lines in the m = 1 and m = 2 panels indicate the frequencies ω = −2Ω_{eq}/(m + 1). The frequency axes of the m = 1 and m = 2 panels are shifted with respect to the other panels. 

Open with DEXTER 
The wave power side lobes due to the window function explain why the ℓ = m = 6 side lobe in Fig. 2 even exceeds the main signal: the adjacent ℓ = m = 7 mode has a higher relative power (see Fig. 3). Systematic effects that are fixed in the Stonyhurst reference frame can be easily misinterpreted as an ℓ = m = 1 Rossby wave signal (see the LCT curve in Fig. 3), as their frequency (the rotation rate) is equal to the ℓ = m = 1 Rossby wave frequency.
We assume that there is background power contributing to the observed power at the Rossby peak, but measuring its contribution directly at the peak is impossible. Since we are limited by the side lobes, we use a region halfway between the peak and the next side lobe, i.e., shifted from the peak by half the rotation rate. We checked that the shift direction does not matter much, so for the central background frequency, we use the Rossby wave frequencies from LGBD19 and LGBS18 for m = 3 and m ≥ 4, respectively, plus half the rotation frequency ν_{eq}. We use the full widths at half maximum γ^{ref} = Γ^{ref}/2π from LGBD19 and LGBS18 for m = 3 and m ≥ 4, respectively, and perform a leastsquares secondorder polynomial fit in m to obtain a smoothed linewidth γ_{smooth}. We use γ_{smooth} for the width of the peak and background frequency intervals. Thus our peak and background frequency intervals at each m are
These definitions are used in the analysis of latitudinal eigenfunctions in Sect. 3.3.
In Fig. 2, we see that the peak interval (dashed red lines) typically captures the main wave power well. The 1D power spectra, however, reveal that the background interval (not shown in Fig. 3, but see dashed black lines in Fig. 2), however, is potentially contaminated by scattered signal power, for example, for ℓ = m = 5, 6, and 14. To check how the frequency interval definition affects our results, we performed our analysis for several different peak and background intervals. The results are consistent, thus we adopt Eq. (3) for the peak and background intervals.
Unlike LGBS18, we see evidence for nonsectoral components of the Rossby waves. For ℓ = m + 2 the 2D power spectrum shows, for 6 ≤ m ≤ 13, a ridge of power at very similar frequencies to those of the ℓ = m Rossby waves seen in Fig. 2, apart from a higher relative noise level and side lobes. This is confirmed by the 1D cuts at fixed values of m. We do not see structure in the power spectra for ℓ = m + k other than for k = 0 and k = 2. In Sect. 3.3.3, we indeed show that the latitudinal eigenfunctions of Rossby waves are not sectoral spherical harmonic functions (in agreement with LGBS18).
3.3. Latitudinal eigenfunctions of Rossby waves
To estimate the latitudinal eigenfunctions, we first remove smallscale convection from the LCT maps via smoothing with a 6° Gaussian in latitude. Next we compute the Fourier transform of the radial vorticity maps ζ(t, r, λ, φ) from LCT and RDA in time and longitude as follows:
The variables are discrete and take values at time steps t_{j} = jT/N_{t} (integer 0 ≤ j < N_{t}), longitudes φ_{k} = 2πk/N_{φ} (integer 0 ≤ k < N_{φ}), frequencies ν_{s} = s/T (integer s, with −N_{t}/2 ≤ s ≤ N_{t}/2 − 1 for even N_{t}), and azimuthal orders m (integer, with −N_{φ}/2 ≤ m ≤ N_{φ}/2 − 1 for even N_{φ}). In this case, T, N_{t}, and N_{φ} are the observation period and the number of data points in time and longitude, respectively. We apply a filter to select the Rossby waves one m at a time, i.e.,
The filter F_{m}(ν) is equal to one within the Rossby wave ridge and zero elsewhere. Since ζ(t, r, λ, φ) is real, the symmetry applies. We then transform back to time to obtain
In this way we obtain filtered timelatitude vorticity maps for every m. Because there is no symmetry , the filtered vorticity maps are in general complex.
LGBS18 do a similar analysis for LCT data, in particular for rotationaveraged maps and filtering within ±30 nHz around the central mode frequencies. We do the entire latitudinal eigenfunction analysis for LCT and RDA, for full timeresolution maps and maps averaged in time within individual solar rotations, and for a ±27 nHz (five frequency pixels) and a linewidth filter (Eq. (3)) around the central mode frequencies. The different timeresolution and filtering cases yield consistent results; we thus show only the outcome for the full timeresolution and linewidth filtering. However, LGBS18 take the real part of the complex . This is equivalent to assuming that the phase of the eigenfunction is independent of latitude. We address the implications of this in Sect. 3.3.3 in more detail.
To estimate uncertainties for all results in this paper, we split the data into equalsize time intervals, apply our analysis to each chunk, and calculate the standard deviation over the results (for complex quantities separately for the real and imaginary part). Appendix C gives more details on error estimation and validation. Because of the small number of chunks, lownumber statistics are an issue and the reported error bars are relatively uncertain.
For the sake of clarity, for the simple case of a single m Rossby wave with a single frequency ν_{m} and an eigenfunction C_{m}(r, λ), the vorticity field would be given by
We apply two different methods to obtain the eigenfunctions C_{m}(r, λ) near the surface, the covariance method (Sect. 3.3.1), and the SVD method (Sect. 3.3.2). The former is used also by LGBS18.
3.3.1. Covariance
We calculate, at each m, the temporal covariance of the vorticity between the equator near the surface (target depth R = R_{⊙} − 0.7 Mm for RDA) and all other latitudes and depths, normalized by the variance at the equator near the surface
where the angle brackets ⟨ ⋅ ⟩_{t} denote a temporal average and is the centered vorticity. The function C_{m}(r, λ) is complexvalued since is in general complex. By construction C_{m}(r = R, λ = 0° ) is unity. Appendix D shows that C_{m} can also be obtained by a linear fit to the vorticity. The same covariance can be computed with the LCT data.
3.3.2. Singular value decomposition
We present a second, new method to obtain latitudinal eigenfunctions. We want to separate the filtered vorticity at each azimuthal order m and depth r, i.e., a 2D matrix, into a latitude and a time dependence, i.e.,
Applying a singular value decomposition (SVD), we can decompose the vorticity as
where s_{(r, m), j} is the singular value of index j with left and right singular vectors U_{(r, m), j} and V_{(r, m), j} and k is the minimum between the number of grid points in time and latitude. The square of s_{(r, m), j} measures the variance captured by its singular vectors. By convention the singular values are sorted in descending order, thus the first singular vector contains more variance than any other individual singular vector.
Assuming that there is only one nonzero singular value, s_{(r, m), 0}, the SVD gives the desired decomposition of the vorticity into one time and one latitude function. Our observations indeed have one clearly dominant singular value.
Given that the noise at high latitudes increases steeply, we crop our vorticity maps for the SVD to latitudes within ±50° of the equator. Also, the SVD does not account for the varying noise of the remaining latitudes. To ensure that latitudes with larger uncertainties are given less weight, we filter the original vorticity maps once more in Fourier space for the noise, calculate the temporal standard deviation σ_{m} of the noisefiltered maps, and compute . We filter for the noise by taking either all frequencies except for five pixels around the peak or all frequencies within the background interval (see Eq. (3)). The two different filters give consistent results. At each m, the SVD is performed on the weighted maps and the resulting latitude vectors are multiplied by σ_{m} again to undo the weighting. We apply the weighting only to LCT, since the ringdiagram data are already apodized (see Sect. 2). We select the first latitude singular vector near the surface and normalize it by its value at the equator.
3.3.3. Results for the latitudinal eigenfunctions
Figures 4 and 5 show the real and imaginary parts of the horizontal eigenfunctions of Rossby waves versus latitude for different m. The real part is consistent with the findings from LGBS18. The imaginary part, however, was not discussed by LGBS18.
Fig. 4. Real part of C_{m}(λ) for different azimuthal orders m and four different methods (see legend in panel m = 4). The shaded areas indicate the 1σ error estimates. 

Open with DEXTER 
Fig. 5. Imaginary part of C_{m}(λ) for different azimuthal orders m and four different methods (see legend in panel m = 4). The shaded areas indicate the 1σ error estimates. For comparison, the solid gray curves show the real part of C_{m} for the LCT covariancebased data. The plotting ranges are the same as in Fig. 4. 

Open with DEXTER 
In the current paper, we find that the LCT and the RDA results are mostly consistent for the nearsurface layers. Also, almost all m show agreement between the covariance and SVD results. This in particular holds for the modes with the largest amplitudes, i.e., for 7 ≤ m ≤ 10. On the other hand, the modes m = 4 and to a lesser extent m = 15, where Rossby wave measurements become difficult, display larger errors but nonetheless consistent results. The m = 3 results for the different techniques disagree and are noisy. The m = 5 and m = 6 results for the real part differ slightly between the covariance and SVD methods. While the covariance yields a real part of the eigenfunction quite similar to those of other modes, the SVDbased results show maxima around latitudes of ±10 − 15°. Apparently, there the SVD picks up some variance that is uncorrelated with the equator. It is unclear whether it is just noise, or a real signal of a different kind of latitudinal eigenfunctions.
The eigenfunction shape is similar for different modes. The real part decreases away from the equator, flips sign, and approaches zero after going through a local minimum. The imaginary part is much noisier than the real part, as indicated by the error estimates. For most m, it is close to zero and flat near the equator, but reaches minima at high latitudes. The latitude of the minima appears to move toward the equator with increasing m.
As can be seen from, for example, the red curves in Fig. 5, the imaginary part appears to be mostly positive for 3 ≤ m ≤ 6. For 7 ≤ m ≤ 9 the sign of the imaginary part is unclear. For 10 ≤ m ≤ 15, the imaginary part is predominantly negative. The presence of an imaginary part induces a phase for the latitudinal eigenfunctions that can be interpreted as a latitudedependent shift of the sinusoid in longitude. A positive sign of the imaginary part means that the horizontal eigenfunctions at high latitudes are leading in the retrograde direction with respect to the equator. Conversely, a negative sign would indicate that the eigenfunctions at high latitudes are trailing with respect to the equator. This may provide important constraints on the theory of latitudinal eigenfunctions of Rossby waves.
Figure 4 suggests that the real part of the eigenfunctions is more confined to low latitudes for higher values of m. We study the mdependence of several characteristic parameters illustrated in Fig. 6, i.e., the width at an eigenfunction real part of Re(C) = 0.5 (a half width at half maximum; HWHM), the latitude of the eigenfunction real part sign reversal, λ_{0}, and the latitude and value of the minimum, λ_{min} and C_{min}, respectively. To reduce the noise level we derive the eigenfunctions from maps symmetrized in latitude before measuring these parameters.
Fig. 6. Schematic description of the real part of C_{m}(λ) for a given m. The various parameters that describe the curve are the HWHM, the latitude at zero crossing (λ_{0}), the latitude at minimum (λ_{min}), and the minimum value (C_{min}). 

Open with DEXTER 
To obtain the latitude widths at Re(C) = 0.5 and Re(C) = 0, we linearly fit the two points closest to these Re(C) values. The latitude and value of the minimum are obtained by quadratically fitting three points around the minimum derived without fitting. We do this to avoid oscillating RDA results due to the coarse 7.5° latitude sampling. For LCT, the effects of fitting the minimum (or not) are minimal. There are no results for m = 3 and m = 4 because of the poor quality and different shape of the eigenfunctions. We already stated the difficulties in characterizing these modes. As described at the beginning of Sect. 3.3 and in Appendix C, to derive uncertainties, we compute the standard deviation over the results for different time chunks, separately for the real and the imaginary part.
Table 1 shows how these parameters, measured for the LCT data from the covariance method, depend on m. Although not given in the table, we also measure the parameters for the RDA and the SVD results. We thus also discuss the mdependence for those measurements; this dependence is mostly consistent with that of the LCT covariancebased parameters.
The latitude width at Re(C) = 0.5 indeed decreases with m, quasilinearly between m = 7 and m = 13. The slope is roughly −1° per m. The decrease might flatten off at high m, but this could also be caused by noise. We observe slightly different latitude widths between the covariance and SVD eigenfunctions at low m for Re(C) = 0.5, but similar widths at Re(C) = 0. Toward higher m, λ_{0} is consistent with a flat profile, until around m = 13 the eigenfunction widths become smaller. The latitude of the minimum, λ_{min}, shows an m(in)dependence similar to λ_{0}. There is a strong discrepancy for m = 13 between LCT and RDA, indicating that this mode is not trivial to characterize. This could be caused by noise. To some extent we could already see this in the power spectrum in Fig. 3, where the m = 13 linewidth is large compared to all other m. The error on λ_{min} might be underestimated here, since as seen in the asymmetric eigenfunctions in Fig. 4 the minimum of the LCT data is more poorly defined for m = 13 than for other modes. Finally, the value of the minimum, C_{min}, is different between the different analysis methods at m = 5 and m = 6, as seen before. Otherwise, it is quasi mindependent and has at most a very mild increase with m from m = 7 onward, which is likely covered by noise, however.
As mentioned before, the latitudinal eigenfunctions appear to have two nodes (zero crossings) at latitudes ±λ_{0}. This in combination with Fig. 2 and the subsequent discussion indicates that the eigenfunctions have significant contributions from ℓ = m and ℓ = m + 2 components. To quantify these contributions, we project the symmetric eigenfunctions C_{m}(λ) onto associated Legendre polynomials , to obtain the coefficients
The sum goes over all latitudes λ = kπ/N_{λ} (integer −N_{λ}/2 ≤ k < N_{λ}/2), where N_{λ} is the number of data points in latitude. The are normalized such that . The associated Legendre polynomials used in the decomposition are not orthogonal over the limited observed latitude range. However, we do not expect this to be a problem since we see later in this section that the nearsectoral associated Legendre polynomials, whose amplitude is concentrated toward the equator, are the dominant contributions to the latitudinal eigenfunctions. Because of the symmetry of the eigenfunctions, only c_{ℓm} with even ℓ − m ≥ 0 are nonzero. We find that the real and the imaginary parts of the eigenfunctions for almost all m can be approximated well (within 1σ) when using only the contributions from c_{ℓm} for m ≤ ℓ ≤ m + 6, except for m = 3, which is very noisy. The approximation also does not work well at the high latitudes (beyond ±40°) for the real part (for some modes) and at the nearequatorial latitudes for the imaginary part (for almost all modes).
Table 2 shows the coefficients c_{ℓm} for m ≤ ℓ ≤ m + 6 for the LCT covariancebased latitudinal eigenfunctions. As usual the uncertainties are calculated from the standard deviation over the coefficients for different time chunks (Appendix C), separately for the real and the imaginary part. The real part of the eigenfunctions is clearly dominated by the ℓ = m component. The contribution from the ℓ = m + 2 component is significant as well and has a relative strength of 30−50%. This is consistent with the observations from the 2D and 1D power spectra in Sect. 3.2. The real part of the c_{mm} and c_{m+2, m} each depend weakly on m. The real part of several of the coefficients with larger ℓ is insignificant. The imaginary part, on the other hand, has significant, dominant contributions at ℓ = m + 4 for m ≤ 10 and at ℓ = m + 2 for m ≥ 11, whereas the ℓ = m and ℓ = m + 6 components are often insignificant. The term insignificant refers to an absolute value of c_{ℓm} of less than 1σ. Nonetheless, independent of the estimated error bars, 11 out of 12 modes within 4 ≤ m ≤ 15 have the same sign for c_{m + 4, m}, suggesting that the ℓ = m + 4 contribution to the real part is significant, despite being within 1σ from zero. A similar argument holds for the imaginary part of the latitude dependence of the Rossby wave eigenfunctions.
Coefficients c_{ℓm} for the LCT covariancebased data.
For the latitudinal eigenfunctions of Rossby waves there are so far only a few theoretical studies such as Lee & Saio (1997) and Townsend (2003). These studies typically gave either analytic (asymptotic) expressions and/or numeric calculations, but the former expressions do not agree well with the latter calculations for Rossby waves (Townsend 2003). Although both studies indicate that the latitudinal eigenfunctions are not concentrated near the equator, we cannot sensibly compare their findings to our measurements. In particular these models assume a uniform rotation rate. Also these authors used the traditional approximation, i.e., they neglected the horizontal component of the rotation vector. This approximation requires the squared BruntVäisälä frequency N^{2} to be much higher than both the squared oscillation frequency ω^{2} and the squared rotation rate Ω^{2}. The validity of the traditional approximation thus has to be critically examined within the convection zone of the Sun.
3.4. Radial eigenfunctions of Rossby waves
3.4.1. Depthdependent ringdiagram systematics
To study the Rossby wave depth dependence, we must check to which depths RDA is reliable. For this we compare the solar rotation profile from ringdiagram velocities with the results from SDO/HMI global modes from the JSOC data series hmi.V_sht_2drls (Larson & Schou 2018). The global modes have a 72day time sampling from April 30, 2010 to June 4, 2017, a 1.875° latitude sampling, and a nonlinear depth grid with many more points near the surface than at larger depths. Global modes are expected to give a precise and accurate solar rotation profile. We interpolate the global mode results to the ringdiagram latitudedepth grid via 2D bicubic splines, which is reasonable because the globalmode inversions do not vary on scales of their original grid; we then average the 72day chunks over time. The chunk scatter of the rotation rate is used to estimate the uncertainty. We divide the ringdiagram data into five intervals of length 480 time steps (20 rotations), average the chunks over time and estimate the error from the scatter, convert the velocities into rotation rates, and add the sidereal Carrington rate to correct for the ringdiagram tracking.
Figure 7 shows the equatorial rotation rate versus depth from global modes and ringdiagram velocities, both averaged over longitude and at Stonyhurst longitudes of ±52.5° (the outermost longitudes in our vorticity maps). The global modes yield a smooth profile with extremely small errors. The ringdiagram data show a small offset at small depths, but, more importantly, inconsistency with the global modes at large depths. Of course, it is difficult to judge how well the results should agree because of the different kernels of the datasets and thus different depth (and latitude) sensitivities. The −52.5° longitude curve has a small local maximum around 8 Mm. For longitudes even further east (not shown), the rotation rate has an even stronger excess (a bump) there.
Fig. 7. Solar equatorial rotation rate as a function of depth. The globalmode helioseismology result is given by the black curve. The blue curve is for RDA after averaging over all longitudes. The green and red curves show the ringdiagram results at Stonyhurst longitudes ±52.5°. The shaded areas give the 1σ error estimates. The observations cover more than seven years. The dashed black lines indicate the depth range within which the ringdiagram results are in best agreement with each other. 

Open with DEXTER 
The most worrisome point is the disagreement between different ringdiagram longitudes themselves and also with the longitude average, below roughly 8 Mm (indicated by the left dashed line in Fig. 7). Because we averaged the data over more than seven years, any shortlived flows and even longerlived structures should be filtered away and the longitude gradient from east to west should thus not exist. This points to a deeper problem with the ring fits and the pipeline processing that generated these fits. The presence of systematic effects in HMI ringdiagram data has also been extensively discussed in Komm et al. (2015).
Finally, we note that Fig. 7 is affected by an issue related to the ringdiagram inversion, since the inversion does not account for the quantity β_{nℓ}. A discussion of this issue and a brief check of the magnitude of the effect is given in Appendix A. The latter showed that the main effect is a depthindependent underestimation of the ringdiagram velocities by 1 − 2 m s^{−1} or equivalently of the rotation rate by less than 0.5 nHz. This does not affect our main conclusions. The small, but significant difference between the rotation rates from global modes and ring diagrams cannot be caused by the β_{nℓ} issue (it has the wrong sign), but may possibly instead be due to different averaging kernel widths, systematics, or other unknown effects.
3.4.2. Determining the Rossby wave depth dependence
In this section, we discuss only the sectoral (ℓ = m) component of the power spectrum of radial vorticity. The Rossby wave power P_{m}(ν, r) and phase Φ_{m}(ν, r) thus depend on frequency, depth, and azimuthal order. Based on the assumption of damped oscillations and stochastic wave excitation, we perform a maximumlikelihood Lorentzian fit (Anderson et al. 1990) to the power spectra for the longer ringdiagram period, separately at each m. We use the functional form
We fit all the depths (except for the surface, i.e., r = 0.0 Mm, where the ringdiagram data are unreliable) at once, with a common central frequency ν_{0, m} and linewidth γ_{m}, but with individual amplitudes A_{m}(r) and backgrounds B_{m}(r). The Lorentzian fit of the power spectra, in most cases, fits well to the observations. As seen in Fig. 3, the ℓ = m = 6 and ℓ = m = 13 modes have large linewidths and their power spectra show fine structure. The ℓ = m = 3 mode has been fit by LGBD19, but not by LGBS18.
To determine error bars for the amplitudes and backgrounds via chunked data (Appendix C), we also perform the Lorentzian fit for each chunk separately, fitting again all depths together, but keeping the central frequency and linewidth fixed at the fit results of ν_{0, m} and γ_{m} from Eq. (12) for the full time period (to prevent unstable fits). Because we keep these parameters fixed for each chunk, we cannot derive their uncertainties based on the standard deviation over the chunks. We thus do a Monte Carlo simulation and generate 1000 realizations of synthetic power spectra according to Eq. (C.1) (Appendix C) and perform the Lorentzian fit for each realization analogously to the fit for the observations. The median of the parameters over the Monte Carlo realizations is consistent with the fit parameters for the observations. While the error based on the Monte Carlo simulation contains realization noise, the model we use (Lorentzian and stochastic excitation) does not include all features of the observed power spectra. The chunkbased error likely describes the physical system more accurately, by also including other variance contributions, such as from temporal effects on the Rossby waves; we could imagine, for example, solar cycle effects. This may also explain the discrepancy between the two types of errors of order 30% (Appendix C). This disagreement is, however, small enough to not affect the significance of the results for the radial eigenfunctions. We thus use the uncertainties based on the Monte Carlo simulation for the central frequency and the linewidth.
Table 3 compares the fit parameters from this study with the results from LGBD19 and LGBS18. As in LGBD19 and LGBS18, the upper and lower errors give the difference between the quantiles comprising the central 68.3% (the distributions are nonGaussian) and the fit parameters for the observations. We also calculate the uncertainties on the central frequency following Libbrecht (1992) and find that those (symmetric) errors typically underestimate the Monte Carlo quantile errors by roughly 1 nHz. A possible reason for this could be that we use a finite frequency fitting interval.
Measured frequencies and linewidths of the Rossby waves from RDA sectoral power spectra with azimuthal orders in the range 3 ≤ m ≤ 15.
The fit parameters for the observations and those from LGBD19 and LGBS18 typically agree within 1σ or better. The central frequencies and the linewidths for the ℓ = m = 6 mode differ by 10 and 37 nHz, but the fit is sensitive to the fitting range. The ℓ = m = 3 fit parameters do not agree. In LGBD19, the authors, using 21 years of data, observed that the multipeak structure of the ℓ = m = 3 power spectrum (Fig. 3) seen in data with shorter periods collapses to a narrow single peak, which indicates that the discrepancy of the fit parameters is explained by stochastic excitation of the Rossby waves and vanishes when fitting data with a longer time period. Our errors are typically more symmetric and often smaller than those of LGBD19 and LGBS18. The lower errors for the linewidth often agree better than the upper errors, indicating a tail of high values (skewness) in the LGBD19 and LGBS18 estimate distributions. Reasons for the differences in the error estimates may lie in the simultaneous fitting of all depths at once or in the different observation periods of our datasets and those of LGBS18.
To determine the power depth dependence, we use the amplitude A derived from the Lorentzian fit (see Eq. (12)) and we define the normalized power of the signal as
We thus normalize by the depth average of the amplitude of the Lorentzian. The normalized power is independent of temporal amplitude variations due to Rossby wave excitation.
3.4.3. Results for the radial eigenfunctions
Figure 8 shows the depth dependence of the ℓ = m = 8 phase, but the behavior is similar for other m. For easier comparison, we remove phase jumps of 360° and move the depth average to zero. The phase at the frequency of maximum power is almost constant with depth, within roughly ±3°. The phase at the background, at the center of the background interval, varies much more strongly with depth, within roughly ±100°, although the phases at other background frequencies sometimes show much less variation. The background phase is not random in depth (see Appendix C). In particular, phase changes are gradual and smooth; the depths are correlated. This could indicate a significant contribution from scattered signal power to the background. Nonetheless, peak and background display distinctly different depth dependences. We also find that phases at different frequencies across the peak and background are different. Frequency averages of phases are thus not useful. However, as seen, for single frequencies the phase at the peak is nearly constant with depth, while the background phase varies with depth.
Fig. 8. Phase at a frequency of roughly −87.4 nHz corresponding to the peak of power for the Rossby mode ℓ = m = 8, as a function of depth (blue curve). The green line shows the phase of the background at the center of the background interval (see Eq. (3)). The depth refers to the median of the ringdiagram averaging kernels (orange dots), which corresponds to certain target depths (black dots). 

Open with DEXTER 
Figure 8 also shows the main parameters of the ringdiagram averaging kernels for a few target depths, i.e., the first, second (median) and third quartiles and the width (interquartile range). The flow measurements are welllocalized near the surface, but smeared out over a broad depth range at large depths. The ringdiagram depth covariance matrix (not shown) indicates a similar behavior and shows that different depths are mostly independent near the surface, while at the largest depths there is high correlation and they thus do not give independent results. This could maybe also explain why the background phase is not random in depth. At large depths, the center of the averaging kernels (second quartile) moves away from the target depth, but the averaging kernels are relatively symmetric.
The background power for different m (not shown) generally increases with depth and at least for some modes there could be a minimum at 8 − 9 Mm, albeit with little significance given the large errors.
Figure 9 shows the signal power (Eq. (13)). The quantity 𝒫_{signal} typically decreases from the surface toward a depth of 8 Mm, significantly as shown by the errors. Even further inside the Sun the power often increases again and reaches nearsurface or even higher values. The 1σ errors shown in this plot give the standard deviation, but they do not indicate 68.3% probability intervals, since the power distribution is nonGaussian (power cannot be negative). More information about error estimation can be found in Appendix C.
Fig. 9. Blue lines show the Rossby wave power 𝒫_{signal} as a function of depth for different values of m. The blue shaded areas indicate the 1σ errors. The orange curves are fits of the form const. × r^{2α} over depths between 0.7 and 7.4 Mm (between the vertical dashed lines). The orange dashed curves are extrapolations to larger depths. 

Open with DEXTER 
Provost et al. (1981) presented a theoretical argument that the Rossby wave eigenfunctions for the horizontal displacement are proportional to r^{m} under the assumption that the modes are incompressible. Thus, under this theory, the radial vorticity is expected to be proportional to r^{m−1}. To compare this to our observations, we perform a fit of the form const. × r^{2α} to 𝒫_{signal} within the dashed black lines (0.7 to 7.4 Mm) where the RDA is more reliable (see Fig. 7). We assume that the data points are uncorrelated in depth. Obviously, the fit does not reproduce the increase of power at large depths.
Figure 10 compares the observed and theoretical exponent α. The fitted exponent has very large error bars. It is consistent with the theoretical model from Provost et al. (1981), but also with the absence of any trend with m. Although the exponent depends strongly on the fit range because of the kink at roughly 5 Mm in Fig. 9, we also do not find inconsistency with a flat dependence on m within other fit intervals. Thus the current error estimates do not allow a definitive statement on the radial dependence of Rossby waves.
Fig. 10. Exponent α as a function of m, measured in the top 7.4 Mm (blue line) and 1σ error (blue shaded area). The dashed black line corresponds to the model α = m − 1, obtained under the assumption of nondivergent motions (Provost et al. 1981). 

Open with DEXTER 
4. Summary
We build on LGBS18, who investigated Rossby waves mostly using granulation tracking, by studying several Rossby wave properties via the analysis of radial vorticities computed from RDA at different depths and LCT at the surface. We obtained several new results: independently the latitudinal eigenfunctions with RDA (including a more complete, complexvalued description of the eigenfunctions), and the Rossby wave power and phase depth dependence.
We calculated latitudinal eigenfunctions of Rossby waves from the radial vorticity maps via the covariance between the equator and different latitudes and from the singular vectors of an SVD. We confirmed the shape of the real part of the eigenfunction from LGBS18, who used the covariance method on symmetrized LCT data. We also saw consistency between covariance and SVD results, except for m = 5 and m = 6, where the SVD eigenfunctions had maxima around ±10−15° instead of at the equator. The shape of the real part of the latitudinal eigenfunctions seen for most m indicates that the Rossby waves have maximum amplitudes near the equator, as found by LGBS18. The imaginary part appears to be mostly positive for low m (3 ≤ m ≤ 6); this part varies around zero for intermediate m (7 ≤ m ≤ 9) and is mostly negative for high m (10 ≤ m ≤ 15). A nonzero imaginary part may be due to attenuation and to the interaction of the waves with largescale flows. In particular, the interaction of viscous Rossby waves with latitudinal differential rotation leads to the formation of critical layers at intermediate latitudes (Gizon & Fournier, priv. comm.).
We defined and measured characteristic parameters for the real part of the eigenfunctions and we found that the width at an eigenfunction value of 0.5 (the HWHM) decreased with m, in contrast to the mindependent width at a value of 0 and the latitude and value of the eigenfunction minimum. We also decomposed the eigenfunctions into associated Legendre polynomials and saw that the real part is dominated by ℓ = m and ℓ = m + 2 contributions, while the imaginary part consists mostly of ℓ = m + 4 and ℓ = m + 2 contributions for low and high m, respectively.
We compared rotation rates from ringdiagram data and global modes as functions of depth and saw a small offset at small depths and disagreement at large depths, but most importantly inconsistency of different ringdiagram longitudes. This indicated systematic effects in the ringdiagram pipeline (see also Komm et al. 2015).
We studied the Rossby wave power and phase depth dependence in detail for the first time. The phase at the peak is stable with depth, in contrast to the phase of the background. The background power almost monotonically increases with depth, while the signal power decreases toward a depth of 8−9 Mm and then increases again. The radial eigenfunctions of the Rossby waves are (at small depths) consistent with a powerlaw decrease, in particular both with the theoretical Provost et al. (1981) model (exponent m − 1) and an mindependent exponent. However, the Provost et al. (1981) model is based on assumptions that are not exactly correct for the Sun (e.g., uniform rotation). We can constrain the radial dependence of the eigenfunctions only very weakly owing to the high uncertainties on the observed exponents.
The analysis presented in this paper implicitly makes the assumption that the Rossby wave eigenfunctions are separable in depth and latitude. Our data show a similar latitude dependence for the different depths, separability thus appears to be a reasonable assumption. The results shown in this work motivate further research on Rossby wave eigenfunctions, which is a necessary condition for the interpretation of the measured mode frequencies.
Acknowledgments
B. Proxauf is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen; he conducted the data analysis and contributed to the interpretation of the results and to the writing of the manuscript. We thank Z.C. Liang for providing help with the fitting of the modes in frequency space and V. Böning for providing beta values for Appendix A. The HMI data are courtesy of NASA/SDO and the HMI Science Team. The data were processed at the German Data Center for SDO funded by the German Aerospace Center DLR. We acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218.
References
 Aerts, C., ChristensenDalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Springer) [Google Scholar]
 Alshehhi, R., Hanson, C. S., Gizon, L., & Hanasoge, S. 2019, A&A, 622, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anderson, E. R., Duvall, Jr., T. L., & Jefferies, S. M. 1990, ApJ, 364, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Baldner, C. S., & Schou, J. 2012, ApJ, 760, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & RabelloSoares, M. C. 2011a, J. Phys. Conf. Ser., 271, 012008 [Google Scholar]
 Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & RabelloSoares, M. C. 2011b, J. Phys. Conf. Ser., 271, 012009 [Google Scholar]
 Bogart, R. S., Baldner, C. S., & Basu, S. 2015, ApJ, 807, 125 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2008, Ap&SS, 316, 113 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2011, Astrophysics Source Code Library [record ascl:1109.002] [Google Scholar]
 Duvall, Jr., T. L., Jefferies, S. M., Harvey, J. W., & Pomerantz, M. A. 1993, Nature, 362, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Fisher, G. H., & Welsch, B. T. 2008, in Subsurface and Atmospheric Influences on Solar Activity, eds. R. Howe, R. W. Komm, K. S. Balasubramaniam, & G. J. D. Petrie, ASP Conf. Ser., 383, 373 [NASA ADS] [Google Scholar]
 Haber, D. A., Hindman, B. W., Toomre, J., et al. 2000, Sol. Phys., 192, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S., & Mandal, K. 2019, ApJ, 871, L32 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Woodard, M., Antia, H. M., Gizon, L., & Sreenivasan, K. R. 2017, MNRAS, 470, 1404 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, F. 1988, ApJ, 333, 996 [NASA ADS] [CrossRef] [Google Scholar]
 Komm, R., González Hernández, I., Howe, R., & Hill, F. 2015, Sol. Phys., 290, 1081 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
 Liang, Z.C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Libbrecht, K. G. 1992, ApJ, 387, 712 [NASA ADS] [CrossRef] [Google Scholar]
 Lisle, J., & Toomre, J. 2004, in SOHO 14 Helio and Asteroseismology: Towards a Golden Future, ed. D. Danesy, ESA Special Publication, 559, 556 [NASA ADS] [Google Scholar]
 Löptien, B., Birch, A. C., Duvall, T. L., Gizon, L., & Schou, J. 2016, A&A, 590, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löptien, B., Birch, A. C., Duvall, T. L., et al. 2017, A&A, 606, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [NASA ADS] [CrossRef] [Google Scholar]
 Mandal, K., & Hanasoge, S. 2019, ArXiv eprints [arXiv:1908.05890] [Google Scholar]
 November, L. J., & Simon, G. W. 1988, ApJ, 333, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 182, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Provost, J., Berthomieu, G., & Rocca, A. 1981, A&A, 94, 126 [NASA ADS] [Google Scholar]
 Saio, H. 1982, ApJ, 256, 717 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Smeyers, P., Craeynest, D., & Martens, L. 1981, Ap&SS, 78, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D. 2003, MNRAS, 340, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Welsch, B. T., Fisher, G. H., Abbett, W. P., & Regnier, S. 2004, ApJ, 610, 1148 [NASA ADS] [CrossRef] [Google Scholar]
 Wolff, C. L., & Blizard, J. B. 1986, Sol. Phys., 105, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Woodard, M. F. 1989, ApJ, 347, 1176 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, J., Nagashima, K., Bogart, R. S., Kosovichev, A. G., & Duvall, Jr., T. L. 2012, ApJ, 749, L5 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Issues of the ringdiagram inversions
In this Appendix, we discuss two issues regarding the ringdiagram pipeline inversions. In order to obtain local velocities at a certain measurement depth r, the reported pipeline velocities must be multiplied by r/R_{⊙}. Additionally the pipeline inversion does not take the quantity β_{nℓ} (see, e.g., Aerts et al. 2010, Eq. (3.357)) into account and thus the reported inversion velocities u_{x} are slightly incorrect.
To see this, we study a simple case. For now, let us assume that the ring diagrams are not tracked. The frequency perturbation δω_{nℓm} of the mode indexed by radial order n, angular degree ℓ, and azimuthal order m due to a radial differential rotation rate Ω(r) is
where K_{nℓ} is the normalized rotation kernel for that mode, i.e., (see, e.g., Aerts et al. 2010, Eq. (3.358)). On the other hand, ring diagrams assume that the velocity mode fits U_{x, nℓ} are equal to a radial integral over the true velocity flow field u_{x}(r) weighted by flow sensitivity kernels. Based on inspection of the pipeline, we think that the used HMI kernels are normalized rotation kernels K_{nℓ} from Eq. (A.1). Thus
To connect the two equations in a simple case, consider the Doppler shift of a sectoral (ℓ = m) mode as seen by a ring diagram at the equator, i.e.,
In this equation, k_{x} is the wavenumber in the prograde direction, which is related to m via k_{x} = m/R_{⊙}. We conclude that
This is not consistent with Eq. (A.2) since β_{nℓ} is missing from Eq. (A.2). Additionally we see that u_{x}(r) should be interpreted as R_{⊙}Ω(r) and not as the local linear velocity rΩ(r).
To see what happens if the tracking rate is not zero, we now suppose that we track at rotation rate Ω_{T}. Equation (A.4) then becomes
where is the ring measurement in the rotating frame. We now define the local deviation from the tracking rate δΩ(r) = Ω(r)−Ω_{T}. Then we obtain
The first term is the same form as in Eq. (A.4), while the second term is an offset that depends on n and ℓ.
The conversion factor r/R_{⊙} is multiplied onto the data before any analysis is performed for this paper; see Sect. 2.2. The offset due to β_{nℓ} depends on the set of mode ring fits, but it should be mostly timeindependent, since the ringdiagram mode set does not vary much with time. Thus the timedependent Rossby waves should not be sensitive to this effect and the only affected result in this paper should be the comparison of rotation rates in Fig. 7.
To estimate the effect of β_{nℓ} on the inversion result, for a given input flow u_{x} we generate artificial ring fits via Eq. (A.6), on which we run the ringdiagram inversion module to retrieve the output velocities. To compute , we assume a depthindependent flow equivalent to the tracking rate (sidereal Carrington rate), thus δΩ(r) = 0. We thus check only the second term of Eq. (A.6) and neglect that β_{nℓ} also appears in the first term as a scaling factor. However, the effect due to the second term should be much larger than that due to the first term, as Ω_{T} is much larger than δΩ(r) for the ring diagrams.
We use β_{nℓ} values provided by V. Böning (priv. comm.). These were computed from eigenfunctions obtained from the Aarhus adiabatic oscillation package (ADIPLS; ChristensenDalsgaard 2008, 2011). We lose roughly 25% of the original ringfit modes, as we only have β_{nℓ} values up to frequencies of 5 mHz. However, this does not critically change the mode set used during the inversion. We replace the actual pipeline ring fits with the artificial data. We leave all other data, including uncertainties on modefit velocities, as is and perform the inversion. The aforementioned conversion factor of r/R_{⊙} is multiplied onto the output velocities u_{x}.
We see that the effect of β_{nℓ} does not depend much on depth and that the retrieved u_{x} are on the order of only 1.5 m s^{−1} (equivalent to roughly 0.1% of the tracking rate). The reason for this is that the inversion gives much more weight to the high ℓ modes for which the uncertainties are comparatively small. These modes typically have β_{nℓ} values around 0.999, thus 1 − β_{nℓ} ∼ 0.1%. We performed this check exemplarily for a ringdiagram tile at the first time step in our dataset (May, 20, 2010) at the point (λ = 0° ,φ = 0°). However, tests using different tiles show that this result does not depend much on time or disk position.
The effect of the pipeline inversion not accounting for β_{nℓ} is thus an underestimation of the true velocity fields by roughly 1 − 2 m s^{−1}, or equivalently approximately 0.4 nHz. This difference would be visible in Fig. 7, but does not change our main conclusions.
Appendix B: Interpolation and apodization of ringdiagram velocities
We interpolate the ringdiagram velocities separately in time and longitude (see Sect. 2) with different functions, depending on the number of available data points:
Before we interpolate the ringdiagram velocities to the surface equatorial rotation rate, we apodize these velocities with a raised cosine in angular greatcircle distance ρ to the point (λ = 0° ,φ = 0°), see Sect. 2), as follows:
where β defines the steepness of the raised cosine flanks. We choose β = 0.3. The quantity T defines the central position of the flanks. We choose T such that zero is reached at ρ = 67.5° (where there are no more valid pixels). Apodizing the ringdiagram velocities (with different β), or not, gives consistent results.
Appendix C: Error estimation and error validation
C.1. Error estimation via chunked data
The uncertainties on all results are derived by dividing the time series of vorticity maps (in total 2448 time steps, i.e., 102 rotations, for RDA) into equalsize chunks and calculating the scatter over the results for the chunks. We find that for chunks longer than a few months (roughly six rotations), the Rossby wave signature is visible in the power spectra. We make a compromise between noise level and chunk statistics and divide the dataset into five chunks that are 480 time steps long each (20 rotations).
For the latitudinal eigenfunctions, where we only study the shorter LCT period (i.e., 78 rotations), we first used four chunks of length 19 rotations (rotationaveraged maps) and 470 time steps (full time resolution, for RDA), but we obtained very large errors for the SVD method for different single m, where single chunks gave singular vectors different from the usual eigenfunction shape. We think that the noise in our filtered maps might have been detected as the dominant term in the decomposition. We thus use a chunking with three chunks of length 26 rotations and 625 time steps, where all chunks have the expected first singular vectors.
C.2. Error validation via Monte Carlo simulation
We validate the chunking approach for the depth dependence via a Monte Carlo simulation. As a plausible physical model for the Rossby wave power spectrum, we assume a Lorentzian profile and a background (constant in frequency), each with a χ^{2}distributed random variable (stochastic excitation). In analogy to Eq. (12) we generate 1000 realizations of synthetic data for the Fourier transform of the radial vorticity ℱ_{syn} (not the power P), at each m, as
We fix the amplitude A_{m}(r), background B_{m}(r), central frequency ν_{0, m}, and full width at half maximum γ_{m} via the fit parameters for the observations. Furthermore we assume that the random variable 𝒩_{A, m}(ν) is constant with depth, i.e., the signal is fully correlated in depth, while for the background, we take a random variable 𝒩_{B, m}(ν, r) that is uncorrelated in depth. The random variables are complex Gaussian variables (with independent real and imaginary parts) with zero mean and unit variance, independent for each frequency. We analyze the realizations in the same way as the observations.
We observe inconsistency between two Monte Carlo estimates (by roughly 30%): one from chunking (like the observations), averaged over the realizations, and the other from the scatter of the power over the realizations. We find that the discrepancy is due to the temporal correlation of the different chunks; the temporal correlation matrix has values around 0.1 on the first offdiagonal. The two quantities agree when using a weighted average with weights based on the temporal covariance matrix.
Both Monte Carlo estimates disagree with the error for the observed data. This could be due to the depth correlation of the observations. We determine a (noisy) estimate of the depth covariance and correlation matrix of the observed background power from the different background frequencies and find strong correlations between different depths, even between the largest depth and the surface, with values above 0.25 in the offdiagonal corners. We use the observed depth covariance matrix for the Fourier transform (not the power) to correlate the Gaussian background random variable of our realizations, via a Cholesky decomposition.
Correlating the background in depth noticeably improves the agreement between Monte Carlo and observed errors. However, there is still a remaining discrepancy that can be attributed to our model, which is missing some features of the observed power spectrum. We do not account for the window function nor a background decreasing with frequency present in the observations. A detailed analysis goes beyond the scope of this paper.
Appendix D: Relation of covariance to linear fit
The covariance method (Sect. 3.3.1) is conceptually equivalent to a linear fit of the vorticity at each depth r and latitude λ, i.e.,
Let us for simplicity assume that the vorticity is real. The latitude and depth dependence in this vorticity separation ansatz is contained in the fit parameter a_{m}(r, λ). Let us assume that only is uncertain. For zeromean quantities (such as our vorticity maps ), the slope of a linear fit without intercept, i.e., a_{m}(r, λ), is given as
Assuming the time dependence is given by the surface equatorial vorticity time series, i.e., , in Eq. (8) we can identify a_{m}(r, λ) with C_{m}(r, λ). Equation (D.1) implies that a_{m}(r = R, λ = 0° ) is unity.
The main disadvantage of the covariance method is the assumption of a noisefree vorticity at the equator, , required so that the timedependence f_{m}(t) is noisefree and the vorticity is the only uncertain quantity of the fit.
All Tables
Measured frequencies and linewidths of the Rossby waves from RDA sectoral power spectra with azimuthal orders in the range 3 ≤ m ≤ 15.
All Figures
Fig. 1. Radial vorticity maps from LCT at the surface and from RDA at depths 0.7, 9.9, and 16.5 Mm. The radial vorticity is averaged over one rotation (from May 20, 2010 to June 16, 2010). The LCT map is smoothed in latitude and longitude with a Gaussian filter (σ = 6°) to filter out smallscale convection. 

Open with DEXTER  
In the text 
Fig. 2. Sectoral power spectrum (ℓ = m) of the radial vorticity for RDA data at depth 0.7 Mm. The solid red line indicates the Rossby wave frequencies from LGBD19 for m = 3 and from LGBS18 for m ≥ 4. Frequency intervals for the Rossby wave region and the background region are indicated by the red and black dashed lines, respectively. The ridge of power at positive frequencies is due to the first side lobe of the window function. For better visibility of the Rossby waves at low m, the power is normalized at each m by the frequency average over [−300, 100] nHz. The color scale is truncated at 50% of the maximum value (black). 

Open with DEXTER  
In the text 
Fig. 3. Sectoral power spectra (ℓ = m) showing the Rossby wave power in the LCT data (black line) and the RDA data at depths 0.7 and 16.5 Mm (cyan and red lines). The power is normalized by the average of the m = 8 power in the range [−300, 100] nHz. The dashdotted vertical lines in the m = 1 and m = 2 panels indicate the frequencies ω = −2Ω_{eq}/(m + 1). The frequency axes of the m = 1 and m = 2 panels are shifted with respect to the other panels. 

Open with DEXTER  
In the text 
Fig. 4. Real part of C_{m}(λ) for different azimuthal orders m and four different methods (see legend in panel m = 4). The shaded areas indicate the 1σ error estimates. 

Open with DEXTER  
In the text 
Fig. 5. Imaginary part of C_{m}(λ) for different azimuthal orders m and four different methods (see legend in panel m = 4). The shaded areas indicate the 1σ error estimates. For comparison, the solid gray curves show the real part of C_{m} for the LCT covariancebased data. The plotting ranges are the same as in Fig. 4. 

Open with DEXTER  
In the text 
Fig. 6. Schematic description of the real part of C_{m}(λ) for a given m. The various parameters that describe the curve are the HWHM, the latitude at zero crossing (λ_{0}), the latitude at minimum (λ_{min}), and the minimum value (C_{min}). 

Open with DEXTER  
In the text 
Fig. 7. Solar equatorial rotation rate as a function of depth. The globalmode helioseismology result is given by the black curve. The blue curve is for RDA after averaging over all longitudes. The green and red curves show the ringdiagram results at Stonyhurst longitudes ±52.5°. The shaded areas give the 1σ error estimates. The observations cover more than seven years. The dashed black lines indicate the depth range within which the ringdiagram results are in best agreement with each other. 

Open with DEXTER  
In the text 
Fig. 8. Phase at a frequency of roughly −87.4 nHz corresponding to the peak of power for the Rossby mode ℓ = m = 8, as a function of depth (blue curve). The green line shows the phase of the background at the center of the background interval (see Eq. (3)). The depth refers to the median of the ringdiagram averaging kernels (orange dots), which corresponds to certain target depths (black dots). 

Open with DEXTER  
In the text 
Fig. 9. Blue lines show the Rossby wave power 𝒫_{signal} as a function of depth for different values of m. The blue shaded areas indicate the 1σ errors. The orange curves are fits of the form const. × r^{2α} over depths between 0.7 and 7.4 Mm (between the vertical dashed lines). The orange dashed curves are extrapolations to larger depths. 

Open with DEXTER  
In the text 
Fig. 10. Exponent α as a function of m, measured in the top 7.4 Mm (blue line) and 1σ error (blue shaded area). The dashed black line corresponds to the model α = m − 1, obtained under the assumption of nondivergent motions (Provost et al. 1981). 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.