© ESO, 2014

1. Introduction

The purpose of time-distance helioseismology (Duvall et al. 1993; Gizon & Birch 2005, and references therein) is to infer the subsurface structure and dynamics of the Sun using spatial-temporal correlations of the random wave field observed at the solar surface. Wave travel times between pairs of points (denoted τ) are measured from the cross-covariance function. Wave speed perturbations and vector flows are then obtained by inversion of the travel times (e.g. Kosovichev 1996; Jackiewicz et al. 2012). Such inversions require knowledge of the noise covariance matrix Cov[τ,τ]. Typically, noise is very high and strong correlations exist among travel times. Gizon & Birch (2004) studied the noise properties of travel times and derived a simple noise model that successfully explains the observations. The model is based on the assumption that the stochastic noise is stationary and horizontally spatially homogeneous, as a result of the excitation of waves by turbulent convection. In addition to time-distance helioseismology, this noise model has found applications in direct modeling inversions (Woodard 2006, 2009) and ring-diagram analysis (Birch et al. 2007).

Time-distance helioseismology has been successfully applied to map flow velocities, v_j, at supergranulation scales (Kosovichev 1996; Duvall & Gizon 2000; Gizon et al. 2001; Jackiewicz et al. 2008). The statistical properties of convection can further be studied by computing horizontal averages of the turbulent velocities. For example, Duvall & Gizon (2000); Gizon et al. (2010) showed that the horizontal divergence and the vertical vorticity of the flows are correlated through the influence of the Coriolis force on convection. It would be highly desirable to extract additional properties of the turbulent velocities; for example, the (anisotropic) Reynolds stresses ⟨ v_iv_j ⟩ that control the global dynamics of the Sun (differential rotation and meridional circulation, see Kitchatinov & Rüdiger 2005). The noise associated with such measurement involves the fourth order moments of the travel times, Cov [ ττ,ττ ].

Alternatively, we would like to consider spatial averages of products of travel times ⟨ ττ ⟩ as the fundamental data from which to infer the Reynolds stresses (or other second-order moments of turbulence). Spatial averages are meaningful when turbulent flows are horizontally homogeneous over the averaging region. Inversions of average products of travel times are desirable, since input data are fewer and less noisy. Once again, we need to know the noise covariance matrix Cov [ ⟨ ττ ⟩ , ⟨ ττ ⟩ ] to perform the inversion.

In this paper, we study the noise properties of travel times and products of travel times. In Sect. 2, the definitions for the cross-covariance function and the travel times are given. Section 3 presents the assumptions of the noise model generalizing the model of Gizon & Birch (2004). In Sect. 4 and the Appendices, we derive analytical formulae for the noise covariance matrices of travel times and products of travel times. These formulae are confirmed in Sect. 5 by comparison to numerical Monte Carlo simulations and to observations of the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO). The effects of horizontal spatial averaging are considered in Sect. 6.

2. Observables: cross-covariance function, travel times, and products of travel times

The fundamental observation in helioseismology is the filtered line-of-sight Doppler velocity φ(x,t) at points x on the surface of the Sun and at times t. The filter acts by multiplication in the Fourier domain. In this paper, we only consider the p₁-ridge filter as an example. We note that all the results presented in this paper do not depend on the choice of the filter. The signal φ(x,t) is recorded over a duration time T = (2N + 1)h_t, where h_t is the temporal resolution at observation times t_n = nh_t for n = − N,...,N. The observed wavefield during the observation time T is denoted φ_T. We have , where is a window function (equal to 1 if |t| ≤ T/ 2 and 0 otherwise).

Helioseismic analysis is performed in Fourier space. Let us define the temporal Fourier transform of φ_T by $\begin{array}{ccc}{\mathit{\phi }}_{\mathit{T}}\mathrm{\left(}{x}\mathit{,\omega }\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}\sum _{\mathit{n}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}\mathit{\phi }\mathrm{\left(}{x}\mathit{,}{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\exp \mathrm{\left(}\mathit{i\omega }{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\mathit{.}& & \end{array}$The frequencies ω are treated as continuous variables in the remainder of this paper to be able to take the frequency correlations into account (see Sect. 3.3). The cross-covariance function between two points at the surface of the Sun is a multiplication in the Fourier domain (Duvall et al. 1993): $\mathit{C}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{T}}{\mathit{\phi }}_{\mathit{T}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,\omega }\mathrm{\right)}{\mathit{\phi }}_{\mathit{T}}\mathrm{\left(}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\mathit{.}$(1)Working in Fourier space is faster (and easier). In the time-domain, the cross-covariance becomes $\mathit{C}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{N}\mathrm{+}\mathrm{1}}\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\max \mathrm{(}\mathrm{-}\mathit{N,}\mathrm{-}\mathit{N}\mathrm{-}\mathit{n}\mathrm{)}}^{\min \mathrm{(}\mathit{N,N}\mathrm{-}\mathit{n}\mathrm{)}}\mathit{\phi }\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}\mathit{\phi }\mathrm{\left(}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{t}}_{\mathit{j}\mathrm{+}\mathit{n}}\mathrm{\right)}\mathit{,}$(2)where t_n is the correlation time lag.

Cross-covariances are the basic data to compute the travel times. We denote τ₊(x₁,x₂) as the travel time for a wave packet traveling from point x₁ to point x₂ and τ₋(x₁,x₂) as the travel time for a wave packet traveling from x₂ to x₁. In the limit discussed by Gizon & Birch (2004), the incremental travel times can be measured from the estimated cross-covariance using $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{\pm }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathrm{\right)}& \mathrm{:}\mathrm{=}\hspace{0.17em}& {\mathit{h}}_{\mathit{t}}\sum _{\mathit{n}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}{\mathit{W}}_{\mathrm{\pm }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\\ & & \mathrm{\times }\mathrm{(}\mathit{C}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\mathrm{-}{\mathit{C}}^{\mathrm{ref}}\mathrm{(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{t}}_{\mathit{n}}{\mathrm{)}}^{\mathrm{)}}\mathit{,}\end{array}$(3)where C^ref is a deterministic reference cross-covariance coming from spatial averaging or from a solar model. Then the weight function W_± is defined as ${\mathit{W}}_{\mathrm{\pm }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,t}\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathrm{\mp }\mathit{f}\mathrm{(}\mathrm{\pm }\mathit{t}\mathrm{)}{\mathit{\partial }}_{\mathit{t}}{\mathit{C}}^{\mathrm{ref}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,t}\mathrm{\right)}}{{\mathit{h}}_{\mathit{t}}{\sum }_{\mathit{n}}\mathit{f}\mathrm{(}\mathrm{\pm }{\mathit{t}}_{\mathit{n}}\mathrm{)}\mathrm{\left[}{\mathit{\partial }}_{\mathit{t}}{\mathit{C}}^{\mathrm{ref}}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{t}}_{\mathit{n}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}$(4)with f a window function used to select an interval of time around the first arrival time of the wave packet (for example, a cut-off function). For the spatially homogeneous noise, notice that we generally choose that C^ref(x₁,x₂,t) = C^ref(x₂ − x₁,t), which implies that W(x₁,x₂,t) = W(x₂ − x₁,t). However, this assumption is not necessary in the remainder of this paper.

We write τ_α, where the subscript $\begin{array}{ccc}\mathit{\alpha }\mathrm{\in }\mathrm{\{}\mathrm{+}\mathit{,}\mathrm{-}\mathit{,}\mathrm{diff}\mathit{,}\mathrm{mean}\mathrm{\}}& & \end{array}$denotes the type of travel time and the corresponding weight function W_α. The mean and difference travel times τ_diff and τ_mean can be obtained from the one-way travel times by τ_diff = τ₊ − τ₋ and τ_mean = (τ₊ + τ₋) / 2.

In this paper, we are interested in the noise covariance matrix for travel times τ_α1(x₁,x₂) and products of travel times τ_α1(x₁,x₂)τ_α2(x₃,x₄), where τ is defined by Eq. (3). To simplify the notations, let $\begin{array}{ccc}& & {\mathit{\tau }}_{\mathrm{1}}\mathrm{:}\mathrm{=}{\mathit{\tau }}_{{\mathit{\alpha }}_{\mathrm{1}}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathrm{\right)}\mathit{,} {\mathit{\tau }}_{\mathrm{2}}\mathrm{:}\mathrm{=}{\mathit{\tau }}_{{\mathit{\alpha }}_{\mathrm{2}}}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathrm{\right)}\mathit{,}\\ & & \mathrm{andmoregenerally}{\mathit{\tau }}_{\mathit{i}}\mathrm{:}\mathrm{=}{\mathit{\tau }}_{{\mathit{\alpha }}_{\mathit{i}}}\mathrm{\left(}{{x}}_{\mathrm{2}\mathit{i}\mathrm{-}\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}\mathit{i}}\mathrm{\right)}\mathit{.}\end{array}$(5)

