Free Access
Issue
A&A
Volume 567, July 2014
Article Number A137
Number of page(s) 20
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201423580
Published online 30 July 2014

© 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, vj, 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 vivj 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 p1-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)ht, where ht is the temporal resolution at observation times tn = nht 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 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): (1)Working in Fourier space is faster (and easier). In the time-domain, the cross-covariance becomes (2)where tn is the correlation time lag.

Cross-covariances are the basic data to compute the travel times. We denote τ+(x1,x2) as the travel time for a wave packet traveling from point x1 to point x2 and τ(x1,x2) as the travel time for a wave packet traveling from x2 to x1. In the limit discussed by Gizon & Birch (2004), the incremental travel times can be measured from the estimated cross-covariance using (3)where Cref is a deterministic reference cross-covariance coming from spatial averaging or from a solar model. Then the weight function W± is defined as (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 Cref(x1,x2,t) = Cref(x2x1,t), which implies that W(x1,x2,t) = W(x2x1,t). However, this assumption is not necessary in the remainder of this paper.

We write τα, where the subscript 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(x1,x2) and products of travel times τα1(x1,x2)τα2(x3,x4), where τ is defined by Eq. (3). To simplify the notations, let (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 x1,...,xM are described by a vector-valued stationary Gaussian time series (φ(x1,tn),...(xM,tn)). For the sake of simplicity, we can also assume without loss of generality that E[φ(xm,tn)] = 0 at each xm 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).

thumbnail 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.

Open with DEXTER

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 the expectation value of the cross-covariance (6)

3.2. Independance of the geometry

Gizon & Birch (2004) assumed that the observations φ(xij,tn) are given on a Cartesian grid { xij } 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 (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 (8)was assumed. We have denoted 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 jl 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 and a smoothed version of . 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 (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, (10)where γ is the half width at half maximum of the Lorentzian centered at a frequency ω0. In this case, one can check that the bound in Eq. (9) is equal to 1 / (4π2γ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 /T2 that can be neglected for long observation times.

4. Model noise covariances

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

  • Cov [ τ1,τ2 ] and Cov [ C1,C2 ] which are linked to the fourth order moment of φT,

  • Cov [ τ1τ2,τ3 ] and Cov [ C1C2,C3 ], 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 [ τ1τ2,τ3τ4 ] and Cov [ C1C2,C3C4 ], which depend on the eighth order moment of φT.

For the covariance between two complex random variables X and Y, we use the convention (11)In particular, as the mean value of the observables is zero, we have

We show that all moments of cross-covariance functions depend on only. Because the travel time measurement procedure is linear in C, the moments of the travel-times can be expressed in terms of and of the weight functions Wi (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 (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 (13)where means that the additional terms decay at least as 1 /Tm + 1 (m corresponds to the regularity, which is the number of derivatives of the functions 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 X2 is derived in Appendices B and D. If the observation time T is so small that X2/T2 cannot be neglected, X2 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, (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 xi, 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 vx and vy, 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), (15)where Z1, Z2, and Z3 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 Zi depend only on 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 Zi in Eq. (15). The term of order T-1 is given by (for details, see Appendix E) (16)where the covariance between two travel times is given by Eq. (13), and is the expectation value of the travel time τj. For example, (17)As Cref and 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 (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 Z3 of order T-3, we introduce a function Γα1,α2 such that which according to Eq. (13), is (19)Then the term of order T-3 is given by (20)where μ = { μ1,μ2,··· ,μ6 } and the subset contains all μ satisfying (21) contains 12 elements, so the term Z3 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

As Cref is generally chosen as an average value of the observations, we have or at least . If there is equality, then we can simplify the formula given in the previous section because . It follows that the term Z1 is zero, as are some elements of Z2. Denoting by , the value of Z2 when , we have (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: (23)

4.3.2. Case

Suppose now that we do not have equality but , where ϵ is a small parameter measuring the difference between the reference cross-covariance and their expectation value. In this case, Z1 is of order ϵ2 and the terms that are cancelled out previously in Z2 when are of order ϵ. The numerical tests from Sect. 5.6.1 confirm that these terms of order ϵ and ϵ2 can be neglected, so that Eq. (23) can be used even if we just have .

thumbnail 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.

Open with DEXTER

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 Z2 will dominate for long observation times. In this case, we have the simplified formula: (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 p1 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 p1 ridge in this case. The quantity can be estimated from observations by averaging over a set of (quiet-Sun) filtered power spectra |φ(k)|2. 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 Nx×Nx × N = 512 × 512 × 610 that are centered around the central meridian at a latitude of 40°. The spatial sampling is hx = 0.35 Mm in both directions, and the temporal sampling is ht = 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 hkR = 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 p1 ridge only. In this paper, we consider only one filter for the sake of simplicity. The function is estimated by averaging |φ(k)|2 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 hkR = 24.5 and hω/ 2π = 34.7μHz.

thumbnail 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 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 .

Open with DEXTER

thumbnail 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.

Open with DEXTER

thumbnail 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.

Open with DEXTER

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: (25)where Var [ τ ] = Cov [ τ,τ ] is the theoretical variance for travel times computed by Eq. (13) and Varn [ τ ] is the variance obtained by Monte Carlo simulations with n realizations. Similarly, we define (26)where Var [ τ2 ] = Cov [ τ2,τ2 ] is the theoretical variance for a product of travel times computed by Eq. (15).

Figure 3 shows the errors Err1(n) for Var [ τdiff ] and Err2(n) for for travel times between two points separated by a distance Δ = 10 Mm. As expected, we have (27)with constants depending on the type of measurement. Even if the rate of convergence is the same for τdiff or , 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 . 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 [ τ+(x1,x2),τ+(x3,x4) ] 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 [ τ+(x1,x2)τ+(x3,x4),τ+(x5,x6)τ+(x7,x8) ], as presented in Fig. 6.

thumbnail 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.

Open with DEXTER

thumbnail 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.

Open with DEXTER

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 vxvy, 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).

thumbnail 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.

Open with DEXTER

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 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 Cref

Let us first consider a fixed observation time (T = 8 h for the numerical examples) and look at the dependence on the term Cref. This dependence is due to the term Z1 and one part of Z2, which depends on . Figure 8 makes this comparison for a product of travel times 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 . If we suppose that , then we have (cf. Appendix F):

Thus, even if ϵ is not small, the term Z1 and the second part of the term Z2 are smaller than the other ones in the far field, as . This is confirmed in Fig. 8, where all the terms are plotted in the worst case, when Cref = 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 seem to be always smaller than the other ones, even when Cref = 0. Thus, as discussed in Sect. 4.3, when Cref is close to and T is not too small, it is a good approximation to neglect the terms containing and thus to use Eq. (23) to compute the noise covariance matrix.

thumbnail 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.

Open with DEXTER

thumbnail Fig. 10

Left: 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.

Open with DEXTER

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 Z1 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 (Cref = 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 Z1 is almost zero, so this term can be neglected independently of the observation time. In the first configuration, the term Z3 is always at least two decades smaller than 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 Z3 is dominant. When the observation time is greater than four hours, it then varies as T-2 and 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 Tc = Z3/Z2. For this test case, a good approximation can be found in the far field as presented in Appendix F, where it is shown that Tc ≈ 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, 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 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 (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 = l2 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 182/ (7 / 4)2 = 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.

Table 1

Var [ τ1τ2 ] and Var [ ⟨ τ1τ2A ] (in s4) with l = 18 Mm for the product of a EW and NS travel time (configuration #2 with d = 0).

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 /T2, and 1 /T3, where T is the total observation time. For standard applications of time-distance helioseismology, we showed that the term in 1 /T2 is dominant: 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 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.

Acknowledgments

The authors acknowledge research funding by Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1 “Astrophysical flow instabilities and turbulence” (Project A1, “Solar turbulent convection probed by helioseismology”).

References

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 tj = htj, j ∈Z in this paper, to avoid some technical difficulties. As a consequence, the frequency variable ω is 2π/ht 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 ht → 0: (A.1)We need the following: an orthogonal projection DN of L2( [ − π/ht/ht ]) onto the space ΠN of 2π/ht-periodic trigonometric polynomials with a degree N and the Dirichlet kernel ; the Fejér smoothing operator FN: L2( [ − π/ht/ht ] ) → ΠN with the Fejér kernel N; and the projected periodic Hilbert transform HN: L2( [ − π/ht/ht ] ) → ΠN with kernel N. They are defined by Here, sgn(k): = 1 and sgn(−k): = −1 for k ∈N, and sgn(0): = 0. The transform HN is related to the standard periodic Hilbert transform H with a convolution kernel ℋ(ω) = cot(ω/ 2) by HN = HDN = DNH. With our convention for the Fourier transform, the Fourier convolution theorem is . In particular (with f(ω) = ℱN(htω) and , etc.), we have (A.2)To simplify the notations, the cross-covariance (respectively its expectation value) C(xa,xb) (respectively ) can be simply written as Cab(ω) (resp. ), and similarly, the weight functions W(xa,xb) can be Wab(ω). We show the following theorem on the correlation function Cab(ω1,ω2) defined as (A.3)The covariance between the wavefield at two frequencies ω1 and ω2 can be expressed as The second term is bounded by (A.7)For a stationary Gaussian time series, the error of the approximate noise model (Eq. (8) in Gizon & Birch 2004) is bounded by (A.8)The proof of the above theorem is given below.

By the definition of φT, the covariance betweeen the observations is given by (A.9)where j = kl, m = −l, and gω(j) = ∑ |m|,|jm| ≤ Neiωht(jm) = ∑ |m|,|j + m| ≤ Ne− iωhtm. For ω1 = ω2, we have g0(j) = 2N + 1 − |j|, so Eq. (A.4) for this case follows from Eq. (A.2). We now consider the case ω1ω2. For j> 0, we have If tj< 0, then . Inserting the expression for gω in Eq. (A.9), using the identity sin(xy) = sinxcosy − cosxsiny for x = T(ω2ω1) / 2 and y = ht(ω2ω1) | j | / 2, and finally using Eq. (A.2) leads to To bound IIab, we may assume that | ω2ω1 | ≤ π/ht without loss of generality due to 2π/ht periodicity. Using the mean value theorem, Eq. (A.2), and the inequality for , we obtain (A.10)This yields Eq. (A.7). It also implies Eq. (A.8) for jl, since ; that is Iab(hωj,hωl) = 0. To show Eq. (A.8) for j = l, we use the bound

Appendix B: Frequency correlations for travel times

In this Appendix we derive the noise covariance matrix for the cross-covariance function C and for the travel time τ when the frequency correlations are taken into account. Appendix A has shown that taking the frequency correlations into account leads to an additional term of order 1 /T in the covariance of the observables at the grid points. As the covariance betweeen two travel times is also of order 1 /T, it is of interest to look if this correction should be taken into consideration. This Appendix proves that the extra term in 1 /T of the observable covariance only leads to an additional term in 1 /T2 for the travel times. We also underline the main difficulties that will occur when computing higher order moments of C and τ.

With our convention Eq. (A.1), the Fourier transform is unitary up to the factor . It follows from definition (3)that Therefore, (B.1)The first difficulty is to evaluate the quantity, Cov [ C12(ω1),C34(ω2) ]. For a higher order moment, we also need to evaluate Cov [ C12(ω1)C34(ω2),C56(ω3) ] and Cov [ C12(ω1)C34(ω2),C56(ω3)C78(ω4) ]. The way to deal with these terms is presented in Appendix C, where it is shown that (B.2)It leads to (B.3)The second difficulty comes from the evaluation of these integrals; that is the evalution of linear functionals of the expectation value of the cross-covariance C given by the weight functions W. Similarly, we need to be able to evaluate, for higher order moments, The method to compute these terms is presented in Appendix D. Applying the result for the second order moment presented in Appendix D.1 leads to the following result:

The travel-time covariance for finite T is given by the travel-time covariance for infinite observation time (Eq. (13)) plus a correction that decreases as 1 /T2(B.6)where m corresponds to the regularity (the number of derivatives) of the functions and Wab, and (B.7)

Remark concerning the setting of Gizon & Birch (2004)

In Gizon & Birch (2004), it was assumed that (B.8)so the covariance of C is (B.9)We note that Eq. (B.9) is exact. It differs slightly from Eq. (C8) in Gizon & Birch (2004), which incorrectly contained an additional term. It leads to the covariance between travel times (B.10)We note that Eq. (B.10) is identical to Eq. (28) in Gizon & Birch (2004), as the extra term in the covariance of C was actually neglected by the authors. Taking the frequency correlations into account, Eq. (B.8) is no longer valid, and correction terms have to be added to Eqs. (B.9), (B.10). These correction terms are given in the previous result.

Appendix C: Noise covariance matrix for high order cross-covariances

In this section, we present the way to compute the noise covariance matrices for the cross-covariance function C: As the cross-covariance function can be written as a function of the observables (C.4)Equations (C.1)–(C.3) imply that the moments of 4, 6, and 8 of the observables have to be computed. In the next section we present a formula to compute high order moment of Gaussian variables. Then, we apply this formula to compute Eqs. (C.1)–(C.3).

Appendix C.1: Expectation value of high-order products of Gaussian random variables

We have seen that the moments of order 4, 6 and 8 of the observables have to be computed to find the noise covariance matrix for cross-covariances and products of cross-covariances. A formula to compute the (2J)th-order moment of a multivariate complex normal distribution with zero-mean can be found in Isserlis (1918): (C.5)where μ and ν have distinct values in ˜1,2J¨ and the set J is defined by (C.6)Here, we used the notation ˜1,2J¨ for the set of all integers between 1 and 2J. To better understand Eq. (C.5), let us explain it for the case J = 2. In this case, Eq. (C.5) can be written as (C.7)where the indices i,j,k,l must satisfy i<j, i<k and k<l according to Eq. (C.6). This forces that i = 1. Then, we can have k = 2 or k = 3. If k = 3, then l = 4 and j = 2. If k = 2, then we have again two possibilities: l = 3 and j = 4 or l = 4 and j = 3. So three combinations are possible: (1,2,3,4), (1,4,2,3), and (1,3,2,4). This leads to (C.8)In particular, we have (C.9)which is the formula required to compute the moment of order 4 in Eq. (C.1). For J = 3, Eq. (C.5) becomes (C.10)where the indices i,j,k,l,m,n must satisfy i<k<m (since the sequence (μi) must increase) and i<j, k<l and m<n (since μi<νi), according to Eq. (C.6). Hence, we obtain (C.11)A problem is that the cardinality of the set J is (4J) ! / [ (2J) ! 4J ] (Isserlis 1918) increases exponentially. The sum in Eq. (C.5) contains three terms for J = 2 and 15 for J = 3, as shown above. Unfortunately, it leads to 105 terms for J = 4 so it is not convenient to write them down explicitly, and we just list the main guidelines in Sect. C.4.

Appendix C.2: Second order moment of C

In the original paper, the fourth order moment of the observables was guessed after looking at all the possible cases in the Fourier domain. Using Eq. (C.9) and the definitions of Cab and Cab (Eqs. (C.4), (A.3)) and recalling that as φj(t) is real-valued, the covariance matrix between two cross-covariances is readily computed as follows: (C.12)

Appendix C.3: Third order moment of C

In this section, we compute the sixth order moment of the observables defined by Eq. (C.2). After writing the cross-correlations as a function of the observables, we need to compute the moment of order 6 of the observables. This can be done using Eq. (C.11) with z1 = φ1(ω1), z2 = φ2(ω1), z3 = φ3(ω2), z4 = φ4(ω2), z5 = φ5(ω3), and z6 = φ6(ω3). After integration against weight functions, it will turn out that the order of the different terms in 1 /T depends on their degree of separability. Therefore, we denote by the sum of the terms, which can be written as product of at most N functions of disjoint subsets of the set of variables { ω1,ω2,ω3 }. Then (C.13)where (C.14)and (C.15)

Appendix C.4: Fourth order moment of C

This section is devoted to the computation of the eigth order moment of the observables defined by Eq. (C.3). Writing the cross-correlations as a function of the observables leads to (C.16)In Eq. (C.16) and in the following we omit the argument ωj of the observables φ2j − 1 = φ2j − 1(ωj) and φ2j = φ2j(ωj). As for the moments of order 4 and 6, we can calculate this expression. As explained in Sect. C.1, the moment of order 8 contains 105 terms, so we do not write explicitely all the terms. As for the moments of order 6, we arrange the terms as (C.17)where is the the sum of all terms, which can be written as a product of at most N functions of disjoint subsets of the set of variables { ω1,ω2,ω3,ω4 }. The three terms will be computed below.

Expression for

These terms are the ones from the subset given by Eq. (C.6) from which the observables use the same frequencies in two expectation values, for example . It leads to the following formula: (C.18)Calculating all the expectation values implies which can be written in terms of the covariance between two cross-covariance functions,

Expression for

Two kinds of products in Eq. (C.5) will lead to terms with only two frequency integrals:

  • In two expectation values, the constraints on ω are the same; for example, (they lead to the first two terms in Eq. (C.19));

  • In one expectation value, the observables use the same frequencies; for example, .

Computing all the terms, one can show that (C.19)The terms Cov [ C,C ] and Cov [ CC,C ] that appear in this expression can be computed using Eqs. (C.12), (C.13).

Expression for

All other terms lead to terms that contain only one frequency integral in the covariance of the product of travel times. After reorganizing all the terms, one can show that can be written as

Appendix D: Evaluation of separable linear functionals of nonseparable products of Cab’s

In this section, we derive asymptotic expansions of the terms, in 1 /T as T → ∞ and explicit formulae for the leading order terms. Recall that C defined in Eq. (A.3) depends on T, although this is suppressed in our notation.

Appendix D.1: Functionals of nonseparable products of two Cab functions

In this subsection, we show that (D.4)where is defined by Eq. (B.7) if and have m derivatives and W12 and W34 have m − 1 derivatives.

Plugging Eq. (A.4) into the left hand side of Eq. (D.4), we arrive at a sum (2π)2(X + 2Y + Z) involving the following three terms: We repeatedly use the following transformation of variables formula for functions f(ω1,ω2), which are 2π/ht periodic in both variables: (D.8)Even though the Jacobian of this transformation of variables is 1 / 2, no factor appears, since we integrate over a domain that can be reassembled to two periodicity cells on the right hand side.

Using Eq. (D.8) and noting that , the first term can be written as We want to interpret the inner product as a convolution with 2N evaluated at 0. First, we note that by a change of variables: . Let f(ω1,ω2) be 2π/ht periodic in both arguments and . Then and hence (D.9)where F2N always acts on the second argument. As , it follows that Since , we get an additional if we omit the orthogonal projections D2N in the last equation.

To bound Y (Eq. (D.6)), we again apply the change of variables in Eq. (D.8) to obtain where f has uniformly bounded derivatives of order m − 1. When T tends to infinity, this corresponds to a high order Fourier coefficient and thus can be made as small as desired. In particular, by repeated partial integration, (D.10)The term Z (Eq. (D.7)) can be transformed, in the same way, and after using that , we find that where the higher order term comes from which is similar to Eq. (D.10). As limn → ∞H2Nf = Hf and all the terms in the integrals are bounded it follows that X is of order 1 /T2. Gathering the expressions for the three terms X, Y, and Z leads to Eq. (D.4).

Appendix D.2: Functionals of nonseparable products of three Cab functions

Let C be defined by Eq. (A.3), and Wi represent some functions of ω. Then, we have the following expension: (D.11)Using Eq. (A.4) in the left hand side of Eq. (D.11), four different types of terms have to be studied

where the expressions I and II are given by Eqs. (A.5), (A.6) respectively.

We use the change of variables: (D.16)where the Jacobian 1 / 3 does not appear for the same reason as in Eq. (D.8). Applying this to X, we obtain where ωi can be replaced by the corresponding value in . The role of the Fejér kernel is played by the function: where we have used the change of variables m = jl, n = k + l, and o = l. If denotes the corresponding convolution operator and is the two-dimensional orthogonal projection, we can use the inequality max( | m | , | n | , | mn | ) ≤ | m | + | n | to obtain (D.17)If f(ω1,ω2,ω3) is 2π/ht periodic in all of its arguments and , we find in analogy to Sect. D.1 that and hence . With Eq. (D.17), we obtain (D.18)The terms Y1 is of very high order using the same method than in Sect. D.1. The term Y2 can be treated in the same way as it also contains a cosine that oscillates with T. Finally, Z is of order 1 /T3 using a similar demonstration than in Sect. D.1.

Appendix D.3: Functionals of nonseparable products of four Cab functions

Let C be defined by Eq. (A.3), and Wi represent some functions of ω. Then, we have the following expension: (D.19)As in the previous proof, different terms have to be treated. The terms with combinations of the expressions I and II can be bounded by the same methods as in Sect. D.2, and the term involving only expressions II can be bounded as in Sect. D.1. The only different term is (D.20)Here, small adaptions of the argument in Section D.2 with the change of variables, (D.21)lead to the formula (D.22)

Appendix E: Noise covariance matrix for products of travel times

Appendix E.1: Third order moment of the travel times

Using the definition of the travel times, we obtain that the covariance for the product of travel times is given by: Using Eq. (C.13) and the two results presented in Sects. D.1 and D.2, we can express the covariance for three travel-times as (E.1)where is the expectation value of τj and the covariance involving two travel times can be computed with Eq. (13).

Appendix E.2: Analytic formula for the covariance matrix for products of travel times

In this section, we derive the main result of this paper. It gives an analytic expression for the covariance matrix between a product of travel times. Using the definition of the travel times, one can show that the covariance of the product of travel times is given by (E.2)In Appendix D, we have shown that not all the terms lead to the same number of frequency integrals. This implies that the covariance given by Eq. (E.2) has terms of different order with respect to the observation time T. The terms containing three integrals in ω are of order T-1, while the other ones are of order T-2 and T-3. We write the covariance as the sum between three terms for the different orders: (E.3)The terms of order 1 /T4 come from the correlation between the frequencies in the frequency domain as detailed in Sect. B for the covariance between travel times. The other terms are detailed below.

Term Z1 of order T-1

Looking at Eq. (E.2), one can see that this term is composed of

  • all the terms involving Cov [ C,C ],

  • the terms with two integrals in ω for the terms with Cov [ CC,C ] (term ),

  • the terms with three integrals in ω for the terms with Cov [ CC,CC ] (term ),

where C is a generic cross-covariance. Reorganizing terms leads to the formula Eq. (16) for Z1.

Term Z2 of order T-2

Looking at Eq. (E.2) one can see that this term is composed of

  • the terms with one integral in ω for the terms with Cov [ CC,C ] (term ),

  • the terms with two integrals in ω for the terms with Cov [ CC,CC ] (term ).

Reorganizing terms leads to the formula Eq. (18) for Z2.

Term Z3 of order T-3

The terms of order T-3 come from the terms with only one integral in ω in Cov [ CC,CC ] (term ). This yields Eq. (20) for Z3.

Appendix F: Far-field approximation for

In this section, we give approximate expressions for the different terms that compose Eq. (15) for in the far field (Δ → ∞). We start with the definitions of Z1, Z2, and Z3:

In the far field, we have . If we suppose that , then the global behavior of the four terms is (F.1)We can thus see that the global behavior of the terms is

As we can conclude that the first term in Z2 and the one in Z3 are dominant in this case. We can go further to see for which observation time Tc these two last terms intersect in the case of difference travel times. If the window function f(t) in the definition of Wdiff, as defined by Eq. (4), is a Heavyside function, then we have (Gizon & Birch 2004) (F.2)For a p-mode ridge κr = κr(ω), the function can be written in the far field as Gizon & Birch (2004): (F.3)where κi is the imaginary part of the wavenumber at resonance and represents attenuation of the waves. The sums in Eq. (F.1) can be approximated because the cosine in Eq. (F.3) oscillates many times within the frequency width ξ of the envelope of : Using the numerical value of ξ/ 2π = 1mHz, the observation time Tc at which the two terms are equal is (F.4)For T>Tc, Z2/T2 is the dominant term. As the observation time is traditionally of at least eight hours in helioseismology, the term of order 1 /T3 can be neglected.

All Tables

Table 1

Var [ τ1τ2 ] and Var [ ⟨ τ1τ2A ] (in s4) with l = 18 Mm for the product of a EW and NS travel time (configuration #2 with d = 0).

All Figures

thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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 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 .

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 10

Left: 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.