Free Access
Issue
A&A
Volume 620, December 2018
Article Number A136
Number of page(s) 8
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201833825
Published online 10 December 2018

© ESO 2018

1. Introduction

Local helioseismology aims at studying the solar interior in three dimensions by exploiting the information contained in the waves observed at the solar surface (see, e.g., Gizon & Birch 2005 and references therein). Helioseismic holography is one particular approach of local helioseismology, which images subsurface scatterers by back-propagating the surface wave field to target points in the interior. Helioseismic holography is also known as Lindsey-Braun (LB) holography (Lindsey & Braun 1997, 2000a and references therein). It has been used to study solar convection (Braun et al. 2004, 2007), active region emergence (Birch et al. 2013, 2016), sunspot subsurface structure (Braun & Birch 2008b; Birch et al. 2009), to image wave sources (Lindsey et al. 2006), to study sunquakes caused by solar flares (Zharkov et al. 2013; Besliu-Ionescu et al. 2017), and to detect active regions on the far side of the Sun (Lindsey & Braun 2000b; Liewer et al. 2014).

In acoustics, a well-established version of holography in a medium that contains sources is Porter-Bojarski (PB) holography (Porter & Devaney 1982). PB holography uses both the wave field and its normal derivative at the surface to produce holographic images (Porter 1969). PB holography was introduced in helioseismology by Skartlien (2001, 2002), where deterministic sources and scatterers were recovered in a solar background. Yang (2018) recently studied PB holographic images in a homogeneous medium permeated by localized deterministic sources to study ghost images near the observational boundary.

In this paper we apply PB holography in a realistic helioseismological setting. First we rewrote the wave equation in Helmholtz form, in order to properly define the Green’s functions that are involved in the definition of PB holographic data. The background density and sound-speed are taken from a standard solar model. The model of wave excitation is described by a reasonable source covariance function, which leads to a solar-like power spectrum for acoustic oscillations.

The signal is defined as the expectation value of the holographic image intensity that results from perturbations in sound speed, density, and flows with respect to the reference solar model. The corresponding sensitivity kernels were computed in the first-order Born approximation (Gizon & Birch 2002; Birch & Gizon 2007; Braun et al. 2007; Birch et al. 2011). This signal must take into account the correlations between incident and scattered wave fields, which are both connected to the sources of excitation (turbulent convection).

Random noise in holographic images is due to the stochastic nature of the sources of excitation. While noise can sometimes be estimated from the data (Lindsey & Braun 1990; Braun & Birch 2008a), a theoretical understanding is useful to design holographic experiments. In this paper we have extended the noise model developed in time-distance helioseismology to holography (Gizon & Birch 2004; Fournier et al. 2014). We did not attempt to image individual sources as in Skartlien (2002), which in our view is not a well-posed problem (see also Lindsey et al. 2006), except in the case of sunquakes. Instead we considered sources to be specified through a source covariance function.

2. Reduced wave equation

At angular frequency ω and spatial position r in the computational domain V, the propagation of acoustic waves in a 3D heterogeneous moving medium is described by the displacement vector ξ(r, ω), solution to

( ω + i γ + i u · ) 2 ξ 1 ρ ( ρ c 2 · ξ ) + gravity terms = F , $$ \begin{aligned} - (\omega + \mathrm{i}\gamma + \mathrm{i}\boldsymbol{u}\cdot \nabla )^2\boldsymbol{\xi }- \frac{1}{\rho } \nabla \left( \rho c^2 \nabla \cdot \boldsymbol{\xi }\right) + \text{ gravity} \text{ terms} = \boldsymbol{F}, \end{aligned} $$(1)

where ρ(r) and c(r) are the density and sound speed, and u(r) is a steady vector flow. Wave attenuation is included through the function γ(r, ω). The source term F(r, ω) is a realization from a random process; it describes the stochastic excitation of the waves by turbulent convection. Following Lamb (1909) and Deubner & Gough (1984), we consider the scalar variable

ψ = ρ 1 / 2 c 2 · ξ $$ \begin{aligned} \psi = \rho ^{1/2} c^2 \nabla \cdot \boldsymbol{\xi }\end{aligned} $$(2)

to recast the wave equation into a Helmholtz-like equation

L ψ : = ( Δ + k 2 ) ψ 2 i ω ρ 1 / 2 c ρ u · ( ψ ρ 1 / 2 c ) = S , $$ \begin{aligned} L\psi := - (\Delta + k^2) \psi - \frac{2\mathrm{i}\omega }{\rho ^{1/2}c} \rho \boldsymbol{u}\cdot \nabla \left( \frac{ \psi }{ \rho ^{1/2}c} \right) = S, \end{aligned} $$(3)

where S = ρ1/2c2∇⋅F is a scalar source term. The local wavenumber k(r, ω) is given by

k 2 = ( ω 2 + 2 i ω γ ) ω c 2 c 2 , $$ \begin{aligned} k^2 = \frac{(\omega ^2 + 2\mathrm{i}\omega \gamma ) - \omega _{\rm c}^2}{c^2}, \end{aligned} $$(4)

where the squared acoustic cut-off frequency is

ω c 2 = ρ 1 / 2 c 2 Δ ( ρ 1 / 2 ) . $$ \begin{aligned} \omega _{\rm c}^2 = \rho ^{1/2} c^2 \Delta ( \rho ^{-1/2} ). \end{aligned} $$(5)

In obtaining Eq. (3), we ignored gravity terms and assumed slow variations of c, u, and γ compared to the wavelength (Gizon et al. 2017). The advection term is such that the corresponding operator is Hermitian symmetric for the inner product ψ 1 , ψ 2 = ψ 1 * ψ 2 dV $ \langle {\psi _1},{\psi _2}\rangle = \int {\psi _1^*} {\psi _2}{\mkern 1mu} {\rm{d}}V $ under the conditions that the flow conserves mass and that it does not cross the boundary (un = 0 on ∂V).

The stochastic sources of excitation are assumed to be stationary and spatially uncorrelated. They are described by a source covariance function of the form

E [ S ( r , ω ) S ( r , ω ) ] = M ( r , ω ) δ ( r r ) . $$ \begin{aligned} {\mathbb{E} }[S^*(\boldsymbol{r}, \omega ) S(\boldsymbol{r}^{\prime },\omega ) ] = M(\boldsymbol{r}, \omega ) \delta (\boldsymbol{r}-\boldsymbol{r}^{\prime }). \end{aligned} $$(6)

To solve Eq. (3), one needs to specify a boundary condition at the computational boundary. Here we apply an outgoing radiation boundary condition

n ψ = i k n ψ on V . $$ \begin{aligned} \partial _n \psi = \mathrm{i}k_n \psi \quad \text{ on} \partial V . \end{aligned} $$(7)

We applied the boundary condition (Atmo RBC 1) from Barucq et al. (2018), which assumes an exponential decay of the background density at the boundary of the domain but neglects curvature. Then, the local wavenumber kn from Eq. (7) is given by

k n 2 = ω 2 + 2 i ω γ c 2 1 4 H 2 , $$ \begin{aligned} k_n^2 = \frac{\omega ^2+2\mathrm{i}\omega \gamma }{c^2} - \frac{1}{4H^2}, \end{aligned} $$(8)

where H = −1/(d ln ρ/dr) is the density scale height at the boundary. The last term in Eq. (8) is connected to the cut-off frequency for an isothermal atmosphere (Lamb 1909), thus kn is an approximation of the wavenumber k from Eq. (4). Fournier et al. (2017) discuss several of the boundary conditions derived in Barucq et al. (2018).

3. Holographic image intensity

The following calculations are done at constant ω, thus we drop ω from the list of function arguments when not explicitly needed. The Porter–Bojarski integral is defined by Porter & Devaney (1982):

Φ α ( x , A ) : = A [ ψ ( r ) n H α ( r , x ) H α ( r , x ) n ψ ( r ) ] d r , $$ \begin{aligned} \Phi _\alpha (\boldsymbol{x}, A) := \int _{A} \left[\psi (\boldsymbol{r}) \partial _n H_\alpha (\boldsymbol{r}, \boldsymbol{x}) - H_\alpha (\boldsymbol{r}, \boldsymbol{x}) \partial _n \psi (\boldsymbol{r})\right] \mathrm{d}\boldsymbol{r}, \end{aligned} $$(9)

where Hα is an acoustic wave propagator for the reference medium and A is a surface on the Sun where ψ and ∂nψ are observed. The role of Hα is to propagate the wave field backward (or forward) in time, which leads to the concept of egression (or ingression) in LB holography (Lindsey & Braun 2000a).