3. Generalization of the noise model

3.1. Assumptions

The basic assumption of the noise model is the following: The observations at the relevant spatial points x₁,...,x_M are described by a vector-valued stationary Gaussian time series (φ(x₁,t_n),...,φ(x_M,t_n)). For the sake of simplicity, we can also assume without loss of generality that E[φ(x_m,t_n)] = 0 at each x_m for all n ∈Z. This model is valid in the quiet Sun (away from evolving active regions) but does not assume that the noise is spatially homogeneous contrary to the model of Gizon & Birch (2004) as detailed in Sect. 3.2. This assumption is supported by the observed distribution of the HMI Doppler velocity: Fig. 1 shows the probability density of the filtered line-of-sight velocity. For a Gaussian distribution, the data should line up along a straight line. We can see a very good agreement for probabilities betweeen 5% and 95%. The deviations in the tail of the plot (for probabilities smaller than 5%) may be due to statistical errors (as we have less realizations for these events).

Fig. 1 Probability density plot representing the filtered line-of-sight velocity φ(t) for a p1-ridge with an observation time T = 8 h. For Gaussian observations, all data should be on a straight line.

One may also replace the spatial points by some spatial averages. Such averages are often used to improve the signal-to-noise ratio. We denote $\overline{\mathit{C}}$ the expectation value of the cross-covariance $\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,}{{x}}_{\mathit{b}}\mathit{,\omega }\mathrm{\right)}\mathrm{=}E\left[\mathit{C}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,}{{x}}_{\mathit{b}}\mathit{,\omega }\mathrm{\right)}\right]\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{T}}E\left[{\mathit{\phi }}_{\mathit{T}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,\omega }\mathrm{\right)}{\mathit{\phi }}_{\mathit{T}}\mathrm{\left(}{{x}}_{\mathit{b}}\mathit{,\omega }\mathrm{\right)}\right]\mathit{.}$(6)

3.2. Independance of the geometry

Gizon & Birch (2004) assumed that the observations φ(x_ij,t_n) are given on a Cartesian grid { x_ij } by an approximately flat patch of the Sun’s surface. The discrete Fourier transform of the finite dimensional signal was assumed to be of the form $\mathit{\phi }\mathrm{\left(}{{k}}_{\mathit{ij}}\mathit{,}{\mathit{\omega }}_{\mathit{l}}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{𝒫}\mathrm{\left(}{{k}}_{\mathit{ij}}\mathit{,}{\mathit{\omega }}_{\mathit{l}}\mathrm{\right)}}{\mathrm{𝒩}}_{\mathit{ijl}}\mathit{,}$(7)where is the power spectrum, ω_l: = 2πl/T, and are complex independent and identically distributed Gaussian variables with zero-mean and unit variance. In this case, the frequency correlations were ignored and $\frac{\mathrm{2}\mathit{\pi }}{\mathit{T}}E\mathrm{\left[}{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\right(}{{x}}_{\mathit{a}}\mathit{,}{\mathit{\omega }}_{\mathit{j}}\mathrm{\left)}\mathit{\phi }\mathrm{\right(}{{x}}_{\mathit{b}}\mathit{,}{\mathit{\omega }}_{\mathit{l}}{\mathrm{\left)}}^{\mathrm{\right]}}\mathrm{=}{\mathit{\delta }}_{\mathit{jl}}{\overline{\mathit{C}}}_{\mathrm{GB}}\mathrm{(}{{x}}_{\mathit{b}}\mathrm{-}{{x}}_{\mathit{a}}\mathit{,}{\mathit{\omega }}_{\mathit{l}}\mathrm{)}$(8)was assumed. We have denoted ${\overline{\mathit{C}}}_{\mathrm{GB}}$ the expectation value of the cross-covariance used by Gizon & Birch (2004). Our assumption is more general as it does not require a planar geometry and allows a natural treatment of spatially averaged quantities. It means that all our results are valid in any geometry, and it is in particular the case for the results presented in Gizon & Birch (2004).

3.3. On frequency correlations

As the observation time T is finite, the discrete Fourier transforms φ_T(x,ω_j) and φ_T(x,ω_l) for j ≠ l are no longer uncorrelated because of the window function. The necessity of a correction term for finite T was discussed but not further analyzed in Gizon & Birch (2004). It turns out that there is an explicit formula for this correction term in terms of the periodic Hilbert transform of $\overline{\mathit{C}}$ and a smoothed version of $\overline{\mathit{C}}$. The exact formulation is given in Appendix A, where it is also shown that the error made by considering a finite observation time can be bounded by $\begin{array}{ccc}& & \underset{\mathit{j,l}}{\sup }\left|\frac{\mathrm{2}\mathit{\pi }}{\mathit{T}}E\mathrm{[}{\mathit{\phi }}_{\mathit{T}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{\phi }}_{\mathit{T}}\mathrm{(}{{x}}_{\mathit{b}}\mathit{,}{\mathit{\omega }}_{\mathit{l}}{\mathrm{)}}^{\mathrm{]}}\mathrm{-}{\mathit{\delta }}_{\mathit{jl}}{\overline{\mathit{C}}}_{\mathrm{GB}}\mathrm{(}{{x}}_{\mathit{b}}\mathrm{-}{{x}}_{\mathit{a}}\mathit{,}{\mathit{\omega }}_{\mathit{l}}\mathrm{)}\mathrm{]}\right|\\ & & \mathrm{\le }\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{4}\mathit{T}}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\left|{\mathit{t}}_{\mathit{k}}\right|\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,}{{x}}_{\mathit{b}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathit{.}\end{array}$(9)Note that the right hand side of (9)depends only on T and on a quantity depending on the correlation length of the waves. This can be better seen using an analytic cross-covariance given by a Lorentzian of the form, $\overline{\mathit{C}}\mathrm{\left(}{x}\mathit{,}{x}\mathit{,\omega }\mathrm{\right)}\mathrm{=}\frac{{\mathit{C}}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\omega }\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{\gamma }}^{\mathrm{2}}}\mathit{,}$(10)where γ is the half width at half maximum of the Lorentzian centered at a frequency ω₀. In this case, one can check that the bound in Eq. (9) is equal to 1 / (4π²γT). Therefore, the correlations between frequencies should only be neglected when this bound is small, meaning that the observation time is long enough to represent correctly the mode.

As the covariance between travel times is known to be also of order 1 /T (Gizon & Birch 2004), it is legitimate to wonder if the frequency correlations should be taken into account. It is shown below (see Eq. (13)) that the consideration of frequency correlations only leads to additional terms of order 1 /T² that can be neglected for long observation times.

4. Model noise covariances

In this section and Appendices B–E, we present explicit formulae for the covariance matrices of cross-covariances C, travel times τ, and products of cross-covariances or travel-times:

• Cov [ τ₁,τ₂ ] and Cov [ C₁,C₂ ] which are linked to the fourth order moment of φ_T,

• Cov [ τ₁τ₂,τ₃ ] and Cov [ C₁C₂,C₃ ], which require the knowledge of the sixth order moment of φ_T and are necessary to compute the moment of order four of τ and C,

