© 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

(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

(2)

to recast the wave equation into a Helmholtz-like equation

(3)

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

(4)

where the squared acoustic cut-off frequency is

(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 under the conditions that the flow conserves mass and that it does not cross the boundary (u_n = 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

(6)

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

(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 k_n from Eq. (7) is given by

(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 k_n 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):

(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 (G₀) and incoming () Green’s functions defined in a reference medium with ρ₀, c₀, γ₀, and u₀ = 0:

(10)

(11)

with

(12)

where k₀ is k in the reference medium and k_n is from Eq. (8). The Green’s functions or are backward propagators (cf. egression), while ( is a forward propagator (cf. ingression). When the surface A is closed, it is equivalent to use and Im G₀ (Devaney & Porter 1985). Tsang et al. (1987) proposed 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

(13)

When , we have

(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

(15)

For the case α = β we define

(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).

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

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 ρ₀, c₀, γ₀, and u₀ = 0. The perturbations to the sources of excitation are described through the perturbations to the source covariance,

(17)

where

(18)

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

(19)

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

(20)

(21)

The perturbed wave operator is

(22)

with

(23)

In terms of the Green’s function G₀, we have

(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

(25)

such that the PB integral is given by

(26)

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

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

(27)

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

(28)

where

(29)

4.3. Perturbations to the image intensity

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

(30)

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

(31)

(32)

Using the definition of the source kernels, we obtain

(33)

where

(34)

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

(35)

The perturbation to the image intensity is given by

(36)

where

(37)

(38)

(39)

The kernels for δk² and δk^2* can be combined using Eq. (23) to obtain kernels for sound-speed , density and attenuation . For example,

(40)

with

(41)

4.4. Choice of the source covariance

In order to compute the above kernels, we still need to choose the source covariance function M₀ in order to define the reference cross-covariance C₀ 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

(42)

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

(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:

(44)

where

(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:

(46)

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

(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 χ² 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 [ω₀ − Δω/2, ω₀ + Δω/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

(48)

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

(49)

(50)

(51)

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

The expected S/N is thus

(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 T^1/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 γ = γ₀|ω/ω₀|^5.77, where γ₀/2π = 4.29 μHz and ω₀/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 k_n (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 at each spherical harmonic degree l, to obtain:

(53)

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

7.2. Sound-speed kernels

The sound-speed kernel is computed using Eq. (41) and the definition of . 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 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.

Fig. 2. Meridional slice through the sound-speed kernel (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.

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 along the polar axis, we consider a localized increase in sound speed of 10% over a volume V_s, such that the signal (Eq. (40)) may be written as

(54)

The volume V_s is that of a ball of diameter λ(r_s)/2 = π/[Re k(r_s, ω₀)] with ω₀/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 V_s. 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 z_s = 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.

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 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 <  z_s/R_⊙ <  0.98. Thus helioseismic holography is a diffraction-limited imaging technique as suggested by Lindsey & Braun (1997).

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 C₀ 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 <  z_s/R_⊙ <  0.98 for a pupil of angular size 120°. We note that the ghost at −z_s is much below the noise level.

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 z_s = 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 z_s = 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.

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.

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 z_s = 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).

Fig. 9. Meridional slice of the sound-speed kernel 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/2c²∇⋅ξ, 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ψ = ik_nψ 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 z_s = 0.7 R_⊙ (resp. at 4.3 mHz for z_s = 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

All Tables

All Figures

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

Fig. 2. Meridional slice through the sound-speed kernel (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

	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

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

	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

	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

	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

	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

	Fig. 9. Meridional slice of the sound-speed kernel 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