Several choices for the propagators have been proposed in the literature, as detailed in Table 1. These depend on the outgoing (G0) and incoming ( G 0 $ G_0^ - $) Green’s functions defined in a reference medium with ρ0, c0, γ0, and u0 = 0:

Table 1.

Possible wave propagators.

L 0 [ G 0 ( r , r ) ] = δ ( r r ) and n G 0 = i k n G 0 on V , $$ \begin{aligned} L_0 [ G_0 (\boldsymbol{r}, \boldsymbol{r}^{\prime }) ] = \delta (\boldsymbol{r}- \boldsymbol{r}^{\prime })&\text{ and}&\partial _n G_0 = \mathrm{i}k_n G_0 \text{ on} \partial V,\end{aligned} $$(10)

L 0 [ G 0 ( r , r ) ] = δ ( r r ) and n G 0 = i k n G 0 on V , $$ \begin{aligned} L_0 [ G^-_0 (\boldsymbol{r}, \boldsymbol{r}^{\prime }) ] = \delta (\boldsymbol{r}- \boldsymbol{r}^{\prime })&\text{ and}&\partial _n G^-_0 = - \mathrm{i}k_n G^-_0 \text{ on} \partial V, \end{aligned} $$(11)

with

L 0 = ( Δ + k 0 2 ) , $$ \begin{aligned} L_0 = - (\Delta + k_0^2) , \end{aligned} $$(12)