• Cov [ τ₁τ₂,τ₃τ₄ ] and Cov [ C₁C₂,C₃C₄ ], which depend on the eighth order moment of φ_T.

For the covariance between two complex random variables X and Y, we use the convention $\mathrm{Cov}\mathrm{\left[}\mathit{X,Y}\mathrm{\right]}\mathrm{=}E\left[\mathit{X}{\mathit{Y}}^{\mathrm{\ast }}\right]\mathrm{-}E\left[\mathit{X}\right]E\left[{\mathit{Y}}^{\mathrm{\ast }}\right]\mathit{.}$(11)In particular, as the mean value of the observables is zero, we have $\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{T}}\mathrm{Cov}\mathrm{\left[}{\mathit{\phi }}_{\mathit{T}}\mathrm{\right(}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\left)}\mathit{,}{\mathit{\phi }}_{\mathit{T}}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,\omega }\mathrm{\left)}\mathrm{\right]}\mathit{.}$

We show that all moments of cross-covariance functions depend on $\overline{\mathit{C}}$ only. Because the travel time measurement procedure is linear in C, the moments of the travel-times can be expressed in terms of $\overline{\mathit{C}}$ and of the weight functions W_i (see Eq. (4)).

4.1. Covariance matrix for C and travel times

As a first step, we show in Appendix C that the covariance between two cross-correlations is given by $\begin{array}{ccc}& & {\left(\frac{\mathit{T}}{\mathrm{2}\mathit{\pi }}\right)}^{\mathrm{2}}\mathrm{Cov}\mathrm{\left[}\mathit{C}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}\mathit{,C}\mathrm{\right(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}\\ & & \mathrm{=}E\mathrm{\left[}{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}\mathit{\phi }\mathrm{\right(}{{x}}_{\mathrm{3}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}E\mathrm{\left[}\mathit{\phi }\mathrm{\right(}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\right(}{{x}}_{\mathrm{4}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}\\ & & \mathrm{+}E\mathrm{\left[}{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\right(}{{x}}_{\mathrm{4}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}E\mathrm{\left[}\mathit{\phi }\mathrm{\right(}{{x}}_{\mathrm{2}}\mathit{,}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}\mathit{\phi }\mathrm{\right(}{{x}}_{\mathrm{3}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}\mathit{.}\end{array}$(12)For a comparison with a small correction to the corresponding formula in Gizon & Birch (2004), we refer to Appendix B. The covariance between two travel times is given by $\begin{array}{ccc}& & \mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{2}}\mathrm{\right]}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}{\mathit{T}}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}^{\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}\mathrm{d}\mathit{\omega }{\mathit{W}}_{{\mathit{\alpha }}_{\mathrm{1}}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\\ & & \mathrm{\times }({\mathit{W}}_{{\mathit{\alpha }}_{\mathrm{2}}}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{3}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{4}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\\ & & \mathrm{+}{\mathit{W}}_{{\mathit{\alpha }}_{\mathrm{2}}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)})\mathrm{+}\frac{{\mathit{X}}_{\mathrm{2}}}{{\mathit{T}}^{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{T}}^{\mathit{m}\mathrm{+}\mathrm{1}}}\right)\mathit{,}\end{array}$(13)where means that the additional terms decay at least as 1 /T^{m + 1} (m corresponds to the regularity, which is the number of derivatives of the functions $\overline{\mathit{C}}$ and W). A good agreement between the leading order term in this formula and SOHO MDI measurements was found by Gizon & Birch (2004). An explicit formula for the second order term X₂ is derived in Appendices B and D. If the observation time T is so small that X₂/T² cannot be neglected, X₂ can easily be evaluated numerically.

4.2. Covariance matrix for products of travel times

In this section, we are interested in the covariance matrix for the travel times correlations, which is to evaluate the quantity, $\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathrm{\left)}{\mathit{\tau }}_{\mathrm{2}}\mathrm{\right(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathrm{\left)}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right(}{{x}}_{\mathrm{5}}\mathit{,}{{x}}_{\mathrm{6}}\mathrm{\left)}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right(}{{x}}_{\mathrm{7}}\mathit{,}{{x}}_{\mathrm{8}}\mathrm{\left)}\mathrm{\right]}\mathit{.}$(14)This quantity is the most general we can evaluate for velocity correlations. It will be helpful to derive all the formulae in more specific frameworks. In general, this quantity depends on the eight points x_i, but it is of course possible to look at simpler cases. For example, we may be interested in the correlations between a east-west (EW) and north-south (NS) travel time as presented in Fig. 4. This quantity can give us information about the correlations between the velocities v_x and v_y, which are velocities in the EW and NS directions, respectively.

The formula for the product of cross-covariances is given in Appendix E (Eq. (C.17)) and is not be discussed in the text, where we focus on products of travel times. In Appendix E, we derive the general formula for Eq. (14), $\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{=}\frac{\mathrm{1}}{\mathit{T}}{\mathit{Z}}_{\mathrm{1}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{2}}}{\mathit{Z}}_{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{3}}}{\mathit{Z}}_{\mathrm{3}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{4}}}\right)\mathit{,}$(15)where Z₁, Z₂, and Z₃ are given by Eqs. (16), (18), and (20) and is detailed later after some general remarks on this formula. An important point is that all the terms in Z_i depend only on $\overline{\mathit{C}}$ and on the weight functions W. Thus, it is possible to directly estimate the noise covariance matrix via this formula instead of performing a large number of Monte-Carlo simulations. This strategy is much more efficient as we see in Sect. 5.3, where we demonstrate the rate of convergence of the stochastic simulations.

The terms on the right hand side of the general formula Eq. (15) are of different orders with respect to the observation time. The behavior of these terms is studied in Sect. 5.6.2.

