Exploring the latitude and depth dependence of solar Rossby waves using ring-diagram analysis

Context. Global-scale 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 ring-diagram 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 ring-diagram 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 non-sectoral 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 ring-diagram 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 ring-diagram 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.


Introduction
Recently, Löptien et al. (2018, hereafter LGBS18) discovered global-scale 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 large-scale 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 time-distance 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 ring-diagram analysis (RDA; Hill 1988) via machine learning, also saw global-scale Rossby waves.  and  provide another recent Rossby wave confirmation using a different technique of helioseismology known as normal-mode 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).

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, 2010to March 30, 2016. The second dataset comprises RDA flow maps from the HMI ring-diagram 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, 2010to March 31, 2016, while for all other results we used a longer period of more than seven years from May 19, 2010to December 29, 2017; this corresponds to 102 Carrington rotations (CRs), i.e., CR 2097 -2198.

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 shrinking-Sun 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 equi-spaced longitudelatitude grid with a step size of 0.4 • in both directions.

Overview of ring-diagram data
The ring-diagram 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 latitude-dependent. 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 ring-fit 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 six-parameter 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 ring-diagram 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 re-inverted 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 ring-diagram pipeline that have not yet been solved. Among these are under-regularization 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 ring-diagram 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.

Post-processing of ring-diagram data
The ring-diagram 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 ring-diagram velocities, such as center-to-limb effects that depend on the disk position of the tile (Baldner & Schou 2012;Zhao et al. 2012). There are time-independent 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 u x (t) = a x sin(2πt/(1 yr)) + b x cos(2πt/(1 yr)) + c x , u y (t) = a y sin(2πt/(1 yr)) + b y cos(2πt/(1 yr)) + c y , 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 ring-diagram 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 latitude-independent longitude grid, we interpolated the flow maps in Stonyhurst longitude using splines (Appendix B).
We also interpolated the ring-diagram 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 A44, page 2 of 16 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 ring-diagram pipeline increase strongly toward the limb, in particular beyond an angular great-circle distance of roughly 65 • to the crossing of the central meridian with the equator (λ = 0 • , ϕ = 0 • ). We thus only used ring-diagram data within 65 • of (λ = 0 • , ϕ = 0 • ).

From velocity maps to power spectra of radial vorticity
From this stage onward ring-diagram 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 global-mode analysis of SDO/HMI observations (Larson & Schou 2018). We shifted the LCT data in Fourier space via a time-dependent phase factor, applying the same convention for the Fourier transform as LGBS18. The ring-diagram data are first apodized by a raised cosine in angular great-circle distance to the point (λ = 0 • , ϕ = 0 • ) to suppress near-limb data and are shifted via splineinterpolation (Appendix B).
We next subtracted the longitude mean from the data to remove any remaining large-scale flows. Differential rotation and meridional circulation should have already been subtracted in the RDA or LCT post-processing, 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 half-open 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.  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 small-scale 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.

Radial vorticity maps
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 near-surface (0.7 Mm) ring-diagram map. While large absolute radial vorticities are visible at high latitudes (beyond ±50 • ) in the LCT but not in the ring-diagram 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 ring-diagram 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. Figure 2 shows the Rossby wave power of the = m component for the ring-diagram 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. 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 time-distance 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 multi-peak structure and is thus difficult to measure. We do not observe Rossby waves for = m ≤ 2; the dash-dotted blue lines for = m = 1 and = m = 2 indicate the expected mode frequencies from the textbook dispersion relation.

Power spectra of radial vorticity
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 ν ref 0 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 LGBS18 for m = 3 and m ≥ 4, respectively, and perform a leastsquares second-order 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 peak interval : 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 non-sectoral 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).

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 , 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 ζ m (ν, r, λ) =ζ * −m (−ν, r, λ) applies. We then transform back to time to obtaiñ In this way we obtain filtered time-latitude vorticity maps for every m. Because there is no symmetryζ m (ν, r, λ) =ζ * m (−ν, r, λ), the filtered vorticity mapsζ m (t, r, λ) are in general complex.
LGBS18 do a similar analysis for LCT data, in particular for rotation-averaged maps and filtering within ±30 nHz around the central mode frequencies. We do the entire latitudinal eigenfunction analysis for LCT and RDA, for full time-resolution 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 time-resolution and filtering cases yield consistent results; we thus show only the outcome for the full time-resolution and linewidth filtering. However, LGBS18 take the real part of the complexζ m (t, r, λ). 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 equal-size 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, low-number 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.