where k0 is k in the reference medium and kn is from Eq. (8). The Green’s functions G 0 $ G_0^ * $ or G 0 $ G_0^ - $ are backward propagators (cf. egression), while ( ( G 0 ) $ {(G_0^ - )^ * } $ is a forward propagator (cf. ingression). When the surface A is closed, it is equivalent to use G 0 $ G_0^ * $ and Im G0 (Devaney & Porter 1985). Tsang et al. (1987) proposed H α = G 0 $ {H_\alpha } = G_0^ - $ as a backward propagator to correct for wave attenuation.

If the observations are made at the computational boundary and the wave field satisfies the same boundary condition as the Green’s function, then Eq. (9) reads

Φ α ( x , A ) = A ψ ( r ) [ n H α ( r , x ) i k n H α ( r , x ) ] d r . $$ \begin{aligned} \Phi _\alpha (\boldsymbol{x}, A) = \int _{A} \psi (\boldsymbol{r}) \left[\partial _n H_\alpha (\boldsymbol{r}, \boldsymbol{x}) - \mathrm{i}k_n H_\alpha (\boldsymbol{r}, \boldsymbol{x})\right] \mathrm{d}\boldsymbol{r}. \end{aligned} $$(13)

When H α = G 0 $ {H_\alpha } = G_0^ * $, we have

Φ ( x , A ) = 2 i Re [ k n ] A ψ ( r ) G 0 ( r , x ) d r , $$ \begin{aligned} \Phi (\boldsymbol{x}, A) = -2\mathrm{i}\text{ Re}[k_n] \int _{A} \psi (\boldsymbol{r}) G_0^*(\boldsymbol{r}, \boldsymbol{x}) \mathrm{d}\boldsymbol{r}, \end{aligned} $$(14)

which corresponds to the egression as defined by Lindsey & Braun (2000a). Thus LB and PB integrals are closely related, at least for the upper boundary condition employed here.

In LB holography, information is extracted from the egression-ingression correlation (wave-speed perturbations and flows) and from the egression power (source covariance). Analogously, we define the PB image intensity (or covariance) as

I α β ( x , A , A ) = Φ α ( x , A ) Φ β ( x , A ) . $$ \begin{aligned} I_{\alpha \beta }(\boldsymbol{x}, A, A^{\prime }) = \Phi ^*_\alpha (\boldsymbol{x}, A) \Phi _\beta (\boldsymbol{x}, A^{\prime }) . \end{aligned} $$(15)

For the case α = β we define

I α ( x , A ) = | Φ α ( x , A ) | 2 . $$ \begin{aligned} I_{\alpha }(\boldsymbol{x}, A) = |\Phi ^*_\alpha (\boldsymbol{x}, A) |^2 . \end{aligned} $$(16)

Different choices of pupils and propagators will provide sensitivity to different quantities as shown in Table 2. Scatterers are detected by correlating the forward and backward propagated wavefields. Different pupil shapes give access to different components of the flow (Table 2 and Fig. 1).

thumbnail Fig. 1.

Pupil geometries used to compute sound-speed or source kernels (P) and flow kernels (Qs), also see Table 2.

Table 2.

Proposed propagators and associated pupils (Hα, A) and (Hβ, A′) for different types of perturbations.

4. First-order perturbations with respect to a reference solar model

We wish to study how perturbations to the background medium affect the holographic images. Using the Born approximation, the first step is to express the perturbations to the wavefield and use this expression in Eq. (9) to obtain the perturbations to the PB integral and the image intensity.

4.1. Perturbations to the wavefield

We considered perturbations δc, δρ, δγ, u with respect to the reference medium defined by ρ0, c0, γ0, and u0 = 0. The perturbations to the sources of excitation are described through the perturbations to the source covariance,

M ( r , ω ) = M 0 ( r , ω ) + δ M ( r , ω ) , $$ \begin{aligned} M (\boldsymbol{r},\omega )= M_0 (\boldsymbol{r},\omega )+ \delta M(\boldsymbol{r},\omega ), \end{aligned} $$(17)

where

δ M ( r , ω ) = E [ S 0 ( r ) δ S ( r ) + δ S ( r ) S 0 ( r ) ] = ϵ ( r ) M 0 ( r , ω ) . $$ \begin{aligned} \delta M (\boldsymbol{r},\omega )= {\mathbb{E} }[S_0^*(\boldsymbol{r}) \delta S(\boldsymbol{r}^{\prime })+ \delta S^*(\boldsymbol{r}) S_0(\boldsymbol{r}^{\prime })] = \epsilon (\boldsymbol{r}) M_0 (\boldsymbol{r},\omega ). \end{aligned} $$(18)

Using the Born approximation up to first order, we write the wave field as

ψ ( r , ω ) = ψ 0 ( r , ω ) + δ ψ ( r , ω ) , $$ \begin{aligned} \psi (\boldsymbol{r}, \omega ) = \psi _0(\boldsymbol{r}, \omega ) + \delta \psi (\boldsymbol{r}, \omega ) , \end{aligned} $$(19)

where the zeroth- and first-order wave fields are given by

L 0 [ ψ 0 ] = S 0 , $$ \begin{aligned} L_0 [\psi _0]&= S_0 , \end{aligned} $$(20)

L 0 [ δ ψ ] = δ L [ ψ 0 ] + δ S . $$ \begin{aligned} L_0 [\delta \psi ]&= - \delta L [\psi _0] + \delta S . \end{aligned} $$(21)

The perturbed wave operator is

δ L [ ψ 0 ] = δ k 2 ψ 0 2 i ω ρ 0 1 / 2 c 0 ρ 0 u · ( ψ 0 ρ 0 1 / 2 c 0 ) , $$ \begin{aligned} \delta L [\psi _0] = - \delta k^2 \psi _0 - \frac{2\mathrm{i}\omega }{\rho ^{1/2}_0c_0} \rho _0\boldsymbol{u}\cdot \nabla \left( \frac{ \psi _0 }{ \rho ^{1/2}_0c_0} \right), \end{aligned} $$(22)

with

c 0 2 δ k 2 = 2 i ω δ γ ( ω 2 + 2 i ω γ 0 ) δ c 2 c 0 2 ( ω c 2 ρ ) 0 δ ρ . $$ \begin{aligned} c_0^2\, \delta k^2 = {2\mathrm{i}\omega \; \delta \gamma } - (\omega ^2+2\mathrm{i}\omega \gamma _0) \frac{\delta c^2}{c_0^2} - \left(\frac{\partial \omega _{\rm c}^2}{\partial \rho } \right)_0 \delta \rho . \end{aligned} $$(23)

In terms of the Green’s function G0, we have

ψ ( r ) = ψ 0 ( r ) + V G 0 ( r , r s ) δ k 2 ( r s ) ψ 0 ( r s ) d r s + V G 0 ( r , s ) δ S ( s ) d s + 2 i ω V ρ 0 ( r s ) u 0 ( r s ) · ( G 0 ( r , r s ) ρ 0 1 / 2 ( r s ) c 0 ( r s ) ) ψ 0 ( r s ) ρ 0 1 / 2 ( r s ) c 0 ( r s ) d r s . $$ \begin{aligned}&\psi (\boldsymbol{r}) = \psi _0(\boldsymbol{r}) + \int _V G_0(\boldsymbol{r}, \boldsymbol{r}_{\rm \! s}) \delta k^2(\boldsymbol{r}_{\rm \! s}) \psi _0(\boldsymbol{r}_{\rm \! s}) \mathrm{d}\boldsymbol{r}_{\rm \! s}+ \int _V G_0(\boldsymbol{r}, \boldsymbol{s}) \delta S (\boldsymbol{s}) \mathrm{d}\boldsymbol{s}\nonumber \\&\qquad \quad + 2\mathrm{i}\omega \int _V \rho _0(\boldsymbol{r}_{\rm \! s})\, \boldsymbol{u}_0(\boldsymbol{r}_{\rm \! s}) \cdot \nabla \left( \frac{G_0(\boldsymbol{r}, \boldsymbol{r}_{\rm \! s})}{\rho _0^{1/2}(\boldsymbol{r}_{\rm \! s}) c_0(\boldsymbol{r}_{\rm \! s})} \right) \frac{\psi _0(\boldsymbol{r}_{\rm \! s})}{\rho _0^{1/2}(\boldsymbol{r}_{\rm \! s}) c_0(\boldsymbol{r}_{\rm \! s})} \mathrm{d}\boldsymbol{r}_{\rm \! s}. \end{aligned} $$(24)

where the scattering location r​s spans the entire solar volume V.

4.2. Perturbations to the PB integral

For convenience, we introduce the source kernels

K α ( x , s , A ) = A [ G ( r , s ) n H α ( r , x ) H α ( r , x ) n G ( r , s ) ] d r , $$ \begin{aligned} K_\alpha (\boldsymbol{x},\boldsymbol{s}, A) = \int _{A} \left[G (\boldsymbol{r},\boldsymbol{s}) \partial _n H_\alpha (\boldsymbol{r}, \boldsymbol{x}) - H_\alpha (\boldsymbol{r}, \boldsymbol{x}) \partial _n G(\boldsymbol{r},\boldsymbol{s})\right] \mathrm{d}\boldsymbol{r}, \end{aligned} $$(25)

such that the PB integral is given by

Φ α ( x , A ) = V K α ( x , s , A ) S ( s ) d s . $$ \begin{aligned} \Phi _\alpha (\boldsymbol{x}, A) = \int _V K_\alpha (\boldsymbol{x},\boldsymbol{s}, A) S(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}. \end{aligned} $$(26)

We denote by Kα, 0 the source kernel in the reference medium (when G is replaced by G0).

Replacing the wavefield by its first order expansion (Eq. (24)) in the definition of the PB integral (Eq. (9)), one obtains

Φ α ( x , A ) = Φ α , 0 ( x , A ) + δ Φ α ( x , A ) , $$ \begin{aligned} \Phi _\alpha (\boldsymbol{x}, A) = \Phi _{\alpha ,0} (\boldsymbol{x},A) + \delta \Phi _\alpha (\boldsymbol{x},A), \end{aligned} $$(27)

where Φα, 0 is the value for the reference medium and δΦα is due to perturbations in the medium:

δ Φ α ( x , A ) = V δ K α ( x , s , A ) S ( s ) d s + V K α , 0 ( x , s , A ) δ S ( s ) d s , $$ \begin{aligned} \delta \Phi _{\alpha } (\boldsymbol{x},A) = \int _V \delta K_\alpha (\boldsymbol{x},\boldsymbol{s}, A) S(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}+ \int _V K_{\alpha ,0} (\boldsymbol{x},\boldsymbol{s}, A) \delta S(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}, \end{aligned} $$(28)

where

δ K α ( x , s , A ) = V K α , 0 ( x , r s , A ) δ k 2 ( r s ) G 0 ( r s , s ) d r s + 2 i ω V K α , 0 ( x , r s , A ) ρ 0 1 / 2 ( r s ) c 0 ( r s ) ρ 0 u 0 · ( G 0 ( r s , s ) ρ 0 1 / 2 ( r s ) c 0 ( r s ) ) d r s . $$ \begin{aligned} \delta K_\alpha&(\boldsymbol{x},\boldsymbol{s}, A) = \int _V K_{\alpha ,0}(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A) \delta k^2(\boldsymbol{r}_{\rm \! s}) G_0(\boldsymbol{r}_{\rm \! s},\boldsymbol{s}) \mathrm{d}\boldsymbol{r}_{\rm \! s}\nonumber \\&\qquad \qquad + 2\mathrm{i}\omega \int _V \frac{K_{\alpha ,0}(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A)}{\rho _0^{1/2}(\boldsymbol{r}_{\rm \! s}) c_0(\boldsymbol{r}_{\rm \! s})} \rho _0 \boldsymbol{u}_0 \cdot \nabla \left( \frac{G_0(\boldsymbol{r}_{\rm \! s},\boldsymbol{s})}{\rho _0^{1/2}(\boldsymbol{r}_{\rm \! s}) c_0(\boldsymbol{r}_{\rm \! s})} \right) \mathrm{d}\boldsymbol{r}_{\rm \! s}. \end{aligned} $$(29)

4.3. Perturbations to the image intensity

We write the holographic image intensity to first order in the form

I α β ( x , A , A ) = I α β , 0 ( x , A , A ) + δ I α β ( x , A , A ) . $$ \begin{aligned} I_{\alpha \beta }(\boldsymbol{x}, A, A^{\prime }) = I_{\alpha \beta , 0}(\boldsymbol{x}, A, A^{\prime }) + {\delta {I}}_{\alpha \beta } (\boldsymbol{x}, A, A^{\prime }) . \end{aligned} $$(30)

The expectation values of the zeroth- and first-order image intensities are

E [ I α β , 0 ( x , A , A ) ] = V K α , 0 ( x , s , A ) K β , 0 ( x , s , A ) M 0 ( s ) d s , $$ \begin{aligned} {\mathbb{E} }[I_{\alpha \beta , 0}(\boldsymbol{x}, A, A^{\prime })]&= \int _{V}K^{*}_{\alpha ,0}(\boldsymbol{x},\boldsymbol{s}, A) K_{\beta ,0}(\boldsymbol{x},\boldsymbol{s}, A^{\prime }) M_{0}(\boldsymbol{s}) \mathrm{d}\boldsymbol{s},\end{aligned} $$(31)

E [ δ I α β ( x , A , A ) ] = V K α , 0 ( x , s , A ) δ K β ( x , s , A ) M 0 ( s ) d s + V δ K α ( x , s , A ) K β , 0 ( x , s , A ) M 0 ( s ) d s + V K α , 0 ( x , s , A ) K β , 0 ( x , s , A ) ϵ ( s ) M 0 ( s ) d s . $$ \begin{aligned} {\mathbb{E} }[{\delta {I}}_{\alpha \beta }(\boldsymbol{x}, A, A^{\prime }) ]&= \int _{V} K_{\alpha ,0}^{*}(\boldsymbol{x},\boldsymbol{s}, A) {\delta {K}}_{\beta }(\boldsymbol{x},\boldsymbol{s}, A^{\prime }) M_{0}(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}\nonumber \\&\qquad + \int _{V} {\delta {K}}_{\alpha }^*(\boldsymbol{x},\boldsymbol{s}, A) K_{\beta ,0}(\boldsymbol{x},\boldsymbol{s}, A^{\prime }) M_{0}(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}\nonumber \\&\qquad + \int _V K_{\alpha ,0}^*(\boldsymbol{x},\boldsymbol{s}, A) K_{\beta ,0}(\boldsymbol{x},\boldsymbol{s}, A^{\prime }) \epsilon (\boldsymbol{s}) M_0(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}. \end{aligned} $$(32)

Using the definition of the source kernels, we obtain

E [ I α β , 0 ( x , A , A ) ] = C 0 ( r , r ) α β ( x , A , A ) , $$ \begin{aligned} {\mathbb{E} }[I_{\alpha \beta ,0} (\boldsymbol{x}, A, A^{\prime })] = \langle \!\langle C_0(\boldsymbol{r},\boldsymbol{r}^{\prime }) \rangle \!\rangle _{\alpha \beta }(\boldsymbol{x}, A, A^{\prime }) , \end{aligned} $$(33)

where

C 0 ( r , r ) = V G 0 ( r , s ) G 0 ( r , s ) M 0 ( s ) d s $$ \begin{aligned} C_0(\boldsymbol{r},\boldsymbol{r}^{\prime }) = \int _V G_0^*(\boldsymbol{r},\boldsymbol{s}) G_0(\boldsymbol{r}^{\prime },\boldsymbol{s}) M_0(\boldsymbol{s}) \mathrm{d}\boldsymbol{s} \end{aligned} $$(34)

is the expectation value of the cross-covariance function and, for any function F(r, r′), the double brackets mean

F ( r , r ) α β ( x , A , A ) = A d r A d r [ n H α ( r , x ) F ( r , r ) n H β ( r , x ) H α ( r , x ) n F ( r , r ) n H β ( r , x ) n H α ( r , x ) n F ( r , r ) H β ( r , x ) + H α ( r , x ) n n F ( r , r ) H β ( r , x ) ] . $$ \begin{aligned} \langle \!\langle F(\boldsymbol{r},\boldsymbol{r}^{\prime }) \rangle \!\rangle _{\alpha \beta }(\boldsymbol{x}, A,A^{\prime })&= \nonumber \\ \int _A \mathrm{d}\boldsymbol{r}\int _{A^{\prime }} \mathrm{d}\boldsymbol{r}^{\prime }\; \biggl [&\partial _n H^*_\alpha (\boldsymbol{r},\boldsymbol{x}) F(\boldsymbol{r},\boldsymbol{r}^{\prime }) \partial _{n^{\prime }} H_\beta (\boldsymbol{r}^{\prime },\boldsymbol{x}) \nonumber \\&- H^*_\alpha (\boldsymbol{r},\boldsymbol{x}) \partial _n F(\boldsymbol{r},\boldsymbol{r}^{\prime }) \partial _{n^{\prime }}H_\beta (\boldsymbol{r}^{\prime },\boldsymbol{x}) \nonumber \\&- \partial _{n}H^*_\alpha (\boldsymbol{r},\boldsymbol{x}) \partial _n F(\boldsymbol{r},\boldsymbol{r}^{\prime }) H_\beta (\boldsymbol{r}^{\prime },\boldsymbol{x}) \nonumber \\&+ H^*_\alpha (\boldsymbol{r},\boldsymbol{x}) \, \partial _{n} \partial _{n^{\prime }}F(\boldsymbol{r},\boldsymbol{r}^{\prime }) \, H_\beta (\boldsymbol{r}^{\prime },\boldsymbol{x}) \biggr ] . \end{aligned} $$(35)

The perturbation to the image intensity is given by

E [ δ I α β ( x , A , A ) ] = V ϵ ( s ) K α β ϵ ( x , s , A , A ) d s + V ( δ k 2 ( r s ) K α β k ( x , r s , A , A ) + δ k 2 ( r s ) K β α k ( x , r s , A , A ) ) d r s + 2 i ω V u ( r s ) · ( K α β u ( x , r s , A , A ) K β α u ( x , r s , A , A ) ) d r s , $$ \begin{aligned}&{\mathbb{E} }[{\delta {I}}_{\alpha \beta } (\boldsymbol{x}, A, A^{\prime })]\nonumber \\&\qquad = \int _V \epsilon (\boldsymbol{s}) \mathcal{K} _{\alpha \beta }^\epsilon ( \boldsymbol{x},\boldsymbol{s}, A, A^{\prime }) \, \mathrm{d}\boldsymbol{s}\nonumber \\&\qquad \quad + \int _V \left( \delta k^{2*}(\boldsymbol{r}_{\rm \! s}) \, \mathcal{K} ^k_{\alpha \beta }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A, A^{\prime }) + \delta k^{2}(\boldsymbol{r}_{\rm \! s}) \, \mathcal{K} ^{k*}_{\beta \alpha }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A^{\prime }, A) \right) \, \mathrm{d}\boldsymbol{r}_{\rm \! s}\nonumber \\&\qquad \quad + 2\mathrm{i}\omega \int _V \boldsymbol{u}(\boldsymbol{r}_{\rm \! s}) \cdot \left( \boldsymbol{ \mathcal{K} }^u_{\alpha \beta }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A, A^{\prime }) - \boldsymbol{\mathcal{K} }^{u *}_{\beta \alpha }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A^{\prime }, A) \right) \, \mathrm{d}\boldsymbol{r}_{\rm \! s}, \end{aligned} $$(36)