Let us now give the expressions for the different terms Z_i in Eq. (15). The term of order T^-1 is given by (for details, see Appendix E) $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{T}}{\mathit{Z}}_{\mathrm{1}}& \mathrm{=}& {\overline{\mathit{\tau }}}_{\mathrm{2}}({\overline{\mathit{\tau }}}_{\mathrm{4}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{3}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]})\\ & & \mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{1}}({\overline{\mathit{\tau }}}_{\mathrm{4}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{3}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]})\mathit{,}\end{array}$(16)where the covariance between two travel times is given by Eq. (13), and ${\overline{\mathit{\tau }}}_{\mathit{j}}$ is the expectation value of the travel time τ_j. For example, ${\overline{\mathit{\tau }}}_{\mathrm{1}}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}^{\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}\mathrm{d}\mathit{\omega }{\mathit{W}}_{\mathrm{1}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\mathrm{\left(}\overline{\mathit{C}}\mathrm{\right(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{)}\mathrm{-}{\mathit{C}}^{\mathrm{ref}}\mathrm{(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{.}$(17)As C^ref and $\overline{\mathit{C}}$ are generally close or even equal it is possible that this quantity is close to 0 or even exactly 0. This simplification is discussed in Sect. 4.3. We note that the time dependence (in T^-1) in Eq. (16) is hidden on the right hand side in the covariance between two travel times (cf. Eq. (13)).

The term of order T^-2 is given by $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{2}}}{\mathit{Z}}_{\mathrm{2}}& \mathrm{=}& \mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\\ & & \mathrm{-}{\overline{\mathit{\tau }}}_{\mathrm{1}}\left(\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{3}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{4}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\right)\\ & & \mathrm{-}{\overline{\mathit{\tau }}}_{\mathrm{2}}\left(\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{3}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{4}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\right)\\ & & \mathrm{-}{\overline{\mathit{\tau }}}_{\mathrm{3}}\left(\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{1}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}{\overline{\mathit{\tau }}}_{\mathrm{2}}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\right)\\ & & \end{array}$(18)where the covariance involving three travel times is given in the Appendix E.2 by Eq. (E.1) and the one between two travel times by Eq. (13). As we see in Sect. 5, the first line of this term is dominant in most of the applications.

Before writing down the term Z₃ of order T^-3, we introduce a function Γ_α1,α2 such that $\begin{array}{ccc}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{2}}\mathrm{\right]}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}{\mathit{T}}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}^{\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}\mathrm{d}\mathit{\omega }{\mathrm{\Gamma }}_{{\mathit{\alpha }}_{\mathrm{1}}\mathit{,}{\mathit{\alpha }}_{\mathrm{2}}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{T}}^{-2}\mathrm{\right)}\mathit{,}& & \end{array}$which according to Eq. (13), is $\begin{array}{ccc}& & {\mathrm{\Gamma }}_{{\mathit{\alpha }}_{\mathrm{1}}\mathit{,}{\mathit{\alpha }}_{\mathrm{2}}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\\ & & \mathrm{=}{\mathit{W}}_{{\mathit{\alpha }}_{\mathrm{1}}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathrm{\right)}({\mathit{W}}_{{\mathit{\alpha }}_{\mathrm{2}}}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{3}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{4}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)}\\ & & \mathrm{+}{\mathit{W}}_{{\mathit{\alpha }}_{\mathrm{2}}}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{4}}\mathit{,\omega }\mathrm{\right)}\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathrm{3}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,\omega }\mathrm{\right)})\mathit{.}\end{array}$(19)Then the term of order T^-3 is given by $\begin{array}{ccc}{\mathit{Z}}_{\mathrm{3}}& \mathrm{=}& \frac{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{7}}}{{\mathit{T}}^{\mathrm{3}}}\sum _{\mathit{\mu }\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\mathrm{ℳ}}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}^{\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}\mathrm{d}\mathit{\omega }{\mathrm{\Gamma }}_{{\mathit{\alpha }}_{\mathrm{1}}\mathit{,}{\mathit{\alpha }}_{{\mathit{\mu }}_{\mathrm{1}}}}\mathrm{\left(}{{x}}_{\mathrm{1}}\mathit{,}{{x}}_{\mathrm{2}}\mathit{,}{{x}}_{{\mathit{\mu }}_{\mathrm{1}}}\mathit{,}{{x}}_{{\mathit{\mu }}_{\mathrm{2}}}\mathit{,\omega }\mathrm{\right)}\\ & & \mathrm{\times }{\mathrm{\Gamma }}_{{\mathit{\alpha }}_{{\mathit{\mu }}_{\mathrm{3}}}\mathit{,}{\mathit{\alpha }}_{{\mathit{\mu }}_{\mathrm{5}}}}\mathrm{\left(}{{x}}_{{\mathit{\mu }}_{\mathrm{3}}}\mathit{,}{{x}}_{{\mathit{\mu }}_{\mathrm{4}}}\mathit{,}{{x}}_{{\mathit{\mu }}_{\mathrm{5}}}\mathit{,}{{x}}_{{\mathit{\mu }}_{\mathrm{6}}}\mathit{,\omega }\mathrm{\right)}\mathit{,}\end{array}$(20)where μ = { μ₁,μ₂,··· ,μ₆ } and the subset ℳ contains all μ satisfying $\{\begin{array}{c}\\ {\mathit{\mu }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathit{<}{\mathit{\mu }}_{\mathrm{2}}\mathrm{if}{\mathit{\mu }}_{\mathrm{1}}\mathrm{odd}\\ {\mathit{\mu }}_{\mathrm{1}}\mathit{<}{\mathit{\mu }}_{\mathrm{2}}\mathrm{if}{\mathit{\mu }}_{\mathrm{1}}\mathrm{even}\\ {\mathit{\mu }}_{\mathrm{3}}\mathit{<}{\mathit{\mu }}_{\mathrm{4}}\mathit{<}{\mathit{\mu }}_{\mathrm{5}}\mathit{<}{\mathit{\mu }}_{\mathrm{6}}\\ {\mathit{\mu }}_{\mathrm{1}}\mathit{,}\mathit{...}{\mathit{\mu }}_{\mathrm{6}}\mathrm{\in }\mathrm{3}\mathit{,}\mathit{...}\mathit{,}\mathrm{8}\mathit{.}\end{array}$(21)ℳ contains 12 elements, so the term Z₃ consists of a sum of 12 terms containing a product of the functions Γ defined by Eq. (19).

4.3. Important special cases

4.3.1. Case ${{C}}^{\mathrm{ref}}\mathrm{=}\overline{{C}}$

As C^ref is generally chosen as an average value of the observations, we have ${\mathit{C}}^{\mathrm{ref}}\mathrm{=}\overline{\mathit{C}}$ or at least ${\mathit{C}}^{\mathrm{ref}}\mathrm{\approx }\overline{\mathit{C}}$. If there is equality, then we can simplify the formula given in the previous section because $\overline{\mathit{\tau }}\mathrm{=}\mathrm{0}$. It follows that the term Z₁ is zero, as are some elements of Z₂. Denoting by $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}$, the value of Z₂ when ${\mathit{C}}^{\mathrm{ref}}\mathrm{=}\overline{\mathit{C}}$, we have $\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{2}}}\begin{array}{c}\mathrm{˜}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}\mathrm{=}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{+}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathit{.}$(22)This term is of order T^-2, as each of the covariance in Eq. (22) are of order T^-1. The noise covariance matrix is now given by the sum of two terms of order T^-2 and T^-3: $\begin{array}{ccc}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}& & \mathrm{=}\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{2}}}\begin{array}{c}\mathrm{˜}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}\mathrm{+}\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{3}}}{\mathit{Z}}_{\mathrm{3}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{4}}}}^{\mathrm{)}}\mathrm{·}\end{array}$(23)

4.3.2. Case ${{C}}^{\mathrm{ref}}\mathrm{\approx }\overline{{C}}$

Suppose now that we do not have equality but ${\mathit{C}}^{\mathrm{ref}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{ϵ}\mathrm{)}\overline{\mathit{C}}$, where ϵ is a small parameter measuring the difference between the reference cross-covariance and their expectation value. In this case, Z₁ is of order ϵ² and the terms that are cancelled out previously in Z₂ when ${\mathit{C}}^{\mathrm{ref}}\mathrm{=}\overline{\mathit{C}}$ are of order ϵ. The numerical tests from Sect. 5.6.1 confirm that these terms of order ϵ and ϵ² can be neglected, so that Eq. (23) can be used even if we just have ${\mathit{C}}^{\mathrm{ref}}\mathrm{\approx }\overline{\mathit{C}}$.

Fig. 2 Average p1 power spectrum obtained from SDO/HMI dopplergrams. The sampling is given by hkR⊙ = 24.5 and hω/ 2π = 34.7μHz. Left panel: cut at frequency ω/ 2π = 3.4 mHz. Right panel: cut at ky = 0. The dark parts correspond to large values of the power spectrum and the white ones to small values.

4.3.3. Simplified formula

We have now defined all the terms involved in Eq. (15) to compute the covariance of a product of travel times. As one term is of order T^-2 and the other one of order T^-3, it follows that Z₂ will dominate for long observation times. In this case, we have the simplified formula: $\begin{array}{ccc}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}& \mathrm{=}\hspace{0.17em}& \mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\\ & & \mathrm{+}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathit{.}\end{array}$(24)In the next section, we show applications of this formula, which validate the model and the simplified formula. In particular, the numerical tests tell us that Eq. (24) can be used if the observation time is more than roughly a few hours.

5. Examples and comparisons

5.1. SDO/HMI power spectrum for p₁ ridge