Covariance
We calculate, at each m, the temporal covariance of the vorticitỹ ζ 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 ζ =ζ − ζ t is the centered vorticity. The function C m (r, λ) is complex-valued sinceζ m 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.

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 noise-filtered maps, and computeζ nw, m (t, r, λ) =ζ m (t, r, λ)/σ m (r, λ). 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ζ nw, m and the resulting latitude vectors are multiplied by σ m again to undo the weighting. We apply the weighting only to LCT, since the ring-diagram 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.

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.
In the current paper, we find that the LCT and the RDA results are mostly consistent for the near-surface 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 SVD-based 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 latitude-dependent 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  Figure 4 suggests that the real part of the eigenfunctions is more confined to low latitudes for higher values of m. We study the m-dependence 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.
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 m-dependence for those measurements; this dependence is mostly consistent with that of the LCT covariance-based parameters.
The latitude width at Re(C) = 0.5 indeed decreases with m, quasi-linearly 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 m-independent 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 P m (sin λ), 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 P m (sin λ) are normalized such that π/N λ λ (P m (sin λ)) 2 cos λ = 2. 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 near-sectoral 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 near-equatorial latitudes for the imaginary part (for almost all modes). Table 2 shows the coefficients c m for m ≤ ≤ m + 6 for the LCT covariance-based 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 A44, page 9 of 16  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.
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 Brunt-Vä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.

Depth-dependent ring-diagram 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 ring-diagram 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 72-day 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 ring-diagram latitude-depth grid via 2D bicubic splines, which is reasonable because the global-mode inversions do not vary on scales of their original grid; we then average the 72-day chunks over time. The chunk scatter of the rotation rate is used to estimate the uncertainty. We divide the ring-diagram 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 ring-diagram tracking.   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 ring-diagram 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 ring-diagram results are in best agreement with each other. Figure 7 shows the equatorial rotation rate versus depth from global modes and ring-diagram 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 ring-diagram 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.
The most worrisome point is the disagreement between different ring-diagram 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 short-lived flows and even longer-lived A44, page 10 of 16 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 ring-diagram 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 depth-independent underestimation of the ring-diagram 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.

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 maximum-likelihood Lorentzian fit (Anderson et al. 1990) to the power spectra for the longer ring-diagram 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 ring-diagram 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 chunk-based 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 LGBS18 Notes. Previous measurements (with superscript "ref") are also listed for comparison.
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 non-Gaussian) 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.
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 multi-peak 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  (3)). The depth refers to the median of the ring-diagram averaging kernels (orange dots), which corresponds to certain target depths (black dots).
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. 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. Figure 8 also shows the main parameters of the ring-diagram 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 well-localized 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.

Results for the radial eigenfunctions
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 P 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 near-surface 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 non-Gaussian (power cannot be negative). More information about error estimation can be found in Appendix C. 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 P 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.

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, complex-valued 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 large-scale 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 m-independent width at a value of 0 and A44, page 12 of 16   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 ring-diagram 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 ring-diagram longitudes. This indicated systematic effects in the ring-diagram 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 power-law decrease, in particular both with the theoretical Provost et al. (1981) model (exponent m − 1) and an m-independent 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.

Appendix A: Issues of the ring-diagram 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., R 0 K n (r)dr = 1 (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., 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Ũ x, n is the ring measurement in the rotating frame. We now define the local deviation from the tracking rate δΩ(r) = Ω(r) − Ω T . Then we obtaiñ 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 time-independent, since the ring-diagram mode set does not vary much with time. Thus the time-dependent 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Ũ x, n via Eq. (A.6), on which we run the ring-diagram inversion module to retrieve the output velocities. To computeŨ x, n , we assume a depth-independent 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; Christensen-Dalsgaard 2008, 2011. We lose roughly 25% of the original ring-fit 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 mode-fit 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 ring-diagram velocities
We interpolate the ring-diagram velocities separately in time and longitude (see Sect. 2) with different functions, depending on the number of available data points: − ≥ 4 data points : cubic splines −3 data points : quadratic splines −2 data points : linear splines. (B.1) Before we interpolate the ring-diagram velocities to the surface equatorial rotation rate, we apodize these velocities with a raised cosine in angular great-circle distance ρ to the point (λ = 0 • , ϕ = 0 • ), see Sect. 2), as follows: (where there are no more valid pixels). Apodizing the ringdiagram velocities (with different β), or not, gives consistent results.

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 equal-size 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 (rotation-averaged 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 χ 2distributed random variable (stochastic excitation). In analogy to Eq. (12) we generate 1000 realizations of synthetic data for the Fourier transform of the radial vorticity F syn (not the power P), at each m, as F syn, m (ν, r) = A m (r) 4(ν − ν 0, m ) 2 /γ 2 m + 1 N A, m (ν) + B m (r)N B, m (ν, r). (C.1) 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 N A, m (ν) is constant with depth, i.e., the signal is fully correlated in depth, while for the background, we take a random variable N 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 off-diagonal 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.