where

K α β ϵ ( x , s , A , A ) = M 0 ( s ) G 0 ( r , s ) G 0 ( r , s ) α β ( x , A , A ) , $$ \begin{aligned} \mathcal{K} _{\alpha \beta }^\epsilon&(\boldsymbol{x},\boldsymbol{s}, A, A^{\prime }) =M_0(\boldsymbol{s}) \left\langle \!\left\langle G_0^*(\boldsymbol{r},\boldsymbol{s}) G_0(\boldsymbol{r}^{\prime },\boldsymbol{s}) \right\rangle \!\right\rangle _{\alpha \beta }(\boldsymbol{x},A,A^{\prime }), \end{aligned} $$(37)

K α β k ( x , r s , A , A ) = C 0 ( r s , r ) G 0 ( r , r s ) α β ( x , A , A ) , $$ \begin{aligned} \mathcal{K} ^k_{\alpha \beta }&(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A, A^{\prime }) = \left\langle \!\left\langle C_0(\boldsymbol{r}_{\rm \! s},\boldsymbol{r}^{\prime }) G_0^*(\boldsymbol{r}, \boldsymbol{r}_{\rm \! s}) \right\rangle \!\right\rangle _{\alpha \beta }(\boldsymbol{x},A,A^{\prime }), \end{aligned} $$(38)

K α β u ( x , r s , A , A ) = ( C 0 ( r s , r ) ρ 0 ( r s ) 1 / 2 c 0 ( r s ) ) ρ 0 ( r s ) 1 / 2 G 0 ( r , r s ) c 0 ( r s ) α β ( x , A , A ) . $$ \begin{aligned} \boldsymbol{\mathcal{K} }_{\alpha \beta }^{u}&(\boldsymbol{x}, \boldsymbol{r}_{\rm \! s},A, A^{\prime }) = \nonumber \\&\left\langle \!\left\langle \nabla \left( \frac{C_0(\boldsymbol{r}_{\rm \! s},\boldsymbol{r}^{\prime })}{\rho _0(\boldsymbol{r}_{\rm \! s})^{1/2} c_0(\boldsymbol{r}_{\rm \! s})} \right) \frac{\rho _0(\boldsymbol{r}_{\rm \! s})^{1/2} G_0^*(\boldsymbol{r}, \boldsymbol{r}_{\rm \! s}) }{c_0(\boldsymbol{r}_{\rm \! s})} \right\rangle \!\right\rangle _{\alpha \beta }(\boldsymbol{x},A,A^{\prime }). \end{aligned} $$(39)

The kernels for δk2 and δk2* can be combined using Eq. (23) to obtain kernels for sound-speed K α β c $ {\cal{K}}_{\alpha \beta}^{c} $, density K α β ρ $ {\cal{K}}_{\alpha\beta}^\rho $ and attenuation K α β γ $ {\cal{K}}_{\alpha\beta}^\gamma $. For example,

E [ δ I α β ( x , A , A ) ] = V δ c ( r s ) K α β c ( x , r s , A , A ) d r s $$ \begin{aligned} {\mathbb{E} }[{\delta {I}}_{\alpha \beta } (\boldsymbol{x}, A, A^{\prime })] = \int _V \delta c (\boldsymbol{r}_{\rm \! s}) \, \mathcal{K} _{\alpha \beta }^c(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}, A, A^{\prime }) \, \mathrm{d}\boldsymbol{r}_{\rm \! s} \end{aligned} $$(40)

with

K α β c ( x , r s , A , A ) = 2 ( ω 2 2 i ω γ ) c 0 3 ( r s ) K α β k ( x , r s , A , A ) 2 ( ω 2 + 2 i ω γ ) c 0 3 ( r s ) K β α k ( x , r s , A , A ) . $$ \begin{aligned} \mathcal{K} ^c_{\alpha \beta } (\boldsymbol{x},\boldsymbol{r}_{\rm \! s},A, A^{\prime })=&- \frac{2 (\omega ^2-2\mathrm{i}\omega \gamma )}{c_0^3(\boldsymbol{r}_{\rm \! s})} \mathcal{K} ^{k}_{\alpha \beta }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s},A, A^{\prime }) \nonumber \\&-\frac{2 (\omega ^2+2\mathrm{i}\omega \gamma )}{c_0^3(\boldsymbol{r}_{\rm \! s})} \mathcal{K} ^{k*}_{\beta \alpha }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s},A^{\prime }, A) . \end{aligned} $$(41)

4.4. Choice of the source covariance

In order to compute the above kernels, we still need to choose the source covariance function M0 in order to define the reference cross-covariance C0 using Eq. (34). One possibility is to place the sources at a single depth, a few hundred kilometers below the solar surface.

Another possibility is to choose a source covariance of the form