In this section, we validate the analytic formulae for the noise by comparing with Monte Carlo simulations. We choose to use a homogeneous noise, so the model depends only on the expectation value of the power spectrum, . This expectation value is computed in the Fourier domain to perform filtering to only keep the p₁ ridge in this case. The quantity can be estimated from observations by averaging over a set of (quiet-Sun) filtered power spectra |φ(k,ω)|². Here, we consider observations of the line-of-sight Doppler velocity from the HMI instrument on board of the SDO spacecraft (Schou et al. 2012) between 6 April 2012 and 14 May 2012. We prepare Postel-projected datacubes of size N_x×N_x × N = 512 × 512 × 610 that are centered around the central meridian at a latitude of 40^°. The spatial sampling is h_x = 0.35 Mm in both directions, and the temporal sampling is h_t = 45 s. The physical size of the data-cube is L × L × T = 180 Mm × 180 Mm × 8 h. The sampling in Fourier space is given by h_kR_⊙ = 24.5 and h_ω/ 2π = 34.7μHz.

The filtered wave field, φ, is obtained by applying a filter in 3D Fourier space that lets through the p₁ ridge only. In this paper, we consider only one filter for the sake of simplicity. The function is estimated by averaging |φ(k,ω)|² over forty 8 h data cubes separated by one day. In Fig. 2, we show cuts through the average power spectrum.

5.2. Monte Carlo simulations

We use the expectation value of the observed power spectrum defined above as input to the noise model. To validate the theoretical model, we run Monte Carlo simulations by generating many realizations of the wave field in Fourier space using Eq. (7). The normal distributions are generated with the ziggurat algorithm of MATLAB (Marsaglia & Tsang 1984). All realizations have the same dimensions as above, which are h_kR_⊙ = 24.5 and h_ω/ 2π = 34.7μHz.

Fig. 3 Convergence of the numerical simulations to the model for a p1-ridge with an observation time T = 8 h. The errors Erri(n) defined by Eqs. (25, 26) are represented for Var [ τdiff ] and $\mathrm{Var}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{diff}}^{\mathrm{2}}\mathrm{\right]}$ for travel times between two points separated by a distance Δ = 10 Mm. The dashed lines have a slope of 1/2 and show that the error decays as ${\mathit{n}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}$.

Fig. 4 Geometrical configuration #1: geometry used for the covariance between a east-west and north-south travel time Cov [ τ1,τ2 ], where τ1 = τα1(x1,x2) and τ2 = τα2(x3,x4). The distance between x1 and x2 and between x3 and x4 is Δ = 10 Mm.

Fig. 5 Cov [ τ1,τ2 ] (in s2) for a p1-ridge at a latitude of 40° with an observation time T = 8 h in configuration #1 given by Fig. 4. A travel time τ+ is used for τ1 and τ2. Left: SDO/HMI observations, middle: Monte Carlo simulation, and right: analytic formula.

5.3. Rate of convergence toward the analytic formula

To show the importance of having an explicit formula for the noise, we look at the convergence of Monte Carlo simulations to the analytic formula. For this, we define the following measure of the error: ${\mathrm{Err}}_{\mathrm{1}}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{=}\frac{\left|\mathrm{Var}\mathrm{\left[}\mathit{\tau }\mathrm{\right]}\mathrm{-}{\mathrm{Var}}_{\mathit{n}}\mathrm{\left[}\mathit{\tau }\mathrm{\right]}\right|}{\mathrm{Var}\mathrm{\left[}\mathit{\tau }\mathrm{\right]}}\mathit{,}$(25)where Var [ τ ] = Cov [ τ,τ ] is the theoretical variance for travel times computed by Eq. (13) and Var_n [ τ ] is the variance obtained by Monte Carlo simulations with n realizations. Similarly, we define ${\mathrm{Err}}_{\mathrm{2}}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{=}\frac{\left|\mathrm{Var}\mathrm{\left[}{\mathit{\tau }}^{\mathrm{2}}\mathrm{\right]}\mathrm{-}{\mathrm{Var}}_{\mathit{n}}\mathrm{\left[}{\mathit{\tau }}^{\mathrm{2}}\mathrm{\right]}\right|}{\mathrm{Var}\mathrm{\left[}{\mathit{\tau }}^{\mathrm{2}}\mathrm{\right]}}\mathit{,}$(26)where Var [ τ² ] = Cov [ τ²,τ² ] is the theoretical variance for a product of travel times computed by Eq. (15).

Figure 3 shows the errors Err₁(n) for Var [ τ_diff ] and Err₂(n) for $\mathrm{Var}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{diff}}^{\mathrm{2}}\mathrm{\right]}$ for travel times between two points separated by a distance Δ = 10 Mm. As expected, we have ${\mathrm{Err}}_{\mathit{i}}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{\approx }{\mathrm{const}}_{\mathit{i}}\hspace{0.17em}{\mathit{n}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}$(27)with constants depending on the type of measurement. Even if the rate of convergence is the same for τ_diff or ${\mathit{\tau }}_{\mathrm{diff}}^{\mathrm{2}}$, the constant is much smaller for a travel time than for a product of travel times. The variance of a product of travel times converges much slower than the travel time variance. For example, an accuracy of 5% is reached with about n = 1000 realizations for τ_diff but around n = 5000 for ${\mathit{\tau }}_{\mathrm{diff}}^{\mathrm{2}}$. This underlines the importance of having an analytic formula to obtain the correct limit when n → ∞, especially in the case of products of travel times.

5.4. Noise of travel times: Comparison with Monte-Carlo simulations and SDO/HMI observations

To show the level of noise in the data, we compare the noise matrix with HMI data from 6 April 2012 until 14 May 2012. The point-to-point travel times are obtained for a distance Δ = 10 Mm in the x and y direction so that we can compare Cov [ τ₊(x₁,x₂),τ₊(x₃,x₄) ] in the configuration given by Fig. 4. The comparison between the data, Monte Carlo simulation, and the explicit formula is given in Fig. 5. As expected, data contain mainly noise as we are looking only at point-to-point travel times, and a good agreement is found between stochastic simulations and the analytic formula.

5.5. Noise of products of travel times: Comparison with Monte-Carlo simulations and SDO/HMI observations

In the previous section, we show that the data are dominated by noise in the case of point-to-point travel times, so it is legitimate to ask if there is information in a product of travel times. We look at the covariance between two products of EW and NS travel times, Cov [ τ₊(x₁,x₂)τ₊(x₃,x₄),τ₊(x₅,x₆)τ₊(x₇,x₈) ], as presented in Fig. 6.

Fig. 6 Geometrical configuration #2: geometry used for the covariance between a product of east-west and north-south travel times Cov [ τ1τ2,τ3τ4 ], τi = ταi(x2i − 1,x2i) are defined in Eq. (5) . The travel distance between pairs of points is Δ = 10 Mm.

Fig. 7 Cov [ τ1τ2,τ3τ4 ] (in s4) for a p1-ridge at a latitude of 40° with an observation time T = 8 h in configuration #2 given by Fig. 6. A travel time τ+ is used for the four travel times. Left: SDO/HMI observations, middle: theory, and right: cut through dy = 0 to compare SDO/HMI observations, theory, and Monte Carlo simulations.

The results are given in Fig. 7. As previously for travel times, we note a good agreement between the analytic formula and the Monte Carlo simulation. In this case, one can see the differences between the noise and the data, which are separated by around 2σ. To confirm that this difference is due to the presence of physical signal (supergranulation) and not to a problem in the model, we show the same covariance in Fig. 9, at the equator, instead of at a latitude of 40^°. In this case, data, analytic formula and Monte Carlo simulations fit perfectly. Since the product ⟨ τ_xτ_y ⟩ (configuration with d = 0) measures the Reynolds stress ⟨ v_xv_y ⟩, it is expected to be zero at the equator and non-zero away from the equator (as we observe).

For both latitudes, the correlation length is identical and equal to λ/ 4 where λ = 7 Mm is the dominant wavelength of the filtered wave field. This is half of the correlation length for travel times, as one can see with the simplified formula Eq. (24).

Fig. 8 Comparison of the three terms in Eq. (15) for the variance of a product of travel times separated by a distance Δ. The comparison is done for a p1-ridge and an observation time of T = 8 h.

5.6. Test of simplified formula for products of travel times using Monte Carlo simulations

