Issue 
A&A
Volume 649, May 2021



Article Number  A59  
Number of page(s)  16  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202039311  
Published online  11 May 2021 
Helioseismological determination of the subsurface spatial spectrum of solar convection: Demonstration using numerical simulations
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: boening@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, PO Box 129188, Abu Dhabi, UAE
Received:
1
September
2020
Accepted:
18
January
2021
Context. Understanding convection is important in stellar physics, for example, when it is an input in stellar evolution models. Helioseismic estimates of convective flow amplitudes in deeper regions of the solar interior disagree by orders of magnitude among themselves and with simulations.
Aims. We aim to assess the validity of an existing upper limit of solar convective flow amplitudes at a depth of 0.96 solar radii obtained using timedistance helioseismology and several simplifying assumptions.
Methods. We generated synthetic observations for convective flow fields from a magnetohydrodynamic simulation (MURaM) using traveltime sensitivity functions and a noise model. We compared the estimates of the flow amplitude with the actual value of the flow.
Results. For the scales of interest (ℓ < 100), we find that the current procedure for obtaining an upper limit gives the correct order of magnitude of the flow for the given flow fields. We also show that this estimate is not an upper limit in a strict sense because it underestimates the flow amplitude at the largest scales by a factor of about two because the scale dependence of the signaltonoise ratio has to be taken into account. After correcting for this and after taking the dependence of the measurements on direction in Fourier space into account, we show that the obtained estimate is indeed an upper limit.
Conclusions. We conclude that timedistance helioseismology is able to correctly estimate the order of magnitude (or an upper limit) of solar convective flows in the deeper interior when the vertical correlation function of the different flow components is known and the scale dependence of the signaltonoise ratio is taken into account. We suggest that future work should include information from different target depths to better separate the effect of nearsurface flows from those at greater depths. In addition, the measurements are sensitive to all three flow directions, which should be taken into account.
Key words: convection / hydrodynamics / instabilities / Sun: helioseismology / Sun: granulation / Sun: photosphere
© V. G. A. Böning et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://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
Helioseismic inferences of convective motions in the solar interior yield consistent results in nearsurface layers using different techniques (e.g., De Rosa et al. 2000; Braun & Lindsey 2003; Hindman et al. 2004; Gizon & Birch 2004; Langfellner et al. 2014, 2015). One important statistical property of convective motions is their spatial power spectrum, which quantifies the scaledependent amplitude of the velocities. Inferences of the power spectrum of convective motions in deeper regions of the solar interior have led to conclusions that differ by more than an order of magnitude (Hanasoge et al. 2012, hereafter HDS2012; Greer et al. 2015). Understanding these differences and obtaining reliable estimates of the spectrum of convective motions in the Sun is important for several processes in stellar interiors, such as the maintenance of differential rotation, and for understanding the role of convection in the emergence of active regions. This is even more urgent because convective velocities from numerical simulations of solar convection exhibit largescale convective motions that are far too strong, even compared to surface observations (Gizon & Birch 2012; Lord et al. 2014). This last issue is often referred to as the convective conundrum.
To date, the only methods that have been used to constrain the spectrum of convective motions in the deeper solar interior are timedistance helioseismology (Hanasoge et al. 2010, 2012) and ringdiagram analysis (Greer et al. 2015). As helioseismic techniques, both infer solar interior flows from surface observations of stochastically excited sound waves and surfacegravity waves, which propagate through the interior. To do so, timedistance helioseismology (Duvall et al. 1993) uses travel times of waves that travel different distances and hence probe different depths. In ringdiagram analysis (Hill 1988), observations of the local power spectrum of the surface oscillations are used and the Dopplershifted frequency of the seismic waves is measured by way of a background flow field. Gizon & Birch (2005) and Gizon et al. (2010) reviewed these local helioseismic techniques, and Hanasoge et al. (2016) reviewed their application to probing interior convection. In addition, global helioseismology can be used to infer properties of timedependent flows such as convection (e.g., Roth & Stix 2003; Woodard 2016; Mani & Hanasoge 2020). This method relies on global properties of the oscillations such as mode frequencies or amplitudes. Timedependent flows result in mode coupling, which can be measured in global observations and be used to reveal flow properties.
For the case of timedistance helioseismology, Hanasoge et al. (2010, 2012) estimated an upper limit of convective motions based on several simplifying assumptions, notably that the signaltonoise ratio (S/N) of the helioseismic measurements is independent of the horizontal spatial scale and that the signal in traveltime measurements is predominantly due to the flow in one direction at or near a targeted depth.
We here therefore aim to assess the validity of the method suggested by HDS2012 to infer an upper limit on convective flow amplitudes in the solar interior. To do this, we generated a number of realizations of synthetic traveltime maps using background flow fields from existing numerical magnetohydrodynamic (MHD) simulations of quietSun convection (Lord et al. 2014). As we know the input flow, we can assess the performance of the resulting estimate and the underlying assumptions. To generate the synthetic observations, we employed a forward model of the signal and a noise model. In the forward model, the signal in the traveltime maps is computed using the first Born approximation (Böning et al. 2016), which takes singlescattering perturbations to the wave field into account (Gizon & Birch 2002). We obtain realizations of the stochastic noise from a Gaussian noise model in Fourier space (Gizon & Birch 2004; Fournier et al. 2014), which is based on the assumption of a large number of random wave sources homogeneously distributed in space and time. These are welltested methods used in the community (e.g., Böning et al. 2017; Gizon et al. 2020).
2. Methods
We consider the problem of estimating subsurface convective flow amplitudes from simulations in Cartesian geometry. We assume that the statistical properties of the flow do not depend on location within a horizontal layer.
Clearly, the use of Cartesian geometry is a simplification, especially when the interest is on larger horizontal scales. The main effects of this simplification are that (i) the resolution in Fourier space is reduced compared to covering the entire sphere, and (ii) effects of the spherical geometry such as rotation and curvature on the flow are not modeled. These are important shortcomings that should be addressed in future studies. However, given the orderofmagnitude difference between HDS2012 and Greer et al. (2015), our main interest in this work is to verify whether one of the proposed methods (HDS2012) is working at all, and potentially, to give clues on the cause of the differences.
In addition, the essential problem that convective flow amplitudes are larger in simulations than in observations is also present for Cartesian simulations of convective flows (Lord et al. 2014, see the comparison to surface observations in Fig. 5). In Cartesian geometry, more realistic simulations of convection are available than in spherical geometry because radiative transfer can be better included and smaller turbulent scales can be resolved (Kupka & Muthsam 2017).
In this section, we first summarize the methods that we used to model the method of HDS2012 and to create synthetic data. We then summarize the assumptions that were made by HDS2012 to obtain an upper limit on the spectrum of convection to be able to study the validity of these assumptions in the following sections.
2.1. Convective flow fields from MURaM
As input for generating synthetic observations, we used 66 consecutive snapshots of MHD simulations of solar convection. These simulations were obtained by Lord (2014) using the MURaM code (Vögler et al. 2005) with the vertical magnetic flux fixed to zero. The details for these simulations are given in Lord (2014, see the case “zero net flux” in Tables 2.1 and 5.1, as well as Sects. 2 and 5)^{1}. Results for the purely hydrodynamic case of this simulation setup were reported in Lord et al. (2014). After the hydrodynamic simulations had reached a statistically steady state, Lord (2014) initialized the magnetic field with a root mean square magnetic field strength of 0.5 G while keeping the horizontally averaged vertical magnetic flux at zero. The dynamorelated magnetic field amplification resulted in a photospheric root mean square magnetic field strength of 140 G and a photospheric mean unsigned vertical magnetic flux of 37.2 G.
Lord (2014) showed snapshots of the vertical velocity u_{z} of the hydrodynamic case in Figs. 2.2 and 2.3 and horizontally averaged density, temperature, and pressure profiles in Fig. 2.5, which are close to a standard solar model (ChristensenDalsgaard et al. 1996). Snapshots and horizontal averages of background quantities are very similar for the magnetic case (see Lord 2014). We show a snapshot of the horizontal velocity u_{x} in Fig. C.1.
We used 66 consecutive snapshots that were archived by Lord (2014). They were extracted at a cadence of 3.8 h of solar time, and rebinned to 512 × 512 × 384 grid points with a resolution of 384 × 384 × 128 km^{3} in a Cartesian box sized 196 × 196 × 50 Mm^{3} in the (x, y, z) directions, where z points upward and L = 196 Mm is the horizontal extent of the box. The original spatial resolution had twice as many grid points in all directions. The reduction in resolution only changes the smaller scales in the simulation, while in this study, we are interested in the very largest scales, which are relatively unaffected by the reduction in resolution. The simulation extended to a depth of z = −48.5 Mm (r = 0.93 R_{⊙}, or about 16 cumulative pressure scale heights from the surface). This study is targeted at flows around r = 0.96 R_{⊙}, which is about 2.2 local pressure scale heights from the bottom boundary.
The helioseismic data in HDS2012 were averaged over time periods of one day. The data are thus sensitive to temporal averages of convective motions. We here mimicked this temporal averaging by subdividing the simulated cubes into 11 chunks of 6 × 3.8 h = 22.8 h and by averaging the flow field over the six snapshots for each chunk. This is further discussed in Sects. 2.4.3 and 4.
2.2. Forward model in Cartesian geometry
As the simulated flow fields are given in Cartesian geometry and are horizontally periodic, we adopt a Cartesian geometry for the rest of this paper. The Cartesian box was centered at the equator. In order to convert from Cartesian coordinates (x, y, z) into spherical coordinates (r, θ, ϕ) and vice versa, where r is the distance from solar center, θ is the colatitude and ϕ is the longitude, we used the relations x = ϕ R_{⊙}, y = (π/2 − θ)R_{⊙}, z = r − R_{⊙}. To facilitate comparison with HDS2012, we continue to refer to the z coordinate by the corresponding value for r.
In Cartesian geometry, the convolution theorem can be used, see Appendix A, to write the relation between travel times τ and flows u = (u_{x}, u_{y}, u_{z}), which is mediated by the sensitivity function K = (K_{x}, K_{y}, K_{z}) in the Fourier domain,
where we mark the horizontal Fourier transform only by the use of the horizontal Fourier variable or wave vector, K, see Appendix A, and we used ℓ = k R_{⊙} to convert from horizontal wave number k = K into harmonic degree. In Eq. (1), we have summed over all vertical grid points z and horizontal grid points x′, we multiplied by the uniform vertical grid spacing h_{z} and horizontal grid spacing h_{x} = h_{y}, and we used x − x′=(x − x′,y − y′) to define a horizontal convolution. In our computations, we used the discretization in Eq. (1), but for brevity in notation, we write in the remainder of this paper for the vertical integral
We here considered the deepfocus traveltime geometry used by HDS2012 for the target depth of r = 0.96 R_{⊙} and computed traveltime sensitivity kernels K for this measurement using the Born approximation and the method of Böning et al. (2016), see Fig. 1. This method assumes that the flow is constant in time. Sensitivity functions obtained using the Born approximation compare well in different approaches (e.g., Birch et al. 2007; Burston et al. 2015; Böning et al. 2016; Mandal et al. 2017). To compute the sensitivity functions, we used a spatial grid that was four times as wide as the MURaM grid to capture the entire sensitivity function. It was then Fourierprojected on the MURaM grid using the periodicity of the flow field. Furthermore, we used only half the spatial resolution for computing kernels. This is sufficient because the smallest seismic wavelength is still more then ten times larger than the spatial resolution. Appendix B gives further details on the deepfocus traveltime geometry, which was modeled using the kernels. The traveltime geometry was designed by HDS2012 to be predominantly sensitive to u_{x} flows (the eastwest component of the flow) close to a target depth of 0.96 R_{⊙}. Figure 1 shows that the traveltime measurements are sensitive to flows in all spatial directions and especially to nearsurface flows, although there is some degree of focussing for K_{x}.
Fig. 1.
Sensitivity of deepfocus travel times (target depth 0.96 R_{⊙}) to flows. The observational setup is taken from HDS2012 and is described in detail in Appendix B. We computed the sensitivity functions using the first Born approximation (Böning et al. 2016). 
Using these kernels and the flow fields from MURaM, we obtained forwardmodeled traveltime maps using Eq. (2). These travel times are free of seismic noise (ϵ = 0) and are thus referred to as noiseless travel times. The bottom left panel in Fig. 2 shows an example of these maps for one realization of the flow field. The top panels of Fig. 2 show that u_{x} indeed contributes most to the travel times, as intended by HDS2012. However, the u_{y} and u_{z} flow components also contribute significantly, and a large fraction of the u_{x} contribution comes from nearsurface flows (bottom middle panel in Fig. 2), which was unintended. On the other hand, the bottom right panel in the figure shows that the signal from the flows in the deeper layers is somewhat correlated with the total signal, which shows that the deep flows might be detected using these data.
Fig. 2.
Example noiseless travel times τ_{NR} (bottom left panel) corresponding to the right panel in Fig. 3 and contributions from the different flow components (top panels) as well as different depth ranges (remaining panels) obtained using Eq. (2) with ϵ = 0. 
2.3. Noise model and covariance of travel times
In Cartesian geometry and in the quiet Sun, it is reasonable to assume that the correlation in the noise of two traveltime measurements only depends on their separation (e.g., Gizon & Birch 2004),
so that, see Appendix A,
where Λ is the noise covariance matrix and h_{k} is the resolution in Fourier space. Here, we assumed that averages ⟨q⟩=⟨q(t)⟩_{t} for any quantity q are taken over many independent realizations q(t), which are observed for different nonoverlapping periods indexed by t. We assumed that enough data were averaged so that the temporal averages approximate their expectation values.
As a consequence of the horizontal translation invariance, the problem decouples in K space, that is,
We compute the noise covariance matrix for the deepfocus measurement using the method of Gizon & Birch (2004) and Fournier et al. (2014) and draw realizations of the noise as described in Appendix A.
Figure 3 shows an example realization of the synthetic noise. It is qualitatively very similar to the actual data while being a different stochastic realization. This can be seen by comparing the left panel in Fig. 3 to Fig. 5 in the supplementary material of HDS2012. Our synthetic noise has a standard deviation of 10.3 s and the contribution of the traveltime signal to the right panel of Fig. 3 has a standard deviation of 1.9 s. Both compare reasonably well with the data (HDS2012), where the noise has a standard deviation of 11.8 ± 0.3 s and the data 2.5 ± 0.02 s (Birch et al., in prep.). We note that this noise model has been tested successfully against data (e.g., Figs. 5, 6, 8, 10, and 11 in Gizon & Birch 2004; Fig. 5 in Fournier et al. 2014).
Fig. 3.
Example realizations of the noise ϵ(x) from our noise model on the kernel (left) and MURaM grids (middle), and synthetic noisy travel times τ(x) (right). The noise has a standard deviation of 10.3 s and the traveltime signal in the right panel of 1.9 s, which compares well with the data (Hanasoge et al. 2012). Left panel: can be compared to Fig. 5 in the supplementary material of Hanasoge et al. (2012). Figure 2 shows the contribution of the traveltime signal to the right panel in more detail. 
However, we also note that it is not clear why Fig. 3 in HDS2012 shows a very similar traveltime map, but with features aligned in northsouth direction, as opposed to the features aligned in eastwest direction in Fig. 5 in the supplementary information of HDS2012 and in our Fig. 3. Except for the difference in the orientation, the three figures are qualitatively very similar. Because the measurements were averaged (see also Fig. 1), it is expected that traveltimes are more correlated in eastwest direction, as is shown in our noise maps and in Fig. 5 in the supplementary material of HDS2012. On the other hand, the orientation of the noise does not affect the results in our study or in HDS2012. Therefore we do not address this issue further.
2.4. Assumptions for estimating an upper limit
In this subsection, we summarize the essential assumptions leading to the estimate obtained by HDS2012 and introduce their method. We also briefly outline the motivation for these assumptions. After summarizing the assumptions that lead to this estimate, we assess their validity in the following section. We therefore do not discuss their validity in this section. We here explicitly mark the most important assumptions made by HDS2012 to obtain an upper limit.
However, we stress already here that rather than an inequality sign, an approximate equality sign is a more adequate assumption in many cases. In the following, we therefore include alternative formulations of the assumptions as approximations for future reference in parentheses.
2.4.1. Noise correction
A major step is the socalled noise reduction, which aims at estimating the noiseless (or noisereduced) travel times τ_{NR} = τ − ϵ, which is the signal of interest. In order to estimate the amplitude of the signal from the noisy traveltime measurements, HDS2012 used a relation from Gizon & Birch (2004, Sect. 6). Correcting an obvious error in Gizon & Birch (2004), Hanasoge et al. (2012) used the correct version
where σ^{2} is the variance of the data, T is the observation time, is the variance of the noise for an observation time T, and S^{2} is the contribution of the signal to the total variance. This equation is valid in the case of timeindependent flows and omits a small 1/T^{2} contribution to the noise (see Eqs. (14) and (B.6) in Fournier et al. 2014). From Eq. (7), it follows that
and by assuming that the S/N is independent of the spatial scale, HDS2012 obtained
Here, q^{fft}(ℓ) indicates an azimuthally averaged spectrum that was obtained from a horizontal Fourier transform of the quantity q (e.g., q = τ or q = u_{x}) to be compared with the spherical notation of HDS2012. The quantity q^{fft}(ℓ) is defined as a mode amplitude per spherical harmonic degree (same unit as q). As usual, an appropriate summation of q^{fft}(ℓ)^{2} can be used to estimate the variance in real space, see Appendix A for the exact normalization. Starting from Eq. (9), we use a notation like “(A1)” to mark inequalities or equalities that correspond to the most important assumptions from HDS2012 that we assess in this study. These assumptions are summarized in Sect. 2.4.3.
The main aspect of assumption (A1) that we assess in this study is the assumed independence of k or ℓ of the signaltonoise ratio, which is implicit in Eq. (9). HDS2012 fit the S/N in Eq. (7) to the data. We do not study the validity of this fitting procedure here.
In order to estimate the signal from the noisy data, HDS2012 used assumption (A1), which says that the noisereduced traveltime mode amplitude is bounded by (or can be approximated by) the S/N times the observed noisy traveltime mode amplitude τ^{fft}(ℓ). Here, is a spectral decomposition of the signal variance, and similarly, τ^{fft}(ℓ)^{2} is a decomposition of the total variance σ^{2}.
For synthetic data, we know
and we can thus compare this value with the value of S/N used by HDS2012.
2.4.2. Calibration curve
When the amplitude of the signal has been estimated from the data or a noiseless traveltime spectrum is given for synthetic data, the aim is to estimate the amplitude of the convective motions at the depth that generates the signal. HDS2012 proposed to do this using a socalled calibration curve. This calibration curve simply converts the noisereduced traveltime spectrum into a velocity spectrum.
To introduce the derivation of the calibration curve proposed by HDS2012, including the underlying assumptions, we start from Eq. (6). Assuming that the statistical properties of the data are translation invariant, we can only consider correlations ⟨τ^{*}(K) τ(K′)⟩ for K = K′. One then obtains
where, following HDS2012, we assumed that the spectrum of the traveltime signal is bounded below by (or can be estimated by) the spectrum of traveltimes stemming only from the xcomponent of the flow in assumption (A3), or in an intermediate step as a sum of spectra from the three components, ignoring correlations of different flow components, in assumption (A2). These assumptions are motivated by the idea that the measurement geometry is designed to be mostly sensitive to the x direction. u_{x} contributes most to the signal, see the top left and bottom left panel in Fig. 2. However, the measurements are also sensitive to the other flow components, see Fig. 1 and the top middle and top right panels in Fig. 2. We evaluate and discuss the validity of the assumptions made by HDS2012 and the effect on the inferred flow spectrum in detail in Sect. 3.
HDS2012 assumed that the flow near the target depth of the traveltime measurement bounds the spectrum from below (contributes dominantly to the spectrum of travel times). In addition, they assumed that the correlation of the flow at different depths may be bounded (or modeled) by a Gaussian, more specifically that some Gaussian is centered at the target depth z_{T} = (0.96 − 1)R_{⊙} (z = 0 is at r = R_{⊙}) with a width of σ_{z′} = 16.77 Mm = 1.8H_{P}, where 1.8 pressure scale heights H_{P} is the mixing length. It is required here that the Gaussian falls off faster than the correlation function, such that
and consequently, if ,
where
We did not confirm whether the assumption of a Gaussian shape of the correlation function is an appropriate choice. This is left to future studies.
The final assumption is
where I_{ℓ} = {K such that K−h_{k}/2 ≤ ℓ/r < K + h_{k}/2} selects all horizontal wave vectors of similar length, N_{ℓ} is the number of elements in I_{ℓ},
and ℓ = ℓ_{j} takes the values ℓ_{j} = j h_{k}r with integer j in order to be compatible with the discrete Fourier grid. This introduces the socalled calibration curve, . By making assumption (A5) in Eq. (17), we assume that the azimuthal dependence of the flow and the kernel spectra are not anticorrelated, such that the azimuthal average of the kernels and flow spectra can be taken independently, without increasing the result. This was assumed but not explicitly mentioned by Hanasoge et al. (2012).
By taking the azimuthal average of the traveltime power spectrum in Eq. (11) and by using the above assumptions, one obtains the final inequality on which the estimate of the upper limit of the flow amplitude by HDS2012 is based. At the same time, we convert from Cartesian into spherical geometry in the Fourier domain using ℓ = kr and h_{ℓ} = h_{k}r. Here and in the definition of I_{ℓ} above, we used r ≈ const = R_{⊙} for the entire domain, which is equivalent to the assumption of Cartesian geometry. In addition, we employed the convention for converting spectra from Fourier space into spherical harmonic space outlined in Appendix A, which is inspired by Birch et al. (in prep.),
In order to facilitate the numerical computation of Eq. (16), HDS2012 instead chose a smaller width , where D = 9.64, see also Birch et al. (in prep.). As HDS2012 derived 𝒦 or from numerical simulations of wave propagation, shrinking the correlation width by the factor D permitted them to deduct a larger number of independent measurements from one simulation, and thus to reduce the numerical cost. However, the quantity D at the same time acts as an additional free parameter. We show below that its introduction has important effects on the estimated spectrum.
After obtaining 𝒦(K, σ_{z}) and the alternative calibration curve
HDS2012 then approximated by
The integration around the target depth with a Gaussian weighting of width is equivalent to assuming that the flows are vertically correlated with each other for about a scale height. This is reasonable for the simulated data (Lord et al. 2014), although a scale dependence of this width may be a more adequate assumption. At the same time, the value of the vertical correlation length of the flows is not known for the Sun.
2.4.3. Summary of assumptions
As a result of the above assumptions, we obtain, see Eqs. (9), (21), and (23),
where the equivalent of the righthand side was used by HDS2012 to estimate the flow near the target depth. We write for the estimate of u_{x}. We now briefly summarize the most important assumptions that were explicitly marked above and that were proposed by HDS2012. We evaluate them in the following section.
(A1) The variance in the traveltime maps stemming from the signal and due to the flow is bounded above by (or can be estimated by) obtaining the S/N from a fit to the data (see Hanasoge et al. 2012). The S/N is independent of spatial scale, see Eq. (9).
(A2) The noisereduced spectrum of travel times is bounded from below by (or can be estimated by) taking only correlations of a flow component with itself into account, see Eq. (12).
(A3) The spectrum of travel times is further bounded from below by (or can be estimated by) taking only correlations of the xcomponent of the flow with itself into account, see Eq. (13). This assumption stems from the observation that the traveltime geometry is chosen such that the major contribution is in fact from this component.
(A4) The spectrum of travel times is further bounded from below by (or can be estimated by) the contribution of the dominant depth range, which is assumed to be a pressure scale height around the target depth of the deepfocus measurement, see Eq. (15).
(A5) The azimuthal dependence of u_{x} and K_{x} is independent, that is, it does not matter whether these quantities are first multiplied and are then averaged over the azimuth of k, or whether the average is taken first and is then multiplied, see Eq. (17).
(A6) To obtain the calibration curve, similar results are obtained regardless of whether the kernels are integrated over about a pressure scale height, or over the smaller width and then the result is multiplied by D, that is, , see Eq. (23).
Furthermore, a few more assumptions were (at least implicitly) made by HDS2012. We state these assumptions here for completeness; these assumptions are not be studied here but should be addressed in future work:
(B1) In fitting the S/N, HDS2012 followed the assumption of Gizon & Birch (2004) that the noise decreases with . A potential 1/T term (see Eqs. (13) and (B.6) in Fournier et al. 2014) is thus neglected.
(B2) The effect of a time dependence of the convective flows, either on the fit of the S/N through a time dependence of the signal, or on the travel times beyond a Bornapproximation model for constant flows (Gizon & Birch 2002), can be neglected. Because convective flows are timedependent, timedependent perturbation theory should be used in an ideal scenario.
(B3) The convolution theorem can be used. This is straightforward for Cartesian geometry, see Eq. (2), but might be inaccurate for spherical geometry.
(B4) The vertical correlation of the flow components can be approximated (or bounded) by a Gaussian, see Eq. (14). The possibility of an anticorrelation of divergent flow components (due to mass conservation) is thus not taken into account, but should be modeled in future work.
3. Results
3.1. Scaledependent noise correction is important
For MURaMtype flows, we find that the S/N significantly depends on the spatial scale and is not independent of scale, as assumed by HDS2012, see the top panel in Fig. 4 and assumption (A1). The true S/N for the Sun may be different to what we report here if the subsurface flows in the Sun are substantially different from MURaM flows. We strongly recommend that the noise correction be made as a function of harmonic degree (or wave number), possibly even as a function of azimuth (i.e., m or k_{x}/K).
Fig. 4.
Inverse of the S/N (top) and its effect on the noise reduction (middle) as well as the estimate of the spectrum of horizontal flows (bottom, as a velocity amplitude per mode) for different approaches. The S/N is clearly scale dependent, in contradiction to assumption (A1). We use ℓ = k R_{⊙} to convert from wave number into harmonic degree. We do not show results for ℓ = 0 because this corresponds to an arbitrary change in reference frame. 
The middle and bottom panels in Fig. 4 show the effect of this scaledependent noise correction on the spectrum of travel times and the estimated flow amplitudes. For MURaM flows, a constant S/N for all scales clearly leads to an underestimation of the flow amplitudes at the largest scales.
3.2. Travel times are sensitive to nearsurface flows and to all flow directions
The deepfocus traveltime geometry of HDS2012 is designed to be sensitive mostly to u_{x} flows (the eastwest component of the flow) around the target depth of 0.96 R_{⊙}. Figure 1 shows that the travel times are in principle sensitive to all flow components and to a broad range of depths, including flows near the surface. However, there is some focus on the target depth for the sensitivity to u_{x} flows, see Fig. 1 (bottom left panel).
This is also reflected in the contribution of the different flow components and depth layers to the traveltime signal, which we show in Fig. 2. Although u_{x} (top left) contributes most to the signal (bottom left), u_{y} and u_{z} also contribute significantly (top center and right). The contribution from the surface layers (bottom center) is also large. However, the signal from the deeper u_{x} flows (bottom right) is visible in the signal, which means that u_{x} flows around the target depth might be detected using this method.
3.3. Calibration curve can be reproduced, but depends on assumed flow properties
Second, we find that the computation of the calibration curve, which was done by HDS2012 using simulations of wave propagation, can be reproduced, see the blue line in Fig. 5. The agreement is quite good when compared to Fig. 2 of the supplementary material of HDS2012, even when we computed it on the coarsergrained wave number grid of the MURaM simulation, see the green dashed line in Fig. C.2.
Fig. 5.
Calibration curve of Hanasoge et al. (2012, Fig. 2 in the supplementary material) can be reproduced (blue solid curve). When a more realistic width σ_{z′} = 16.8 Mm = 1.8 H_{P} of flow features is used to compute the calibration curve, a change is seen (dashed orange curve), which would lead to an increase in the estimated upper limit, see assumption (A6). The crosses and pluses indicate the grid points in harmonic degree, which were obtained by converting horizontal Fourier modes from Cartesian to spherical geometry (see Appendix A.4). The kernels and thus the calibration curves shown here were computed on a spatial grid with a four times larger spatial extent than the MURaM grid, therefore the resolution in Fourier space, , is four times finer than in other plots in this paper. See Fig. C.2 for a more detailed comparison. 
However, HDS2012 computed the calibration curve assuming a width of convective features of σ_{z} = 1.74 Mm and later rescaled the result by a factor of D = 9.64 to mimic an effective width of , which is similar to the mixing length or 1.8 pressure scale heights at 0.96R_{⊙} (see Sect. 2.4.2 and Birch et al., in prep.). We compare this to directly assuming a larger width of convective features of , see in Fig. 5. This alternative calibration curve is lower. Using thus results in a correspondingly higher estimated upper limit, see Eq. (24).
In other words, the use of the width and rescaling the calibration curve by using DC_{ℓ} as done by HDS2012 effectively results in a lower estimate for the upper limit. This difference arises because the kernel is not constant in the area around the target depth, see Fig. C.4. This increases the value of DC_{ℓ} compared to . Consequently, assumption (A6) is thus not met and the resulting inequality is in the wrong sense. As the value of σ_{z} is an assumption on the vertical correlation length of the flow, an unknown property of the Sun enters the method here. In addition, the vertical correlation function is approximated by a Gaussian, although its shape is not known, see the discussion in Sect. 3.7.
The introduction of the rescaling with D acts similar to a free parameter. The overall effect of this is further discussed in Sect. 3.6.
3.4. Concept of a calibration curve is generally valid
Furthermore, we confirmed the assumptions (A2)–(A6) that were made to derive the calibration curve, see Fig. 6. While some of the inequalities were met, some were violated, although not substantially. We conclude from this that the concept of a calibration curve is generally valid when the resulting estimate is rather understood as a rough estimate of the amplitude of the flow and not a strict upper limit. We recall that the calibration curve depends on the assumed vertical correlation function, which is not known for the Sun, see the discussion in Sect. 3.7.
Fig. 6.
Test of the validity of the assumptions that enter the derivation of the calibration curve, shown as a noisereduced signal in travel times (left panel) and flow amplitude (center and right panels) as a function of scale. If all assumptions (A2)–(A5) were satisfied, the different curves in each panel would not cross each other, but be ordered from higher to lower values in the same sense as they are ordered from top to bottom in the legends (e.g., Eq. (12) ≥ Eq. (13) ≥ Eq. (15) ≥ Eq. (17) for the left panel). 
3.4.1. Azimuthal averages have to be treated more carefully
We find that when the effect of offdiagonal flowcorrelation is neglected, assumption (A2) is not justified because the offdiagonal elements make a significant negative contribution, see Fig. C.6. The reason for this is as follows. When the dependence of and on the direction of k is considered, we observe that the two quantities are anticorrelated, see the left panel in Fig. C.5. This is the case for almost all k and results in a negative total effect of the offdiagonal terms on the traveltime spectrum. It might be possible to remove this anticorrelation in future studies by using a pure pointtopoint traveltime geometry, where cross correlations are not averaged over arcs as done by HDS2012, or by considering a decomposition in divergent and vortical flows, rather than x and y components, see also the discussion in Sect. 4.
We find that the diagonal elements have a positive contribution (everywhere but for a small region for d = z, see Fig. C.6). As a consequence, neglecting d = d′=y and d = d′=z, as in assumption (A3), is justified to obtain an upper limit. However, we find that all flow components contribute significantly to the traveltime spectrum, see Fig. 6 (compare the blue, orange, and green curves in the left panel). This is consistent with all flow components contributing to the traveltime signal (see Sect. 3.2). To obtain an accurate estimate of the flow spectrum, the effect of all flow components thus has to be taken into account.
Similar to the offdiagonal case, neglecting the dependence of and on the direction of K introduces an error, see the right panel in Fig. C.5. The magnitude of the error depends on k and r/R_{⊙}, but it is generally about 20% (as the effect on the spectrum, not on the mode amplitude). The resulting effect on the traveltime spectrum is of the wrong sign and artificially increases it, see Fig. 6. Strictly speaking, assumption (A5) is thus not justified.
3.4.2. Major contribution from nearsurface flows, small contribution from target depth
We find that the major contribution to the traveltime spectrum comes from near the surface. Only about 10–30%, rarely up to 40%, of the traveltime spectrum is produced below 0.99 R_{⊙} or in a Gaussian envelope with width around the target depth, see Fig. C.3. Assumption (A4) is thus well satisfied as an inequality, but is far from being a good approximation. This result is consistent with the sensitivity kernels peaking near the surface (Fig. 1) and with the traveltime signal having a dominant contribution from the surface (Sect. 3.2).
3.5. Rather a rough estimate than a strict upper limit
Finally, we compared the resulting estimates of the flow amplitude with the actual estimates in Fig. 7. Taking all previous results into account, we find that it is more adequate to consider the estimates to be rough estimates of the magnitude of the flow than a strict upper limit. Equivalently, we find it to be more accurate to write an approximate equality sign instead of an inequality in most cases for assumptions (A1)–(A6).
Fig. 7.
Comparison of estimates of the spectrum of convective flows obtained from the synthetic data based on flows from a MURaM simulation. The method of HDS2012 applied to this synthetic data (solid blue line) can be improved by using a scaledependent signaltonoise correction (dashed orange line, this work). Both estimates give the correct order of magnitude of the flow (green line), although the estimate from the method of HDS2012 is not a strict upper limit. When the underlying assumptions are improved (solid orange line), our best estimate results in a larger upper limit. 
Figure 7 compares the estimates
with the amplitude of the timeaveraged flow, , which was averaged over six snapshots that were taken within the period of one day (see Sect. 2.1). We averaged all estimates and actual flow amplitudes over the 11 available realizations. These averages were taken over the corresponding spectra before the square root was taken to obtain the flow amplitude.
3.6. Why the method gives the correct order of magnitude as a result
Considering the results shown in Fig. 7, we conclude that for MURaMlike flows, the method of HDS2012 just as applied in that paper (blue curve) gives the correct order of magnitude of the flow within a factor of about two for the larger scales with ℓ ≲ 90. For scales with ℓ < 60, the estimate is not an upper limit because it underestimates the true flow by a factor of about two. For the smaller scales with ℓ > 100, the estimate is too large by nearly an order of magnitude, although it is a true upper limit.
When the noise correction is applied correctly and the free parameter D is eliminated (solid orange curve), the agreement is better. In nearly the entire range (ℓ < 110), the improved estimate constitutes an upper limit and is consistent with the actual flow spectrum within a factor of two.
Given that for both cases some of the inequalities are not met, see Fig. 6 and Sect. 3.4.1, the question remains why the method works, at least for estimating the order of magnitude. To understand this, we first note that considering only the effect of d = x flow correlations roughly gives the right amplitude of the traveltime spectrum. This is because the cumulative effects of the other flow components are partly positive and partly negative, and roughly cancel out, see Fig. 6 (blue versus green curve in left panel) and Fig. C.6 (cumulative effect of assumptions A2 and A3).
When we only take the flow near the target depth into account, the traveltime spectrum is greatly reduced, see Fig. C.3. The flows around the target depth only contribute about 10–40% to the squared signal. See also the (A4) curve in Fig. 6 and the corresponding estimate in the right panel in Fig. 6. This effect is almost completely compensated for by the amplification from the change from the calibration curve to DC_{ℓ}, see middle panel in Fig. 6. By this change, a free parameter D is effectively introduced by HDS2012, which compensates for the inaccuracies that occurred previously. We can only speculate about the reason why the compensation using D is apparently of the correct magnitude to yield a tight upper limit, as can be seen in the middle panel of Fig. 6 and for the dashed line in Fig. 7. It may be due to chance or arise because Hanasoge et al. (2010) and HDS2012 calibrated the inferred spectrum of convection with the flows from simulations, and hence intuitively chose a value for D that gives the correct magnitude of the estimated spectrum for the flow field used in their study.
For smaller scales (ℓ > 110), the estimates become worse, as can be seen toward the right end of Figs. 7 and 6. This is because for increasingly smaller horizontal scales, the sensitivity of the travel times to flows near the target depth decreases and the sensitivity to nearsurface flows increases.
3.7. Importance of the vertical correlation structure and the surface contribution
Because the signal receives a major contribution from nearsurface flows (Sects. 3.2 and 3.4.2, Fig. C.3), a correct conversion of a traveltime spectrum to a flow spectrum always has to take the relation of the surface flows to the deep flows into account, see especially Fig. C.6. It would be possible to measure this correlation function for MURaM flows, but the true vertical correlation function for the Sun is not known. We therefore find as a major conclusion that the vertical correlation structure of the flow is needed as a prior knowledge for this method to be usable, or equivalently, that its result crucially depends on the assumed correlation structure.
This point concerns not only the correlation of the surface with the deeper regions, but also the relative amplitude of the deep flows compared to those closer to the surface. In other words, without knowing the effect of the nearsurface layers to the traveltime spectrum, it is impossible to infer the spectrum at greater depths. Inferences of the spectrum at depth from observations can thus only be trustworthy when the spectrum is inferred at several depths at once, for example, using information from different target depths and an inversion procedure. Using different target depths would in principle also include information on the vertical correlation structure.
Alternatively, it may be possible to improve the observational setup so that the sensitivity kernel is better focused at the target depth (see the bottom left panel in Fig. 1 for the current lack of focus and Sect. 3.2). If the measurement were only sensitive to flows close to the target depth, the problem discussed above would be avoided. Reducing the sensitivity to nearsurface flows would be a first step toward such an improvement.
4. Conclusions
We have performed an independent evaluation of the method suggested by HDS2012 to estimate an upper limit of solar convective flow velocities from helioseismic measurements at a target depth of r = 0.96 R_{⊙}. We did this using synthetic helioseismic data for which the true flow is given by the output of a larger scale MURaM MHD convection simulation of the quiet Sun (Lord 2014) similar to Lord et al. (2014). Furthermore, we used consistent normalizations of the relevant spherical harmonic and Fourier transforms as proposed by Birch et al. (in prep.).
For the case of MURaM flows, we find that the estimates obtained using this method give the correct order of magnitude of the actual flow amplitudes within a factor of two. However, at the largest scales, the actual flow amplitude is underestimated and the upper limit is too low by a factor of about two. This is due to the scale dependence of the S/N, which was not considered in the original study (HDS2012).
Because several simplifying assumptions enter the method we studied, the question arises why it roughly works. One way to look at this is by considering individual effects of every assumption made. When the correction for the S/N is made, we find that when only the x component of the flow is considered, this is a reasonable approximate model for the traveltime spectrum, although the other components contribute significantly as well. Even though the flows around the target depth then only contribute 10–40% to the traveltime spectrum, the effect of this on the estimated upper limit is nearly compensated for by a rescaling of the calibration curve introduced by HDS2012. As a result, the estimated upper limit is just above the actual flow spectrum. If the rescaling of the calibration curve were not done, the upper limit would be higher by a factor of about 1.5. The dilemma is that the answer to how the calibration curve is correctly computed depends on the vertical correlation length of the actual flows in the Sun, the shape of the vertical correlation function, and the ratio of the amplitudes at depth to the surface, none of which are known. Another way to look at this is as follows. For MURaM flows, even though some inequalities that each corresponds to an assumption are not met, others are met with substantial leeway, so that the final estimate still is an upper bound, except at large scales because of the scale dependence of the signaltonoise correction.
Several lessons can be learned from our study that should be taken into account in future inferences of the flow spectrum at depth. The first is the mentioned scale dependence of the signaltonoise ratio.
Second, the uncertainty of how the calibration curve is to be computed is substantial. This is because the correct calibration curve depends on the properties of the actual flows in the Sun, especially the vertical correlation function of the horizontal flows, which is not known. In addition, the traveltime measurements are more sensitive to nearsurface flows than to flows near the target depth. As the ratio between the nearsurface spectrum and the spectrum at depth is not known, the effect of the nearsurface flows on the measured traveltime spectrum is hard to infer from observations of only one target depth. Possible solutions to this may be performing inversions using several target depths to infer the vertical correlation structure to improve the goodness of the focus at depth (e.g., using helioseismic holography, Lindsey & Braun 2000), or to fit several different models for power spectra and vertical correlation lengths to the data.
Third, the dependence of the kernel and the flow on the azimuthal angle in the Fourier domain has to be taken into account. As a result of this dependence, correlations of u_{x} with u_{y} have nonnegligible effects on the traveltime spectrum. A solution may be decomposing the flow into divergent and vortical components, rather than x and y, or using pointtopoint or pointtoannulus geometry rather than an arctoarc traveltime geometry. Alternatively, because the problem decouples in K space, an inversion may solve for every k individually.
This study may be extended in future work. First, the effect of the time dependence of the flow has not been studied, neither on the S/N due to the decrease of the signal with time, nor by forwardmodeling the effect of timedependent flows on the travel times. Second, the potential relevance of a 1/T^{2} term in the noise model, which may arise in addition to the usual 1/T dependence of the variance (see Eqs. (13) and (B.6) in Fournier et al. 2014), may be studied analytically. Third, this study should be extended to spherical geometry and to include different types of flow fields, including simulations with different spectral shapes for the flow field (e.g., Cossette & Rast 2016; Featherstone & Hindman 2016) and different flow topologies (e.g., Spruit 1997; Brandenburg 2016; Bekki et al. 2017; Hotta et al. 2019; Anders et al. 2019).
On the one hand, the results presented here give evidence that the method of HDS2012 correctly estimates the order of magnitude of the flow under the assumption that solar flows were similar to MURaM (if consistent normalizations of spherical harmonics are used, see also Birch et al., in prep.). However, the observed spectrum then leads to the conclusion that the convective conundrum persists at depth. If the vertical correlation function of flow components in the Sun is very different to MURaM, for example, with a very different correlation length, this conclusion may well be different. Determining the correlation function of solar convection may be simpler if a parameterized version of this function were available, for example, obtained by comparing different simulations.
On the other hand, we cannot resolve the inconsistency of these results with those obtained by Greer et al. (2015) using ringdiagram analysis at this stage. A deeper analysis of that method is required to resolve this issue, see also Nagashima et al. (2020). We are confident that by a careful and detailed evaluation of all aspects of both methods this discrepancy can be understood and resolved. In addition, global helioseismic methods might soon be used to infer the spectrum of solar convection at depth (e.g., Roth & Stix 2003; Woodard 2016; Mani & Hanasoge 2020). Helioseismology of convection is thus becoming a promising topic to be studied in the near future.
The Ph.D. Thesis of Lord (2014) “Deep Convection, Magnetism and Solar Supergranulation” may be downloaded from https://ui.adsabs.harvard.edu/abs/2014PhDT.......241L/abstract.
Acknowledgments
We thank Matthias Rempel for providing the MURaM simulations. VB thanks the referee for valuable feedback and Robert Cameron, Markus Roth, and Shravan Hanasoge for valuable discussions. This work was supported in part by the Max Planck Society through a grant on PLATO Science. The computational resources were provided by the German Data Center for SDO through a grant from the German Aerospace Center (DLR). We acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218. This work used NumPy (Oliphant 2006; van der Walt et al. 2011), matplotlib (Hunter 2007), and SciPy (Virtanen et al. 2020).
References
 Anders, E. H., Manduca, C. M., Brown, B. P., Oishi, J. S., & Vasil, G. M. 2019, ApJ, 872, 138 [Google Scholar]
 Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851 [Google Scholar]
 Birch, A. C., Gizon, L., Hindman, B. W., & Haber, D. A. 2007, ApJ, 662, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Böning, V. G. A., Roth, M., Zima, W., Birch, A. C., & Gizon, L. 2016, ApJ, 824, 49 [Google Scholar]
 Böning, V. G. A., Roth, M., Jackiewicz, J., & Kholikov, S. 2017, ApJ, 845, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A. 2016, ApJ, 832, 6 [Google Scholar]
 Braun, D. C., & Lindsey, C. 2003, in GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. SawayaLacoste, ESA Spec. Publ., 517, 15 [Google Scholar]
 Burston, R., Gizon, L., & Birch, A. C. 2015, Space Sci. Rev., 196, 201 [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [Google Scholar]
 Cossette, J.F., & Rast, M. P. 2016, ApJ, 829, L17 [NASA ADS] [CrossRef] [Google Scholar]
 De Rosa, M., Duvall, T. L., Jr, & Toomre, J. 2000, Sol. Phys., 192, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, T. L., Jeffferies, S. M., Harvey, J. W., & Pomerantz, M. A. 1993, Nature, 362, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Featherstone, N. A., & Hindman, B. W. 2016, ApJ, 830, L15 [Google Scholar]
 Fournier, D., Gizon, L., Hohage, T., & Birch, A. C. 2014, A&A, 567, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., & Birch, A. C. 2002, ApJ, 571, 966 [Google Scholar]
 Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [Google Scholar]
 Gizon, L., & Birch, A. C. 2005, Liv. Rev. Sol. Phys., 2, 6 [Google Scholar]
 Gizon, L., & Birch, A. C. 2012, Proc. Nat. Acad. Sci., 109, 11896 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Birch, A. C., & Spruit, H. C. 2010, ARA&A, 48, 289 [Google Scholar]
 Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020, Science, 368, 1469 [Google Scholar]
 Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17 [Google Scholar]
 Hanasoge, S. M., Duvall, Jr., T. L., & DeRosa, M. L. 2010, ApJ, 712, L98 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Nat. Acad. Sci., 109, 11928 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, F. 1988, ApJ, 333, 996 [NASA ADS] [CrossRef] [Google Scholar]
 Hindman, B. W., Gizon, L., Duvall, J. T. L., Haber, D. A., & Toomre, J. 2004, ApJ, 613, 1253 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Iijima, H., & Kusano, K. 2019, Sci. Adv., 5, eaau2307 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comput. Astrophys., 3, 1 [Google Scholar]
 Langfellner, J., Gizon, L., & Birch, A. C. 2014, A&A, 570, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindsey, C., & Braun, D. C. 2000, Sol. Phys., 192, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Lord, J. W. 2014, Ph.D. Thesis, University of Colorado at Boulder [Google Scholar]
 Lord, J. W., Cameron, R. H., Rast, M. P., Rempel, M., & Roudier, T. 2014, ApJ, 793, 24 [Google Scholar]
 Mandal, K., Bhattacharya, J., Halder, S., & Hanasoge, S. M. 2017, ApJ, 842, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Mani, P., & Hanasoge, S. 2020, ApJ, 901, 139 [Google Scholar]
 Nagashima, K., Birch, A. C., Schou, J., Hindman, B. W., & Gizon, L. 2020, A&A, 633, A109 [EDP Sciences] [Google Scholar]
 Oliphant, T. E. 2006, A guide to NumPy (USA: Trelgol Publishing), 1 [Google Scholar]
 Roth, M., & Stix, M. 2003, A&A, 405, 779 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spruit, H. C. 1997, Mem. Soc. Astron. It., 68, 397 [NASA ADS] [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Woodard, M. F. 2016, MNRAS, 460, 3292 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Detailed formulas and derivations
A.1. Fourier transform convention
For the forward and inverse discrete Fourier transforms of 2D spatial variables, we used the convention (see also Gizon & Birch 2004)
where K = (k_{x}, k_{y}) is the Fourier variable (wave vector), x = (x, y) is the horizontal spatial variable, h_{x} = L/N_{x} is the sampling distance in x, h_{k} = 2π/L is the resolution in Fourier space, L the spatial extent of the data in one direction, and N_{x} is the number of points in dimension. We used the same symbol for a quantity and its Fourier transform and indicate the Fourier transform by the use of the Fourier variable.
Then, the Fourier transform of the convolution
is
Parseval’s theorem becomes, with Y(x) = X(−x)
A.2. Travel times
For the noiseless forwardmodeled travel times, this means
A.3. Covariance and noise
The noise covariance matrix is defined as
In order to generate realizations of the noise, we set
or equivalently,
with
where 𝒩(K) are complex Gaussians with independent real and imaginary parts and unit variance that are independent for different kK, except for the condition 𝒩(−K) = 𝒩(K)^{*}.
Then, we can easily verify that Eq. (A.8) is verified,
We use σ^{2} = ⟨ϵ(x)^{2}⟩=Λ(x − x′ = 0).
A.3.1. Relation to Eq. (33) in Gizon & Birch (2004)
As a consequence of the above definitions, we have
which differs by a factor of from Eq. (33) in Gizon & Birch (2004).
A.3.2. Continuous setting
In a continuous setting (Gizon & Birch 2002),
To transform back to a discrete setting, we can substitute
and we also obtain C = 1/h_{k}.
A.4. Power spectra and equivalent mode amplitudes in spherical harmonics
Following Birch et al. (in prep.), we converted the power spectrum of any quantity from Fourier space into an equivalent power spectrum in spherical harmonic space with the same variance at a single location or per area. As our Fourier transform definition is differently normalized than in Birch et al. (in prep.), we here summarize all formulas appropriate for the convention chosen here.
For any quantity q (e.g., q = τ or q = u_{x}), we write
where h_{k}r is the equivalent resolution in harmonic degree, I = {j h_{k}j = 0, …, N_{x} − 1}, the sets I_{ℓ} are defined in Sect. 2.4.2, and we use ℓ = kr to convert wave number into spherical harmonic degree. We then define the amplitude per mode by
so that
where q(ℓ)^{2} is an interpolated power per spherical harmonic mode,
where is an interpolated version of . As we work in Cartesian geometry, this is equivalent to assuming that the radial coordinate r is approximately constant throughout the domain. We thus assume r = const = R_{⊙}. For u_{x} at r = 0.96 R_{⊙}, the difference due to that compared to assuming r = const = 0.96 R_{⊙} is only .
The corresponding energy spectrum E_{q} as defined by Birch et al. (in prep.) satisfies
In summary, the quantity E_{q} is directly comparable to E_{ϕ} from Birch et al. (in prep.) when we here choose q = u_{x}.
Appendix B: Details on traveltime geometry
To compute the kernels, we used the data analysis filters and the geometrical setup of HDS2012, which are summarized as follows. The data were filtered using a highpass filter that includes a raised cosine,
and a wide phasespeed filter that targets r = 0.96 R_{⊙},
where ν_{0} = 20 μHz, HWHM = 6.7μHz, ℓ_{max} = 383, and the filter is applied by multiplying the power spectrum by f_{1}f_{2}, or equivalently, by multiplying the Fouriertransformed data by .
From the filtered data, HDS2012 measured traveltime differences in an arctoarc way. This was done by obtaining crosscorrelation functions for pointtopoint measurements, and then averaging the crosscorrelation functions over arcs before fitting the travel time. Here, they used only opposite points on the arcs. Afterward, HDS2012 averaged the travel times with a simple arithmetic mean over the different arcs. The arcs are on concentric circles, have a width of 90 degrees, and are in the east and west quadrants. Table B.1 gives the distances of the arcs to the central point. In the case of Δ_{1} = Δ_{2}, the setup is symmetric about the central point and only one pair of arcs exists. Otherwise, both pairs of arcs, one with Δ_{1} to the east and Δ_{2} to the west of the central point, and the other vice versa, were used, see Figure 2 in HDS2012.
Details of the deepfocus traveltime geometry.
HDS2012 projected the data using Postel’s projection, centered on the central point between the arcs. To keep the computational burden for the kernels and the noise covariance manageable, we did not use Postel’s projection, but approximated the number of points on each arc by , where we used the floor function. The spatial resolution of the input data is Δx = 0.46875 ° × π/180°, which corresponds to 5.69 Mm at disc center.
Appendix C: More details for the validity of the assumptions
In this appendix, we give more detailed background information for our results. Figure C.1 shows an example snapshot for u_{x} taken at the middle of the time period relevant for Fig. 2.
Fig. C.1.
Example realization of the horizontal flow component u_{x} at the surface (z = 0). The displayed snapshot is taken in the middle of the time period relevant for the travel times shown in Fig. 2 and is saturated at 33% of its maximum absolute value. Example realizations for u_{y} are statistically similar but transposed, and examples for u_{z} can be found in Lord (2014, Figs. 2.2 and 2.3). 
Figure C.2 shows a more detailed comparison of the calibration curves C_{l} and , computed for the kernel grid and after interpolating the kernels on the MURaM grid. Figure C.3 shows the contribution of the region below 0.99 R_{⊙} to the traveltime spectrum.
Fig. C.2.
Comparison of calibration curves in a more detailed way than in Fig. 5. When the kernels are interpolated onto the MURaM grid, the results do not substantially change (dotted curves), except at ℓ ≈ 20. The crosses and plusses indicate the grid points in harmonic degree, which were obtained by converting horizontal Fourier modes from Cartesian into spherical geometry (see Appendix A.4). The kernels were computed on a spatial grid with a four times larger spatial extent than the MURaM grid, therefore the resolution in Fourier space, , is four times finer. 
Fig. C.3.
Contribution to the noiseless traveltime spectrum from below 0.99R_{⊙} (solid blue line) and from a Gaussianweighted region around the target depth 0.96R_{⊙} (dashed blue) and their forward models based on the actual flow and two versions of the calibration curve (solid orange and dashdotted green). 
Finally, we show a few more detailed plots to verify the assumptions for selected spatial scales k R_{⊙} and depths z. Plots for other spatial scales and depths are qualitatively similar. As a shorthand for Figs. C.4–C.6, we write for any quantity Q(K)
and we use ℜ[] and ℑ[] for the real and imaginary parts of a complex number.
Fig. C.4.
Azimuthal averages as defined in Eqs. (C.1) and (C.3) saturated at 10% of the maximum value of the left panel. The units are s^{4} cm^{−4} (i.e., cgs units). The quantities are zero for d ≠ d′, see also Fig. C.5. 
Fig. C.5.
Exemplary dependence of kernels and flow correlations on the direction of the horizontal Fourier vector. KLeft panel: , right panel: . For the plot, each quantity was divided by its maximum value in the given range. Error bars of C^{u} were obtained from the standard deviation of the 11 available simulation snapshots and are similar for C^{K}C^{u}. At other depths and similar wave numbers, the plots are qualitatively similar. Right panel: we also show that taking the mean before or after multiplying and changes the result (thick solid black line). 
Fig. C.6.
Azimuthal averages as defined in Eqs. (C.1) and (C.2) saturated at 10% of the maximum value of the top left panel. The units are s^{2} cm^{2} (i.e., cgs units). 
All Tables
All Figures
Fig. 1.
Sensitivity of deepfocus travel times (target depth 0.96 R_{⊙}) to flows. The observational setup is taken from HDS2012 and is described in detail in Appendix B. We computed the sensitivity functions using the first Born approximation (Böning et al. 2016). 

In the text 
Fig. 2.
Example noiseless travel times τ_{NR} (bottom left panel) corresponding to the right panel in Fig. 3 and contributions from the different flow components (top panels) as well as different depth ranges (remaining panels) obtained using Eq. (2) with ϵ = 0. 

In the text 
Fig. 3.
Example realizations of the noise ϵ(x) from our noise model on the kernel (left) and MURaM grids (middle), and synthetic noisy travel times τ(x) (right). The noise has a standard deviation of 10.3 s and the traveltime signal in the right panel of 1.9 s, which compares well with the data (Hanasoge et al. 2012). Left panel: can be compared to Fig. 5 in the supplementary material of Hanasoge et al. (2012). Figure 2 shows the contribution of the traveltime signal to the right panel in more detail. 

In the text 
Fig. 4.
Inverse of the S/N (top) and its effect on the noise reduction (middle) as well as the estimate of the spectrum of horizontal flows (bottom, as a velocity amplitude per mode) for different approaches. The S/N is clearly scale dependent, in contradiction to assumption (A1). We use ℓ = k R_{⊙} to convert from wave number into harmonic degree. We do not show results for ℓ = 0 because this corresponds to an arbitrary change in reference frame. 

In the text 
Fig. 5.
Calibration curve of Hanasoge et al. (2012, Fig. 2 in the supplementary material) can be reproduced (blue solid curve). When a more realistic width σ_{z′} = 16.8 Mm = 1.8 H_{P} of flow features is used to compute the calibration curve, a change is seen (dashed orange curve), which would lead to an increase in the estimated upper limit, see assumption (A6). The crosses and pluses indicate the grid points in harmonic degree, which were obtained by converting horizontal Fourier modes from Cartesian to spherical geometry (see Appendix A.4). The kernels and thus the calibration curves shown here were computed on a spatial grid with a four times larger spatial extent than the MURaM grid, therefore the resolution in Fourier space, , is four times finer than in other plots in this paper. See Fig. C.2 for a more detailed comparison. 

In the text 
Fig. 6.
Test of the validity of the assumptions that enter the derivation of the calibration curve, shown as a noisereduced signal in travel times (left panel) and flow amplitude (center and right panels) as a function of scale. If all assumptions (A2)–(A5) were satisfied, the different curves in each panel would not cross each other, but be ordered from higher to lower values in the same sense as they are ordered from top to bottom in the legends (e.g., Eq. (12) ≥ Eq. (13) ≥ Eq. (15) ≥ Eq. (17) for the left panel). 

In the text 
Fig. 7.
Comparison of estimates of the spectrum of convective flows obtained from the synthetic data based on flows from a MURaM simulation. The method of HDS2012 applied to this synthetic data (solid blue line) can be improved by using a scaledependent signaltonoise correction (dashed orange line, this work). Both estimates give the correct order of magnitude of the flow (green line), although the estimate from the method of HDS2012 is not a strict upper limit. When the underlying assumptions are improved (solid orange line), our best estimate results in a larger upper limit. 

In the text 
Fig. C.1.
Example realization of the horizontal flow component u_{x} at the surface (z = 0). The displayed snapshot is taken in the middle of the time period relevant for the travel times shown in Fig. 2 and is saturated at 33% of its maximum absolute value. Example realizations for u_{y} are statistically similar but transposed, and examples for u_{z} can be found in Lord (2014, Figs. 2.2 and 2.3). 

In the text 
Fig. C.2.
Comparison of calibration curves in a more detailed way than in Fig. 5. When the kernels are interpolated onto the MURaM grid, the results do not substantially change (dotted curves), except at ℓ ≈ 20. The crosses and plusses indicate the grid points in harmonic degree, which were obtained by converting horizontal Fourier modes from Cartesian into spherical geometry (see Appendix A.4). The kernels were computed on a spatial grid with a four times larger spatial extent than the MURaM grid, therefore the resolution in Fourier space, , is four times finer. 

In the text 
Fig. C.3.
Contribution to the noiseless traveltime spectrum from below 0.99R_{⊙} (solid blue line) and from a Gaussianweighted region around the target depth 0.96R_{⊙} (dashed blue) and their forward models based on the actual flow and two versions of the calibration curve (solid orange and dashdotted green). 

In the text 
Fig. C.4.
Azimuthal averages as defined in Eqs. (C.1) and (C.3) saturated at 10% of the maximum value of the left panel. The units are s^{4} cm^{−4} (i.e., cgs units). The quantities are zero for d ≠ d′, see also Fig. C.5. 

In the text 
Fig. C.5.
Exemplary dependence of kernels and flow correlations on the direction of the horizontal Fourier vector. KLeft panel: , right panel: . For the plot, each quantity was divided by its maximum value in the given range. Error bars of C^{u} were obtained from the standard deviation of the 11 available simulation snapshots and are similar for C^{K}C^{u}. At other depths and similar wave numbers, the plots are qualitatively similar. Right panel: we also show that taking the mean before or after multiplying and changes the result (thick solid black line). 

In the text 
Fig. C.6.
Azimuthal averages as defined in Eqs. (C.1) and (C.2) saturated at 10% of the maximum value of the top left panel. The units are s^{2} cm^{2} (i.e., cgs units). 

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.