M 0 ( r , ω ) = Π ( ω ) γ ( r , ω ) c 0 2 ( r ) , $$ \begin{aligned} M_0(\boldsymbol{r},\omega ) = \Pi (\omega ) \frac{\gamma (\boldsymbol{r}, \omega )}{c_0^2(\boldsymbol{r})}, \end{aligned} $$(42)

where Π(ω) is the source power spectrum (see Gizon et al. 2017). This choice implies

C 0 ( r , r , ω ) = Π ( ω ) 2 ω Im G 0 ( r , r , ω ) + surface term . $$ \begin{aligned} C_0(\boldsymbol{r},\boldsymbol{r}^{\prime },\omega ) = \frac{\Pi (\omega )}{2\omega } \mathrm{Im}\, G_0 (\boldsymbol{r}^{\prime }, \boldsymbol{r}, \omega ) + \text{ surface} \text{ term}. \end{aligned} $$(43)

The surface term depends on the boundary condition. It vanishes for a Dirichlet boundary condition (free surface), while it remains for radiative boundary conditions (e.g., Sommerfeld). Below the acoustic cutoff frequency, modes are trapped well below the observational and computational boundaries and the surface term is negligible. In this paper we use Eq. (42) in the convection zone and switch off the sources above the photosphere. By doing so, the surface term in Eq. (43) vanishes.

5. Noise

To compute the noise level, we computed the variance of the image intensity in the absence of scatterers:

Var [ I α β , 0 ( x ) ] = Var [ K α , 0 ( x , s ) K β , 0 ( x , s ) S ( s ) S ( s ) d s d s ] = V 4 K α , 0 ( x , s 1 ) K β , 0 ( x , s 1 ) K α , 0 ( x , s 2 ) K β , 0 ( x , s 2 ) × M 4 ( s 1 , s 1 , s 2 , s 2 ) d s 1 d s 1 d s 2 d s 2 | V 2 K α , 0 ( x , s ) K β , 0 ( x , s ) E [ S ( s ) S ( s ) ] d s d s | 2 , $$ \begin{aligned} \mathrm{Var}[ I_{\alpha \beta ,0}(\boldsymbol{x}) ]&= \mathrm{Var}\left[ \int K^*_{\alpha ,0}(\boldsymbol{x}, \boldsymbol{s}) K_{\beta ,0}(\boldsymbol{x}, \boldsymbol{s}^{\prime }) S^*(\boldsymbol{s}) S(\boldsymbol{s}^{\prime }) \mathrm{d}\boldsymbol{s}\mathrm{d}\boldsymbol{s}^{\prime }\right] \nonumber \\&= \int _{V^4} K^*_{\alpha ,0}(\boldsymbol{x}, \boldsymbol{s}_1) K_{\beta ,0}(\boldsymbol{x}, \boldsymbol{s}_1^{\prime }) K_{\alpha ,0}(\boldsymbol{x}, \boldsymbol{s}_2) K^*_{\beta ,0}(\boldsymbol{x}, \boldsymbol{s}_2^{\prime }) \nonumber \\&\quad \times M_4(\boldsymbol{s}_1, \boldsymbol{s}_1^{\prime }, \boldsymbol{s}_2, \boldsymbol{s}_2^{\prime }) \, \mathrm{d}\boldsymbol{s}_1 \mathrm{d}\boldsymbol{s}_1^{\prime }\mathrm{d}\boldsymbol{s}_2 \mathrm{d}\boldsymbol{s}_2^{\prime }\nonumber \\&\qquad - \left| \int _{V^2} K^*_{\alpha ,0}(\boldsymbol{x}, \boldsymbol{s}) K_{\beta ,0}(\boldsymbol{x}, \boldsymbol{s}^{\prime }) {\mathbb{E} } \left[S^*(\boldsymbol{s}) S(\boldsymbol{s}^{\prime }) \right] \mathrm{d}\boldsymbol{s}\mathrm{d}\boldsymbol{s}^{\prime }\right|^2, \end{aligned} $$(44)

where

M 4 ( s 1 , s 1 , s 2 , s 2 ) = E [ S ( s 1 ) S ( s 1 ) S ( s 2 ) S ( s 2 ) ] . $$ \begin{aligned} M_4(\boldsymbol{s}_1, \boldsymbol{s}_1^{\prime }, \boldsymbol{s}_2, \boldsymbol{s}_2^{\prime }) = {\mathbb{E} }\left[ S^*(\boldsymbol{s}_1)S(\boldsymbol{s}_1^{\prime })S(\boldsymbol{s}_2)S^*(\boldsymbol{s}_2^{\prime }) \right]. \end{aligned} $$(45)

Under the (very reasonable) assumption that S is a realization drawn from a Gaussian random process, the fourth-order moment is the sum of products of the second-order moments:

M 4 ( s 1 , s 1 , s 2 , s 2 ) = E [ S ( s 1 ) S ( s 1 ) ] E [ S ( s 2 ) S ( s 2 ) ] + E [ S ( s 1 ) S ( s 2 ) ] E [ S ( s 1 ) S ( s 2 ) ] + E [ S ( s 1 ) S ( s 2 ) ] E [ S ( s 2 ) S ( s 1 ) ] . $$ \begin{aligned}&M_4(\boldsymbol{s}_1, \boldsymbol{s}_1^{\prime }, \boldsymbol{s}_2, \boldsymbol{s}_2^{\prime }) = {\mathbb{E} } \left[S^*(\boldsymbol{s}_1) S(\boldsymbol{s}_1^{\prime }) \right] {\mathbb{E} } \left[S(\boldsymbol{s}_2) S^*(\boldsymbol{s}_2^{\prime }) \right] \nonumber \\&\qquad \qquad \qquad \qquad + {\mathbb{E} } \left[S^*(\boldsymbol{s}_1) S(\boldsymbol{s}_2) \right] {\mathbb{E} }\left[S(\boldsymbol{s}_1^{\prime }) S^*(\boldsymbol{s}_2^{\prime }) \right] \nonumber \\ &\qquad \qquad \qquad \qquad + {\mathbb{E} } \left[S^*(\boldsymbol{s}_1) S^*(\boldsymbol{s}_2^{\prime }) \right] {\mathbb{E} }\left[S(\boldsymbol{s}_2) S(\boldsymbol{s}_1^{\prime })\right]. \end{aligned} $$(46)

The first term in M4 cancels out the squared term in Eq. (44). The third term is zero as the frequencies are uncorrelated. Thus,

Var [ I α β , 0 ( x ) ] = V 4 K α , 0 ( x , s 1 ) K β , 0 ( x , s 1 ) K α , 0 ( x , s 2 ) K β , 0 ( x , s 2 ) × E [ S ( s 1 ) S ( s 2 ) ] E [ S ( s 1 ) S ( s 2 ) ] d s 1 d s 1 d s 2 d s 2 = V | K α , 0 ( x , s ) | 2 M ( s ) d s V | K β , 0 ( x , s ) | 2 M ( s ) d s = E [ I α , 0 ( x ) ] E [ I β , 0 ( x ) ] . $$ \begin{aligned} \mathrm{Var}[ I_{\alpha \beta ,0}(\boldsymbol{x}) ]&= \int _{V^4} K^*_{\alpha ,0}(\boldsymbol{x}, \boldsymbol{s}_1) K_{\beta ,0}(\boldsymbol{x}, \boldsymbol{s}_1^{\prime }) K_{\alpha ,0}(\boldsymbol{x}, \boldsymbol{s}_2) K^*_{\beta ,0}(\boldsymbol{x}, \boldsymbol{s}_2^{\prime }) \nonumber \\&\quad \times {\mathbb{E} } \left[S^*(\boldsymbol{s}_1) S(\boldsymbol{s}_2) \right] {\mathbb{E} } \left[S(\boldsymbol{s}_1^{\prime }) S^*(\boldsymbol{s}_2^{\prime }) \right] \mathrm{d}\boldsymbol{s}_1 \mathrm{d}\boldsymbol{s}_1^{\prime }\mathrm{d}\boldsymbol{s}_2 \mathrm{d}\boldsymbol{s}_2^{\prime }\nonumber \\&= \int _V |K_{\alpha ,0}(\boldsymbol{x}, {\boldsymbol{s}}) |^2 M({\boldsymbol{s}}) \mathrm{d}{\boldsymbol{s}} \int _V |K_{\beta ,0}(\boldsymbol{x}, {\boldsymbol{s}}) |^2 M({\boldsymbol{s}}) \mathrm{d}{\boldsymbol{s}} \nonumber \\&= {\mathbb{E} }[ I_{\alpha ,0}(\boldsymbol{x}) ] {\mathbb{E} }[ I_{\beta ,0}(\boldsymbol{x}) ] . \end{aligned} $$(47)