We have shown that some simplifications can be made to the analytic formula for the noise covariance matrix if ${\mathit{C}}^{\mathrm{ref}}\mathrm{=}\overline{\mathit{C}}$ in Sect. 4.3. In this section, we show numerically that these simplifications can be done even if we do not have equality and that Eq. (24) is a good approximation for the noise covariance matrix.

5.6.1. Sensitivity to choice of C^ref

Let us first consider a fixed observation time (T = 8 h for the numerical examples) and look at the dependence on the term C^ref. This dependence is due to the term Z₁ and one part of Z₂, which depends on $\overline{\mathit{\tau }}$. Figure 8 makes this comparison for a product of travel times ${\mathit{\tau }}_{\mathrm{diff}}^{\mathrm{2}}$ between points separated by Δ. In this simple case, it is possible to write down the global behavior of the different terms in the far field, when $\overline{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathit{,\omega }\mathrm{\right)}\mathit{<}\overline{\mathit{C}}\mathrm{\left(}\mathrm{0}\mathit{,\omega }\mathrm{\right)}$. If we suppose that ${\mathit{C}}^{\mathrm{ref}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{ϵ}\mathrm{)}\overline{\mathit{C}}$, then we have (cf. Appendix F):

$\begin{array}{ccc}& & \frac{\mathrm{1}}{\mathit{T}}{\mathit{Z}}_{\mathrm{1}}\mathrm{\sim }\frac{{\mathit{ϵ}}^{\mathrm{2}}}{\mathit{T}}\overline{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathit{,\omega }{\mathrm{)}}^{\mathrm{2}}\overline{\mathit{C}}\mathrm{\right(}\mathrm{0}\mathit{,\omega }{\mathrm{)}}^{\mathrm{2}}\\ & & \frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{2}}}{\mathit{Z}}_{\mathrm{2}}\mathrm{\sim }\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{2}}}\overline{\mathit{C}}\mathrm{(}\mathrm{0}\mathit{,\omega }{\mathrm{)}}^{\mathrm{4}}\mathrm{+}\frac{\mathit{ϵ}}{{\mathit{T}}^{\mathrm{2}}}\overline{\mathit{C}}\mathrm{(}\mathrm{0}\mathit{,\omega }{\mathrm{)}}^{\mathrm{3}}\overline{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathit{,\omega }\mathrm{\right)}\\ & & \frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{3}}}{\mathit{Z}}_{\mathrm{3}}\mathrm{\sim }\frac{\mathrm{1}}{{\mathit{T}}^{\mathrm{3}}}\overline{\mathit{C}}\mathrm{(}\mathrm{0}\mathit{,\omega }{\mathrm{)}}^{\mathrm{4}}\mathit{.}\end{array}$Thus, even if ϵ is not small, the term Z₁ and the second part of the term Z₂ are smaller than the other ones in the far field, as $\overline{\mathit{C}}\mathrm{\left(}\mathrm{\Delta }\mathit{,\omega }\mathrm{\right)}\mathit{<}\overline{\mathit{C}}\mathrm{\left(}\mathrm{0}\mathit{,\omega }\mathrm{\right)}$. This is confirmed in Fig. 8, where all the terms are plotted in the worst case, when C^ref = 0. Results are similar for the test cases using the configuration #2, so we did not plot them. Even if the simplifications presented above are only applicable for this particular test case, the terms containing $\overline{\mathit{\tau }}$ seem to be always smaller than the other ones, even when C^ref = 0. Thus, as discussed in Sect. 4.3, when C^ref is close to $\overline{\mathit{C}}$ and T is not too small, it is a good approximation to neglect the terms containing $\overline{\mathit{\tau }}$ and thus to use Eq. (23) to compute the noise covariance matrix.

Fig. 9 Cov [ τ1τ2,τ3τ4 ] (in s4) for a p1-ridge at the equator with an observation time T = 8 h in configuration #2 given by Fig. 6. This is a cut through dy = 0, which compares SDO/HMI observations, theory, and Monte Carlo simulations.

Fig. 10 Left: $\mathrm{Var}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{diff}}^{\mathrm{2}}\mathrm{\right]}$ as a function of the observation time with Δ = 20 Mm. Right: comparison of the three terms in Eq. (15) for the variance between a product of east-west and north-south travel time with Δ = 20 Mm.

5.6.2. Dependence on observation duration T

The formula giving the covariance for a product of travel times (Eq. (15)) contains three terms that behave differently as a function of the observation time T. It is thus interesting to compare these terms to see if some can be dropped or if some are dominant. The term Z₁ is initially kept to ensure that the dependence on the observation time does not make this term become significant. As previously noted, we suppose that we have no knowledge about a reference cross-covariance (C^ref = 0). Figure 10 makes this comparison for the variance in configuration #1 and the covariance in configuration #2 as a function of T (with Δ = 20 Mm). We see that the contribution of the term Z₁ is almost zero, so this term can be neglected independently of the observation time. In the first configuration, the term Z₃ is always at least two decades smaller than $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}$ and so only this last term can be kept. The situation is sligthly different for the second configuration. When T is smaller than one hour, then the standard deviation varies as T^-3 and the term Z₃ is dominant. When the observation time is greater than four hours, it then varies as T^-2 and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}$ is dominant. If T is very long, then the variations should be in T^-1. This area happens theoretically for an observation time longer than two months, which is not realistic for solar applications and is thus not shown in Fig. 10. The intersection between both terms is given by T_c = Z₃/Z₂. For this test case, a good approximation can be found in the far field as presented in Appendix F, where it is shown that T_c ≈ 100 min, which is confirmed numerically in Fig. 10. These comparisons of the different terms are extremely important, as it implies that we can use the approximation given by Eq. (22) if we consider observation times of a few hours which is generally the case. If the observation time is shorter, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}$ is still a good approximation and gives a good estimate of the noise even if the amplitude is not exact. It is certainly sufficient to use $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{Z}}_{\mathrm{2}}\end{array}$ as a noise covariance matrix to perform an inversion but numerical tests still have to be performed.

6. Spatial averages

We define the average value of a quantity q over an area A as follows $\mathrm{⟨}\mathit{q}{\mathrm{⟩}}_{\mathit{A}}\mathrm{=}\frac{\mathrm{1}}{\mathit{A}}{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}\sum _{{x}\mathrm{\in }\mathit{A}}\mathit{q}\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}$(28)The noise covariance matrix for averaged travel times and products of travel times can be obtained by integrating Eq. (13) and Eq. (15) respectively. Averaging data has the advantage of increasing the signal-to-noise ratio and allows us to deal with fewer data. Table 1 shows the accuracy of the analytic formula and the importance of the averaging. It compares the value of the variance for a product betweeen east-west and north-south travel times (configuration #2 with d = 0) and the same variance when the quantities are averaged over a domain A = l² with l = 18 Mm. First of all, we note good agreement between the analytic formula and the Monte Carlo simulations. Second, the value of the variance is reduced by a factor 100 when we average the product of travel times over the spatial domain. As expected, the variance decreases with the number of independent realizations which is the area A divided by square of the correlation length λ/ 4 (see Sect. 5.5) that is 18²/ (7 / 4)² = 105. Finally, the signal to noise ratio increases with the averaging, and we can see a difference due to physical signal between the observations and the noise model.

7. Conclusions

In this paper, we presented two main generalizations of the noise model of Gizon & Birch (2004) for helioseismic travel times. First, the assumption of spatial homogeneity has been dropped. This is useful in modeling noise in regions of magnetic activity (sunspots and active regions), where oscillation amplitudes are significantly reduced, and in noise across the solar disk as at different center-to-limb distances. Second, we generalized the noise model to higher-order moments of the travel times, in particular products of travel times. We showed that the covariance matrix for products of travel times consists of three terms that scale like 1 /T, 1 /T², and 1 /T³, where T is the total observation time. For standard applications of time-distance helioseismology, we showed that the term in 1 /T² is dominant: $\begin{array}{ccc}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}& \mathrm{=}\hspace{0.17em}& \mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\\ & & \mathrm{+}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{1}}\mathit{,}{\mathit{\tau }}_{\mathrm{4}}\mathrm{\right]}\mathrm{Cov}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{2}}\mathit{,}{\mathit{\tau }}_{\mathrm{3}}\mathrm{\right]}\mathit{.}\end{array}$This very simple formula links the noise covariance of products of travel times to the covariance of travel times and depends only on the expectation value of the cross-covariance $\overline{\mathit{C}}\mathrm{\left(}{x}\mathit{,\omega }\mathrm{\right)}$ and can be obtained directly from the observations. The model is accurate and computationally efficient. It compares very well with Monte Carlo simulations and SDO/HMI observations. The analytic formulae presented in this paper can be used to compute the noise covariance matrices for averaged quantities and thus increase the signal-to-noise ratio. Finally, we would like to emphasize that our results (moments of order 4, 6, and 8 of the wavefield φ(x,ω)) can be extended to modeling noise for other methods of local helioseismology, such as ring-diagram analysis, holography, or far-side imaging.

Online material

Appendix A: Frequency correlations for the observables

In this Appendix, we study the correlations in frequency space that result from a finite observation duration T. First, we collect some definitions.

Since observations are discrete, we only consider discrete time points t_j = h_tj, j ∈Z in this paper, to avoid some technical difficulties. As a consequence, the frequency variable ω is 2π/h_t periodic. However, our definitions of the discrete Fourier transform, and its inverse are chosen such that we obtain the time-continuous case in the limit h_t → 0: $\mathrm{𝒫}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{t}}_{\mathit{j}}}\mathrm{𝒫}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}\mathit{,} \mathrm{𝒫}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}^{\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }{\mathit{t}}_{\mathit{k}}}\mathrm{𝒫}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\omega }\mathit{.}$(A.1)We need the following: an orthogonal projection D_N of L²( [ − π/h_t,π/h_t ]) onto the space Π_N of 2π/h_t-periodic trigonometric polynomials with a degree ≤N and the Dirichlet kernel ; the Fejér smoothing operator F_N: L²( [ − π/h_t,π/h_t ] ) → Π_N with the Fejér kernel ℱ_N; and the projected periodic Hilbert transform H_N: L²( [ − π/h_t,π/h_t ] ) → Π_N with kernel ℋ_N. They are defined by $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$Here, sgn(k): = 1 and sgn(−k): = −1 for k ∈N, and sgn(0): = 0. The transform H_N is related to the standard periodic Hilbert transform H with a convolution kernel ℋ(ω) = cot(ω/ 2) by H_N = HD_N = D_NH. With our convention for the Fourier transform, the Fourier convolution theorem is $\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}{\sum }_{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathit{f}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{𝒫}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{t}}_{\mathit{k}}}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}^{\mathit{\pi }\mathit{/}{\mathit{h}}_{\mathit{t}}}\mathit{f}\mathrm{(}\mathit{\omega }\mathrm{-}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\omega }\end{array}\mathrm{)}\mathrm{𝒫}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\omega }\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\omega }\end{array}$. In particular (with f(ω) = ℱ_N(h_tω) and $\mathit{f}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{h}}_{\mathit{t}}}\frac{\mathit{N}\mathrm{+}\mathrm{1}\mathrm{-}\mathrm{\left|}\mathit{k}\mathrm{\right|}}{\mathit{N}\mathrm{+}\mathrm{1}}$, etc.), we have $\mathrm{\left(}{\mathit{D}}_{\mathit{N}}\mathrm{𝒫}\mathrm{\right)}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}\mathrm{𝒫}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{t}}_{\mathit{k}}}\mathit{,} \mathrm{\left(}{\mathit{F}}_{\mathit{N}}\mathrm{𝒫}\mathrm{\right)}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}\frac{\mathit{N}\mathrm{+}\mathrm{1}\mathrm{-}\mathrm{\left|}\mathit{k}\mathrm{\right|}}{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{𝒫}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{t}}_{\mathit{k}}}\mathit{,} \mathrm{\left(}{\mathit{H}}_{\mathit{N}}\mathrm{𝒫}\mathrm{\right)}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}\frac{\mathit{sgn}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathrm{i}}\mathrm{𝒫}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{t}}_{\mathit{k}}}\mathit{.}$(A.2)To simplify the notations, the cross-covariance (respectively its expectation value) C(x_a,x_b,ω) (respectively $\overline{\mathit{C}}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,}{{x}}_{\mathit{b}}\mathit{,\omega }\mathrm{\right)}$) can be simply written as C_ab(ω) (resp. ${\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$), and similarly, the weight functions W(x_a,x_b,ω) can be W_ab(ω). We show the following theorem on the correlation function C_ab(ω₁,ω₂) defined as $C_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{T}}E\mathrm{\left[}{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\right(}{{x}}_{\mathit{a}}\mathit{,}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}\mathit{\phi }\mathrm{\right(}{{x}}_{\mathit{b}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}\mathit{.}$(A.3)The covariance between the wavefield at two frequencies ω₁ and ω₂ can be expressed as $\begin{array}{ccc}{C}_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \{\begin{array}{c}\\ & \\ & \end{array}\\ \mathrm{where} {\mathit{I}}_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}& \mathrm{:}\mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{T}}{\mathrm{𝒟}}_{\mathit{N}}\mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\left)}\mathrm{\right)}(\mathrm{(}{\mathit{D}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathrm{(}{\mathit{D}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\right)}) \\ \mathit{I}{\mathit{I}}_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}& \mathrm{:}\mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{T}}\frac{\cos \left(\mathit{T}\mathrm{(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\right)}{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}(\mathrm{(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}\mathrm{(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)})\mathit{.} \end{array}$The second term is bounded by $\left|\mathit{I}{\mathit{I}}_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\right|\mathrm{\le }\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{4}}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\frac{\mathrm{\left|}{\mathit{t}}_{\mathit{k}}\mathrm{\right|}}{\mathit{T}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathit{.}$(A.7)For a stationary Gaussian time series, the error of the approximate noise model (Eq. (8) in Gizon & Birch 2004) is bounded by $\left|C_{\mathit{ab}}\mathrm{\left(}{\mathit{h}}_{\mathit{\omega }}\mathit{j,}{\mathit{h}}_{\mathit{\omega }}\mathit{l}\mathrm{\right)}\mathrm{-}{\mathit{\delta }}_{\mathit{j,l}}\mathrm{\left(}{\mathit{D}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\left(\frac{\mathrm{2}\mathit{\pi j}}{\mathit{T}}\right)\right|\mathrm{\le }\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{4}}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\frac{\mathrm{\left|}{\mathit{t}}_{\mathit{k}}\mathrm{\right|}}{\mathit{T}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right| \mathrm{for}\mathit{j,l}\mathrm{\in }Z\mathit{,}\mathrm{\left|}\mathit{j}\mathrm{\right|}\mathit{,}\mathrm{\left|}\mathit{l}\mathrm{\right|}\mathrm{\le }\mathit{N}\mathit{.}$(A.8)The proof of the above theorem is given below.

By the definition of φ_T, the covariance betweeen the observations is given by $\begin{array}{ccc}{C}_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi T}}\sum _{\mathit{l}\mathrm{=}\mathrm{-}\mathit{N}}^{\mathit{N}}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}E\left[{\mathit{\phi }}^{\mathrm{\ast }}\mathrm{\left(}{{x}}_{\mathit{a}}\mathit{,}{\mathit{t}}_{\mathit{l}}\mathrm{\right)}\mathit{\phi }\mathrm{\left(}{{x}}_{\mathit{b}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right]{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{l}}}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{2}}{\mathit{t}}_{\mathit{k}}}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi T}}\sum _{\mathit{l}\mathrm{=}\mathrm{-}\mathit{N}}^{\mathit{N}}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{N}}^{\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{(}{\mathit{t}}_{\mathit{k}}\mathrm{-}{\mathit{t}}_{\mathit{l}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{l}}}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{2}}{\mathit{t}}_{\mathit{k}}} \\ & \mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi T}}\sum _{\mathit{j}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{j}}}\sum _{\mathrm{\left|}\mathit{m}\mathrm{\right|}\mathrm{\le }\mathit{N,}\mathrm{|}\mathit{j}\mathrm{-}\mathit{m}\mathrm{|}\mathrm{\le }\mathit{N}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{j}\mathrm{-}\mathit{m}}}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{2}}{\mathit{t}}_{\mathit{j}\mathrm{-}\mathit{m}}}\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi T}}\sum _{\mathit{j}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}{\mathit{g}}_{{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}}\mathrm{\left(}\mathit{j}\mathrm{\right)}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{j}}}\mathit{,} \end{array}$(A.9)where j = k − l, m = −l, and g_ω(j) = ∑ _{|m|,|j − m| ≤ N}e^{iωh_t(j − m)} = ∑ _{|m|,|j + m| ≤ N}e^{− iωh_tm}. For ω₁ = ω₂, we have g₀(j) = 2N + 1 − |j|, so Eq. (A.4) for this case follows from Eq. (A.2). We now consider the case ω₁ ≠ ω₂. For j> 0, we have $\begin{array}{ccc}{\mathit{g}}_{\mathit{\omega }}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{=}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{N}}^{\mathit{N}\mathrm{-}\mathit{j}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathit{m}}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathit{N}}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}}}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}}}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathrm{(}\mathit{N}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathrm{(}\mathrm{2}\mathit{N}\mathrm{-}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}}}{{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathit{/}\mathrm{2}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathit{/}\mathrm{2}}}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }{\mathit{h}}_{\mathit{t}}\mathit{j}\mathit{/}\mathrm{2}}\frac{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathit{\omega }\mathrm{\right(}\mathrm{2}\mathit{N}\mathrm{-}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}{\sin \mathrm{(}{\mathit{h}}_{\mathit{t}}\mathit{\omega }\mathit{/}\mathrm{2}\mathrm{)}}\mathrm{·}& & \end{array}$If t_j< 0, then ${\mathit{g}}_{\mathit{\omega }}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{=}{\mathit{g}}_{\mathit{\omega }}^{\mathrm{\ast }}\mathrm{(}\mathrm{-}\mathit{j}\mathrm{)}$. Inserting the expression for g_ω in Eq. (A.9), using the identity sin(x − y) = sinxcosy − cosxsiny for x = T(ω₂ − ω₁) / 2 and y = h_t(ω₂ − ω₁) | j | / 2, and finally using Eq. (A.2) leads to $\begin{array}{ccc}{C}_{\mathit{ab}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathit{,}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi T}\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathrm{(}{\mathit{\omega }}_{\mathrm{1}}\mathrm{+}{\mathit{\omega }}_{\mathrm{2}}\mathrm{)}{\mathit{t}}_{\mathit{j}}\mathit{/}\mathrm{2}}\sin \left({\mathit{h}}_{\mathit{t}}\frac{{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}}{\mathrm{2}}\mathrm{(}\mathrm{2}\mathit{N}\mathrm{+}\mathrm{1}\mathrm{-}\left|\mathit{j}\right|\mathrm{)}\right)\\ & \mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi T}\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\left(\sin \left(\frac{{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}}{\mathrm{2}}\mathit{T}\right)\sum _{\mathit{j}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}\mathrm{(}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{j}}}\mathrm{+}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{2}}{\mathit{t}}_{\mathit{j}}}\mathrm{)}\mathrm{-}\cos \left(\frac{{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}}{\mathrm{2}}\mathit{T}\right)\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}\frac{\mathit{sgn}\mathrm{\left(}\mathit{j}\mathrm{\right)}}{\mathit{i}}\mathrm{(}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{2}}{\mathit{t}}_{\mathit{j}}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathrm{1}}{\mathit{t}}_{\mathit{j}}}\mathrm{)}\right)\\ & \mathrm{=}& \frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{T}}\frac{\sin \mathrm{\left(}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{T}\mathit{/}\mathrm{2}\mathrm{)}}{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\mathrm{(}\mathrm{\left(}{\mathit{D}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathrm{\left(}{\mathit{D}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\mathrm{(}{\mathit{\omega }}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{)}}\mathrm{-}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{T}}\frac{\cos \mathrm{\left(}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{T}\mathit{/}\mathrm{2}\mathrm{)}}{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\mathrm{(}\mathrm{\left(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathrm{\left(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\mathrm{(}{\mathit{\omega }}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{)}}\mathit{.}\end{array}$To bound II_ab, we may assume that | ω₂ − ω₁ | ≤ π/h_t without loss of generality due to 2π/h_t periodicity. Using the mean value theorem, Eq. (A.2), and the inequality $\frac{\mathrm{\left|}\mathit{x}\mathrm{\right|}}{\mathrm{\left|}\sin \mathrm{\right(}\mathit{x}\mathrm{\left)}\mathrm{\right|}}\mathrm{\le }\frac{\mathit{\pi }\mathit{/}\mathrm{2}}{\sin \mathrm{(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}$ for $\mathrm{\left|}\mathit{x}\mathrm{\right|}\mathrm{\le }\frac{\mathit{\pi }}{\mathrm{2}}$, we obtain $\begin{array}{ccc}\left|\frac{\mathrm{\left(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathrm{\left(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\right)}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{1}}\mathrm{\right)}}{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\right|& & \mathrm{\le }\frac{\mathrm{|}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{|}}{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\underset{\mathit{\omega }}{\sup }\left|\mathrm{\left(}{\mathit{H}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}{\mathrm{)}}^{\mathrm{\prime }}\mathrm{\right(}\mathit{\omega }\mathrm{)}\right|\\ & & \mathrm{\le }\frac{\mathrm{1}}{\mathit{\pi }}\frac{{\mathit{h}}_{\mathit{t}}\mathrm{|}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{|}\mathit{/}\mathrm{2}}{\sin \mathrm{\left(}{\mathit{h}}_{\mathit{t}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{)}}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\mathrm{\left|}{\mathit{t}}_{\mathit{k}}\mathrm{\right|}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathrm{\le }\frac{\mathrm{1}}{\mathrm{2}}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\mathrm{\left|}{\mathit{t}}_{\mathit{k}}\mathrm{\right|}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathit{.}\end{array}$(A.10)This yields Eq. (A.7). It also implies Eq. (A.8) for j ≠ l, since ${\mathrm{𝒟}}_{\mathit{N}}{}^{\mathrm{\right(}}\frac{\mathrm{2}\mathit{\pi }\mathrm{(}\mathit{j}\mathrm{-}\mathit{l}\mathrm{)}}{\mathrm{2}\mathit{N}\mathrm{+}\mathrm{1}}^{\mathrm{\left)}}\mathrm{=}\mathrm{0}$; that is I_ab(h_ωj,h_ωl) = 0. To show Eq. (A.8) for j = l, we use the bound $\begin{array}{ccc}\left|\mathrm{(}{\mathit{D}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{-}{\mathit{F}}_{\mathrm{2}\mathit{N}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{)}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\right|\mathrm{\le }\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{2}\mathit{\pi }}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\frac{\mathit{k}}{\mathrm{4}\mathit{N}\mathrm{+}\mathrm{2}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathrm{=}\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{4}\mathit{\pi }}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\frac{\mathrm{\left|}{\mathit{t}}_{\mathit{k}}\mathrm{\right|}}{\mathit{T}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathrm{\le }\frac{{\mathit{h}}_{\mathit{t}}}{\mathrm{4}}\left|\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathrm{2}\mathit{N}}^{\mathrm{2}\mathit{N}}\frac{\mathrm{\left|}{\mathit{t}}_{\mathit{k}}\mathrm{\right|}}{\mathit{T}}{\overline{\mathit{C}}}_{\mathit{ab}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\right|\mathit{.}& & \end{array}$