When α = β, the standard deviation of Iα, 0 is equal to its expectation value. This is because the probability density function of Iα, 0 is a χ2 with two degrees of freedom, in other words, an exponential distribution.

6. Average over frequencies

In order to increase the signal-to-noise ratio (S/N), one usually averages the image intensity over a range of frequencies [ω0 − Δω/2, ω0 + Δω/2]. For an observation duration T, this interval will contain N = ΔωT/2π independent frequencies.

The frequency-averaged perturbation to the image intensity (i.e. the signal) is denoted by

δ I α β ( x ) = 1 N i = 1 N δ I α β ( x , ω i ) . $$ \begin{aligned} \langle {\delta {I}}_{\alpha \beta }(\boldsymbol{x})\rangle = \frac{1}{N} \sum _{i=1}^N {\delta {I}}_{\alpha \beta }(\boldsymbol{x}, \omega _i). \end{aligned} $$(48)

The variance of the noise in the average image intensity is then given by

Var I α β , 0 ( x ) = Var ( 1 N i = 1 N I α β , 0 ( x , ω i ) ) $$ \begin{aligned} \mathrm{Var} \left\langle I_{\alpha \beta , 0}(\boldsymbol{x}) \right\rangle&= \mathrm{Var} \left( \frac{1}{N} \sum _{i=1}^N I_{\alpha \beta ,0} (\boldsymbol{x},\omega _i) \right) \end{aligned} $$(49)

= 1 N 2 i = 1 N Var I α β , 0 ( x , ω i ) $$ \begin{aligned}&= \frac{1}{N^2} \sum _{i=1}^N \mathrm{Var}\ I_{\alpha \beta ,0} (\boldsymbol{x},\omega _i) \end{aligned} $$(50)

= 1 N Var I α β , 0 ( x ) , $$ \begin{aligned}&= \frac{1}{N} \left\langle \mathrm{Var} \; I_{\alpha \beta ,0} (\boldsymbol{x}) \right\rangle , \end{aligned} $$(51)

since the noise in the data at different frequencies is uncorrelated.

The expected S/N is thus

SNR ( x ) = | E δ I α β ( x ) | Var I α β , 0 ( x ) = N | E δ I α β ( x ) | E [ I α , 0 ( x ) ] E [ I β , 0 ( x ) ] · $$ \begin{aligned} \mathrm{SNR}(\boldsymbol{x}) = \frac{ \left| {\mathbb{E} } \langle {\delta {I}}_{\alpha \beta }(\boldsymbol{x}) \rangle \right|}{\sqrt{ \mathrm{Var} \langle I_{\alpha \beta , 0} (\boldsymbol{x}) \rangle } } = \frac{ \sqrt{N} \, \left| {\mathbb{E} } \langle {\delta {I}}_{\alpha \beta }(\boldsymbol{x}) \rangle \right|}{\sqrt{ \langle {\mathbb{E} }[I_{\alpha ,0}(\boldsymbol{x})] \,{\mathbb{E} }[I_{\beta ,0}(\boldsymbol{x})]\rangle } } \cdot \end{aligned} $$(52)

The number of available frequencies N within a fixed frequency band Δω is proportional to the observation duration T, hence the noise level goes like T−1/2. Provided that the frequency interval Δω is small with respect to the variations of the signal, then the S/N will increase as T1/2.

7. Example computations

In order to illustrate the theory, we have computed holographic images in the presence of sound-speed perturbations at different depths and calculate the corresponding S/Ns.

7.1. Reference Green’s function

The main input quantity required to compute PB integrals is the reference Green’s function (Eq. (10)). Here it is computed in the frequency domain using the spherically-symmetric standard solar Model S (Christensen-Dalsgaard et al. 1996). The wave attenuation model is taken from Gizon et al. (2017). Below 5.3 mHz, we have γ = γ0|ω/ω0|5.77, where γ0/2π = 4.29 μHz and ω0/2π = 3 mHz. Above 5.3 mHz, γ/2π = 125 μHz is kept constant. The radiation boundary condition defined by Eq. (7) is applied at the computational boundary located 500 km above the solar surface with the local wavenumber kn (where H = 105 km). The wave equation is solved using the finite-element solver Montjoie (Duruflé 2006; Gizon et al. 2017).

The reference Green’s function only depends on the angular distance Θ between the two points at radii r and r′. To speed up the computations, we place one of the points on the polar axis and compute the axisymmetric component of the Green’s function G l m=0 (r,r,ω) $ G_l^{m = 0}(r,r\prime ,\omega ) $ at each spherical harmonic degree l, to obtain:

G 0 ( r , r , ω ) l = 0 l max G l m = 0 ( r , r , ω ) P l ( cos Θ ) , $$ \begin{aligned} G_0(\boldsymbol{r},\boldsymbol{r}^{\prime }, \omega ) \simeq \sum _{l=0}^{l_{\rm max}} G_l^{m=0}(r,r^{\prime }, \omega ) P_l(\cos \Theta ) , \end{aligned} $$(53)

where we used an approximate equality because the sum is truncated at lmax = 300.

7.2. Sound-speed kernels

The sound-speed kernel is computed using Eq. (41) and the definition of K αβ k $ {\cal K}_{\alpha \beta }^k $. We need to evaluate two surface integrals, which can be computed analytically via a decomposition of all quantities into spherical harmonic coefficients (Fournier et al. 2018).

Figure 2 shows a sound-speed kernel K αβ c $ {\cal K}_{\alpha \beta }^c $ at a single frequency of 3 mHz. The pupil P is a polar cap of angular size 120° and the wave propagators Hα and Hβ are given in Table 2. As expected the kernel peaks around the scatterer position at z = 0.7 R. A cut along the polar axis is shown in Fig. 3; the kernel width is about half the local wavelength. In addition, we find ghost values at the antipode.

thumbnail Fig. 2.

Meridional slice through the sound-speed kernel K αβ c $ {\cal K}_{\alpha \beta }^c $(x, r​s) computed in Model S at a single frequency of 3 mHz, in units of 10−34 kg m−3 s−3. Both the real (left panel) and imaginary (right panel) parts of the kernel are shown. The scatterer at zs = 0.7 R is indicated by the crosses. The observation pupil P is a polar cap of full angular size 120°. We notice the “ghost values” at the antipode.

thumbnail Fig. 3.

Cut along the z axis through the sound-speed kernel from Fig. 2. The scatterer is at zs = 0.7 R.

7.3. Signal

At position r s = z s z ̂ $ {\boldsymbol{r}_{\mathrm{\! s}}} = {z_{\mathrm{s}}} {\hat{\mathbf{z}}} $ along the polar axis, we consider a localized increase in sound speed of 10% over a volume Vs, such that the signal (Eq. (40)) may be written as

E [ δ I α β ( x ) ] 0.1 V s c 0 ( r s ) K α β c ( x , r s ) . $$ \begin{aligned} {\mathbb{E} }[{\delta {I}}_{\alpha \beta }(\boldsymbol{x})] \simeq 0.1 V_{\rm s} \, c_0(\boldsymbol{r}_{\rm \! s}) \,\mathcal{K} ^c_{\alpha \beta }(\boldsymbol{x},\boldsymbol{r}_{\rm \! s}). \end{aligned} $$(54)

The volume Vs is that of a ball of diameter λ(rs)/2 = π/[Re k(rs, ω0)] with ω0/2π = 3 mHz. This is an approximate but much faster way to compute the effect of a perturbation of volume comparable to the highest possible holographic resolution. We checked that the answer does not differ significantly from the one obtained by integrating numerically the kernel over the ball of volume Vs. For reference, we note that λ/2 = 38 Mm at r = 0.7 R and λ/2 = 20 Mm at r = 0.9 R.

Figure 4 shows the signal |𝔼[δIαβ(x)]| for sound-speed perturbations located at two different depths zs = 0.7 R (red curve) and 0.9 R (blue curve). The pupil P and the wave propagators are the same as those of Fig. 2. The left panel of Fig. 4 shows the results at a single frequency of 3 mHz. With only one frequency, the signal peaks close to the scattering location but the spatial resolution is relatively poor, with a ghost on the far side. To demonstrate the benefits of averaging, the right panel shows the signal after averaging over 101 frequencies uniformly distributed in the interval 2.75–3.25 mHz. The frequency resolution 5 μHz corresponds to an observation duration T = 55.5 h. Averaging over frequencies improves the spatial resolution which approaches λ/2 and the ghost is suppressed. A horizontal line segment is plotted on the right panel at each depth to mark half the local wavelength.

thumbnail Fig. 4.

Left panel: image intensity |𝔼[δIαβ(x)]| at a single frequency of 3 mHz, as a function of position x along the z-axis. The sound-speed perturbation (see Eq. (54)) is placed at two different positions along the z-axis, zs = 0.7 R (red) and 0.9 R (blue). The standard deviation of the noise Var I α β , 0 ( x ) $ \sqrt{ \mathrm{Var} \langle I_{\alpha\beta, 0} (\boldsymbol{x}) \rangle } $ is given by the black curve. We note that the jagged aspect of the curves is not due to numerical inaccuracies. Right panel: Image intensity and noise level after averaging over 101 frequencies uniformly distributed in the interval from 2.75 to 3.25 mHz. The frequency resolution is 5 μHz, implying an observation duration of T = 55.5 h. A horizontal line segment is plotted at each depth to mark half of the local wavelength.

As seen in Fig. 5 the spatial extent of the frequency-averaged kernels is approximately λ/2 in both the radial and horizontal directions, for all scattering points in the range 0.6 <  zs/R <  0.98. Thus helioseismic holography is a diffraction-limited imaging technique as suggested by Lindsey & Braun (1997).

thumbnail Fig. 5.

Radial and horizontal widths of the frequency-averaged sound-speed kernel |⟨𝒦c⟩| as functions of scattering position. These are close to half the local wavelength at 3 mHz (dotted line).

7.4. Noise

The noise is obtained from Eq. (47), which requires the computation of 𝔼[Iα, 0] and 𝔼[Iβ, 0] using Eq. (33). The reference cross-covariance C0 is precomputed. The double surface integral is evaluated in a similar way as for the kernel computations.

For a frequency of 3 mHz the left panel of Fig. 4 (black curve) shows the noise level, together with the signal described in the previous section. The jagged aspect of the noise variations with position is not due to numerical inaccuracies but to the details of the Green’s function. As shown in the right panel of Fig. 4, the noise level goes down by a factor of about ten after averaging over 101 frequencies, and varies more smoothly with depth.

Braun & Birch (2008a) studied the noise level in observed travel times measured from LB holography. These measurements, however, include contributions from supergranulation and so are not directly comparable to what is shown in Fig. 4.

7.5. Signal-to-noise ratio

Figure 6 shows the S/N as a function of scatterer location. We recall that the sound-speed perturbation is specified by Eq. (54) and is the same as in Sect. 7.3. The results are shown at a single frequency of 3 mHz and after averaging over 101 frequencies in the interval 2.75–3.25 mHz. After averaging, the S/N is above two and is roughly independent of depth in the range 0.6 <  zs/R <  0.98 for a pupil of angular size 120°. We note that the ghost at −zs is much below the noise level.

thumbnail Fig. 6.

S/N in PB image intensity for a 10% sound-speed perturbation over a volume Vs(zs) placed at zs along the polar axis (Eq. (54)). The results are shown at a single frequency of 3 mHz (solid curve) and after averaging over 101 frequencies in the interval from 2.75 to 3.25 mHz (dashed curve).

We find that both signal and noise vary rapidly with frequency for deep located sound-speed scatterers. Figure 7 shows an example of a sound-speed scatterer located at zs = 0.7 R. Strong frequency variations in signal and noise are evident for frequencies below 3.5 mHz. This can be understood as follows. Low-frequency modes have narrowly-peaked power spectra due to their long lifetimes. At these low frequencies, the amplitude of the kernel function and the noise will change rapidly when the frequency coincides with a particular mode frequency. Additionally, the kernel function may not peak exactly at the sound-speed scatterer position when only a few modes contribute to the kernel function. Figure 8 shows the S/N as a function of frequency for a sound-speed scatterer located at zs = 0.9 R. For this target depth closer to the surface, the rapid variations disappear above 3 mHz, due to the larger contribution of high-degree modes which are not resolved in frequency space because of their short lifetimes.

thumbnail Fig. 7.

Signal, noise, and S/N as function of frequency for a sound-speed scatterer located at zs = 0.7 R. Here we show the result for a frequency range of 2–7 mHz. The rapid changes are not due to numerical inaccuracies. The red dots indicate the values at frequency 2.4000 mHz.

thumbnail Fig. 8.

S/N as in Fig. 7, but for a sound-speed scatterer located closer to the surface at zs = 0.9 R.

The kernel function at frequency 2.4000 mHz for zs = 0.7 R is shown in Fig. 9; this particular frequency corresponds to the peak indicated in Fig. 7 with a red dot. We see that this kernel is much less localized around the scattering point than the kernel at 3 mHz (Fig. 2).

thumbnail Fig. 9.

Meridional slice of the sound-speed kernel K αβ c $ {\cal K}_{\alpha \beta }^c $ with zs = 0.7 R computed at the low frequency of 2.4000 mHz, which corresponds to the spike with a red dot in Fig. 7. This kernel displays oscillations and is not peaked as much as the 3 mHz kernel from Fig. 2.

8. Conclusion

We derived a framework for computing the expected signal and the noise level in PB helioseismic holography. The same framework could be used to interpret LB data, including phase-sensitive data.

Porter-Bojarski holography requires knowledge of the wave field, ψ = ρ1/2c2∇⋅ξ, and its normal derivative on the solar surface, ∂nψ. With this definition of ψ, the Green’s function used in the definition of PB integrals solves a well-defined Helmholtz-like equation, which we solve numerically (Gizon et al. 2017). The need for finite-frequency Green’s functions was demonstrated in LB holography by Pérez Hernández & González Hernández (2010). In the numerical examples shown in the previous section, we assumed that we have full knowledge of ∂nψ on the surface. In practice, we do not observe directly the normal derivative of the wave field; it must be approximated. According to complementary simulations (not shown here), this can be achieved by using the approximate outgoing radiation condition ∂nψ = iknψ derived by Barucq et al. (2018).

We found that, for a sufficiently large pupil, scatterers can be imaged at a resolution that is very close to half the local wavelength, λ/2. This confirms the claim by Lindsey & Braun (1997, 2000a) that helioseismic holography is diffraction-limited. In that sense, helioseismic holography is superior to deep-focusing time-distance helioseismology, which gives lower spatial resolution (Munk 2001; Pourabdian et al. 2018). For large pupils, we found that the S/N in PB images does not vary much with depth in the convection zone, when a perturbation in sound-speed fills a volume corresponding to the holographic resolution.

Averaging over frequencies improves the S/N. For a scatterer at the bottom of the convection zone, the signal and the noise vary smoothly with frequency above 4 mHz (see Fig. 7). At lower frequencies, however, the signal varies rapidly with frequency (due to contributions from individual long-lived p modes) and it is not obvious how the signal should be averaged. A specific analysis of low-frequency data is required, especially for deep scatterers.

We found that the S/N in PB holography is maximum around 3.7 mHz for zs = 0.7 R (resp. at 4.3 mHz for zs = 0.9 R). There is no indication in our calculations that there is a benefit in using the frequencies above the acoustic cutoff (unlike predictions by Ruzmaikin & Lindsey 2003). The S/N drops to very small values above 5 mHz. One may ask if this drop is somehow compensated by the increase in spatial resolution at high frequencies. The answer is negative. Our calculations indicate that noise has a horizontal correlation length that is about half the local wavelength. Far too few independent measurements are available at high frequencies to recover an acceptable S/N by horizontal spatial averaging.

Our synthetic data do not contain a convective background. The effect of this background on S/Ns in holography should be studied. Future work should also investigate the performance of PB holography for target locations that are away from the axis of the pupil, especially for farside imaging applications.

Acknowledgments

The theoretical framework was developed by L.G. and D.F. at Mathematisches Forchungsinstitut Oberwolfach in May 2017. The numerical computations were performed by D.Y. using the Montjoie solver. D.Y. is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen. We thank M. Duruflé and J. Chabassier from Inria Magique-3D for the helioseismology-related developments of Montjoie. We also thank Chris Hanson from NYUAD for the fine tuning of the model power spectrum of solar oscillations. L.G. acknowledges support from NYUAD Institute grant G1502. The computing resources were provided in part by the German Data Center for SDO, a project funded by the German Aerospace Center (DLR).

References

  1. Barucq, H., Chabassier, J., Duruflé, M., Gizon, L., & Leguèbe, M. 2018, Math. Modell. Numer. Anal., 52, 945 [Google Scholar]
  2. Besliu-Ionescu, D., Donea, A., & Cally, P. 2017, Sun Geosphere, 12, 59 [NASA ADS] [Google Scholar]
  3. Birch, A. C., Braun, D. C., Hanasoge, S. M., & Cameron, R. 2009, Sol. Phys., 254, 17 [NASA ADS] [CrossRef] [Google Scholar]
  4. Birch, A. C., Braun, D. C., Leka, K. D., Barnes, G., & Javornik, B. 2013, ApJ, 762, 131 [NASA ADS] [CrossRef] [Google Scholar]
  5. Birch, A. C., & Gizon, L. 2007, Astron. Nachr., 328, 228 [NASA ADS] [CrossRef] [Google Scholar]
  6. Birch, A. C., Parchevsky, K. V., Braun, D. C., & Kosovichev, A. G. 2011, Sol. Phys., 272, 11 [NASA ADS] [CrossRef] [Google Scholar]
  7. Birch, A. C., Schunker, H., Braun, D. C., et al. 2016, Sci. Adv., 2, e1600557 [NASA ADS] [CrossRef] [Google Scholar]
  8. Braun, D. C., & Birch, A. C. 2008a, ApJ, 689, L161 [NASA ADS] [CrossRef] [Google Scholar]
  9. Braun, D. C., & Birch, A. C. 2008b, Sol. Phys., 251, 267 [NASA ADS] [CrossRef] [Google Scholar]
  10. Braun, D. C., Birch, A. C., & Lindsey, C. 2004, in SOHO 14 Helio- and Asteroseismology: Towards a Golden Future, ed. D. Danesy, ESA SP, 559, 337 [NASA ADS] [Google Scholar]
  11. Braun, D. C., Birch, A. C., Benson, D., Stein, R. F., & Nordlund, Å. 2007, ApJ, 669, 1395 [NASA ADS] [CrossRef] [Google Scholar]
  12. Christensen-Dalsgaard, J., Däppen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [CrossRef] [Google Scholar]
  13. Deubner, F.-L., & Gough, D. 1984, ARA&A, 22, 593 [NASA ADS] [CrossRef] [Google Scholar]
  14. Devaney, A., & Porter, R. 1985, JOSA A, 2, 2006 [NASA ADS] [CrossRef] [Google Scholar]
  15. Duruflé, M. 2006, PhD Thesis, ENSTA ParisTech, France [Google Scholar]
  16. Fournier, D., Gizon, L., Hohage, T., & Birch, A. C. 2014, A&A, 567, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Fournier, D., Leguèbe, M., Hanson, C. S., et al. 2017, A&A, 608, A109 [CrossRef] [EDP Sciences] [Google Scholar]
  18. Fournier, D., Hanson, C. S., Gizon, L., & Barucq, H. 2018, A&A, 616, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Gizon, L., & Birch, A. C. 2002, ApJ, 571, 966 [NASA ADS] [CrossRef] [Google Scholar]
  20. Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gizon, L., & Birch, A. C. 2005, Liv. Rev. Sol. Phys., 2, 6 [Google Scholar]
  22. Gizon, L., Barucq, H., Duruflé, M., et al. 2017, A&A, 600, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Lamb, H. 1909, Proc. London Math. Soc., 7, 122 [Google Scholar]
  24. Liewer, P. C., González Hernández, I., Hall, J. R., Lindsey, C., & Lin, X. 2014, Sol. Phys., 289, 3617 [Google Scholar]
  25. Lindsey, C., & Braun, D. C. 1990, Sol. Phys., 126, 101 [Google Scholar]
  26. Lindsey, C., & Braun, D. C. 1997, ApJ, 485, 895 [Google Scholar]
  27. Lindsey, C., & Braun, D. C. 2000a, Sol. Phys., 192, 261 [NASA ADS] [CrossRef] [Google Scholar]
  28. Lindsey, C., & Braun, D. C. 2000b, Science, 287, 1799 [NASA ADS] [CrossRef] [Google Scholar]
  29. Lindsey, C., Birch, A. C., & Donea, A. C. 2006, in Solar MHD Theory and Observations: A High Spatial Resolution Perspective, eds. J. Leibacher, R. F. Stein, & H. Uitenbroek, ASP Conf. Ser., 354, 174 [NASA ADS] [Google Scholar]
  30. Munk, J. M. 2001, PhD Thesis, University of Aarhus, Aarhus, Denmark [Google Scholar]
  31. Pérez Hernández, F., & González Hernández, I. 2010, ApJ, 711, 853 [NASA ADS] [CrossRef] [Google Scholar]
  32. Porter, R. P. 1969, Phys. Lett. A, 29, 193 [NASA ADS] [CrossRef] [Google Scholar]
  33. Porter, R. P., & Devaney, A. J. 1982, JOSA, 72, 327 [NASA ADS] [CrossRef] [Google Scholar]
  34. Pourabdian, M., Fournier, D., & Gizon, L. 2018, Sol. Phys., 293, 66 [NASA ADS] [CrossRef] [Google Scholar]
  35. Ruzmaikin, A., & Lindsey, C. 2003, in GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. Sawaya-Lacoste, ESA SP, 517, 71 [NASA ADS] [Google Scholar]
  36. Skartlien, R. 2001, ApJ, 554, 488 [NASA ADS] [CrossRef] [Google Scholar]
  37. Skartlien, R. 2002, ApJ, 565, 1348 [NASA ADS] [CrossRef] [Google Scholar]
  38. Tsang, L., Ishimaru, A., Porter, R. P., & Rouseff, D. 1987, JOSA A, 4, 1783 [NASA ADS] [CrossRef] [Google Scholar]
  39. Yang, D. 2018, Sol. Phys., 293, 17 [NASA ADS] [CrossRef] [Google Scholar]
  40. Zharkov, S., Green, L. M., Matthews, S. A., & Zharkova, V. V. 2013, Sol. Phys., 284, 315 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Possible wave propagators.

Table 2.

Proposed propagators and associated pupils (Hα, A) and (Hβ, A′) for different types of perturbations.

All Figures

thumbnail Fig. 1.

Pupil geometries used to compute sound-speed or source kernels (P) and flow kernels (Qs), also see Table 2.

In the text
thumbnail Fig. 2.

Meridional slice through the sound-speed kernel K αβ c $ {\cal K}_{\alpha \beta }^c $(x, r​s) computed in Model S at a single frequency of 3 mHz, in units of 10−34 kg m−3 s−3. Both the real (left panel) and imaginary (right panel) parts of the kernel are shown. The scatterer at zs = 0.7 R is indicated by the crosses. The observation pupil P is a polar cap of full angular size 120°. We notice the “ghost values” at the antipode.

In the text
thumbnail Fig. 3.

Cut along the z axis through the sound-speed kernel from Fig. 2. The scatterer is at zs = 0.7 R.

In the text
thumbnail Fig. 4.

Left panel: image intensity |𝔼[δIαβ(x)]| at a single frequency of 3 mHz, as a function of position x along the z-axis. The sound-speed perturbation (see Eq. (54)) is placed at two different positions along the z-axis, zs = 0.7 R (red) and 0.9 R (blue). The standard deviation of the noise Var I α β , 0 ( x ) $ \sqrt{ \mathrm{Var} \langle I_{\alpha\beta, 0} (\boldsymbol{x}) \rangle } $ is given by the black curve. We note that the jagged aspect of the curves is not due to numerical inaccuracies. Right panel: Image intensity and noise level after averaging over 101 frequencies uniformly distributed in the interval from 2.75 to 3.25 mHz. The frequency resolution is 5 μHz, implying an observation duration of T = 55.5 h. A horizontal line segment is plotted at each depth to mark half of the local wavelength.

In the text
thumbnail Fig. 5.

Radial and horizontal widths of the frequency-averaged sound-speed kernel |⟨𝒦c⟩| as functions of scattering position. These are close to half the local wavelength at 3 mHz (dotted line).

In the text
thumbnail Fig. 6.

S/N in PB image intensity for a 10% sound-speed perturbation over a volume Vs(zs) placed at zs along the polar axis (Eq. (54)). The results are shown at a single frequency of 3 mHz (solid curve) and after averaging over 101 frequencies in the interval from 2.75 to 3.25 mHz (dashed curve).

In the text
thumbnail Fig. 7.

Signal, noise, and S/N as function of frequency for a sound-speed scatterer located at zs = 0.7 R. Here we show the result for a frequency range of 2–7 mHz. The rapid changes are not due to numerical inaccuracies. The red dots indicate the values at frequency 2.4000 mHz.

In the text
thumbnail Fig. 8.

S/N as in Fig. 7, but for a sound-speed scatterer located closer to the surface at zs = 0.9 R.

In the text
thumbnail Fig. 9.

Meridional slice of the sound-speed kernel K αβ c $ {\cal K}_{\alpha \beta }^c $ with zs = 0.7 R computed at the low frequency of 2.4000 mHz, which corresponds to the spike with a red dot in Fig. 7. This kernel displays oscillations and is not peaked as much as the 3 mHz kernel from Fig. 2.

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.