Quantitative passive imaging by iterative holography: the example of helioseismic holography

In passive imaging, one attempts to reconstruct some coefficients in a wave equation from correlations of observed randomly excited solutions to this wave equation. Many methods proposed for this class of inverse problem so far are only qualitative, e.g. trying to identify the support of a perturbation. Major challenges are the increase in dimensionality when computing correlations from primary data in a preprocessing step, and often very poor pointwise signal-to-noise ratios. In this paper, we propose an approach that addresses both of these challenges: it works only on the primary data while implicitly using the full information contained in the correlation data, and it provides quantitative estimates and convergence by iteration. Our work is motivated by helioseismic holography, a well-established imaging method to map heterogenities and flows in the solar interior. We show that the back-propagation used in classical helioseismic holography can be interpreted as the adjoint of the Fréchet derivative of the operator which maps the properties of the solar interior to the correlation data on the solar surface. The theoretical and numerical framework for passive imaging problems developed in this paper extends helioseismic holography to nonlinear problems and allows for quantitative reconstructions. We present a proof of concept in uniform media.


Introduction
In this paper, we consider passive imaging problems described by a linear time-harmonic wave equation L[q]ψ = s with a random source s and some unknown coefficient q, which is the quantity of interest.We assume that E [s] = 0 such that E ψ = 0 by linearity of L [q].Solutions ψ to this wave equation are observed on part of the boundary Γ = ∂Ω of a domain Ω for many independent realizations of s.Thus we can approximately compute the cross-covariance Our aim is to determine the unknown parameter q given noisy observations of C or the corresponding integral operator (C f )(x 1 ) := Γ C(x 1 , x 2 ) f (x 2 ) dx 2 .If Tr Γ is the trace operator onto Γ, then straightforward calculations show that the forward operator mapping q to C = Cov[Tr Γ ψ] is given by (Recall that the covariance operator Cov[v] ∈ L(X) of a random variable v with values in a Hilbert-space X is defined implicitely by Cov(⟨v, ψ⟩ X , ⟨v, φ⟩ X ) = ⟨Cov[v]ψ, φ⟩ X for all φ, ψ ∈ X.)An early and influential reference on passive imaging is the work of Duvall et al. [1] on time-distance helioseismology.Later, passive imaging has also been used in many other fields such as seismology ( [2]), ocean acoustics ( [3]), and ultrasonics ( [4]).We refer to the monograph [5] by Garnier & Papanicolaou for many further references.Concerning the uniqueness of passive imaging problems, we refer to [6] for results in the time domain and to [7,8,9,10] for results in the frequency domain.For the unique recovery of the source and the potential from passive far-field data, we refer to [11].
Local helioseismology analyzes acoustic oscillations at the solar surface in order to reconstruct physical quantities (subsurface flows, sound speed, density) in the solar interior (e.g.[12] and references therein).Since solar oscillations are excited by near-surface turbulent convection, it is reasonable to assume random, non-deterministic noise terms.In this paper, we will describe sound propagation in the solar interior by a scalar time-harmonic wave equation and study the passive imaging problem of parameter reconstruction from correlation measurements.
Very large data sets of high-resolution solar Doppler images have been recorded from the ground and from space over the last 25 years.This leads to a five-dimensional (2 2 spatial dimensions and 1 temporal dimension) cross-correlation data set on the solar surface, which cannot be stored and analyzed all at once.In traditional approaches, like timedistance helioseismology, the cross-correlations are reduced to a smaller number of observable quantities, such as travel times ( [1]) or cross-correlation amplitudes (e.g.[13], [14], [15]).Since the reduction to these quantities leads to a loss in information, we are interested in using the whole cross-correlation data throughout the inversion procedure, stepping forward to full waveform inversions.
Helioseismic holography, a technique within the field of local helioseismology, has proven to be a powerful tool for studying various aspects of the Sun's interior.It operates by propagating the solar wavefield backward from the surface to specific target locations within the Sun ( [16]).A notable success of helioseismic holography is the detection of active regions on the Sun's far side (e.g.[17], [18], [19]).Furthermore, helioseismic holography is used in many other applications, e.g. to study the subsurface structure of sunspots ( [20], [21], [22]), wave absorption in magnetic regions ( [23], [24], [25]), and seismic emission from solar granules [26].The main idea of helioseismic holography is the back-propagation ("egression") of the wavefield at the solar surface ( [27]).Improvements have been proposed in the choice of backward propagators (e.g. using Porter-Bojarski holograms ( [28], [29]).Helioseismic holography has a strong connection to conventional beam forming, where imaging functionals similar to the holographic back-propagation occur (e.g.[30]).In contrast to these approaches, we will achieve improvements by iterations.
In the present paper, we connect holographic imaging methods to iterative regularization methods.This way, holography can be extended to a full converging regularization method.This approach was successfully applied to inverse source problems in aeroacoustics ( [10]) and is extended in this work to parameter identification problems.
The organization of the paper is as follows.In Section 2 we introduce a generic model for the forward problem.In Section 3 we establish foundations of our functional analytic setting by establishing sufficient conditions under which the diagonal of an integral operator is well defined, using Schatten class properties of embedding operators.With this we compute the Fréchet derivative of the forward operator and its adjoint in section 4. Next, we discuss the algorithm of iterative holography in Section 5. Based on the analysis of Sections 2-Appendix B we then introduce forward operators in local helioseismology, their derivatives and adjoints in Section 6.Then we discuss iterative helioseismic holography as an extension of conventional helioseismic holography in Section 7, and demonstrate its performance in numerical examples with simulated data in Section 8 before we end the paper with conclusions in Section 9. Some technical issues are discussed in three short appendices.

A model problem
We first present the main ideas of this paper for a generic scalar time-harmonic wave equation.
Let Ω 0 ⊂ Ω be a smooth, bounded domain in R d and let Γ ⊂ Ω \ Ω 0 the hypersurface on which measurements are performed.Γ may be part of the boundary ∂Ω or it may be contained in the interior of Ω.Moreover, consider the parameters Assume that the excitation of wavefields ψ in R d by random sources s, which are supported in Ω 0 , is described by the model for the outward pointing normal vector n on ∂Ω and some operator B ∈ L H 1/2 (∂Ω) → H −1/2 (∂Ω) .(Here and in the following L (X, Y) denotes the space of bounded linear operators between Banach spaces X and Y.) Typically, B is some transparent boundary condition, e.g.Bψ = ikψ for ∂Ω = S d−1 .We may also choose Bψ := DtN ψ with an exterior Dirichlet-to-Neumann map for the Helmholtz equation with the Sommerfeld radiation condition.In this case, equation ( 2) is equivalent to a problem posed on R d with the Sommerfeld radiation condition.
Assumption 1. Suppose that for some B 0 ∈ L H 1/2 (∂Ω), The conditions (3c)-(3e) are obviously satisfied for Bζ := ikζ, and they also hold true if B is the exterior Dirichlet-to-Neumann map on a sphere or a circle (see [31,32]).Throughout this paper we denote by H s 0 (Ω) the closure of the space of distributions on Ω in H s (R d ).For a Lipschitz domain, we have the duality H s (Ω) * = H −s 0 (Ω) ( [33,Thm 3.30]).
Proposition 1.Under Assumption 1 the problem (2) is well posed in the sense that for all s ∈ H −1 0 (Ω) there exists a unique ψ ∈ H 1 (Ω) satisfying (2) in the weak sense, and ψ depends continuously on s with respect to these norms.
Proof.We only sketch the proof, which is a straightforward modification of similar proofs in [31,32].The weak formulation of Problem ( 2) is given by To show that for s = 0 this variational problem only has the trivial solution, we choose ϕ = ψ and take the imaginary part.Noting that Im(−2iA 2  and using a partial integration and (3b), we obtain It follows from (3a) and (3c) that both sides must vanish.Hence, Tr ∂Ω ψ = 0.By elliptic regularity, ψ ∈ H 2 (Ω) is also a strong solution to (2) with ∂ψ ∂n = 0 on ∂Ω.Due to vanishing Cauchy data on ∂Ω, ψ may be extended by 0 as a strong solution of the wave equation to the exterior of Ω.Now it follows from unique continuation results (see [34,Theorem 4.2]) that ψ vanishes identically.
Using Assumptions (3d) and (3e), it can be shown that the sesquilinear form of the variational formulation is coercive up to a compact perturbation.Therefore, the operator representing this sesquilinear form is Fredholm of index 0.By uniqueness, it is boundedly invertible.

□
If we write the solution operator is the Green's function, which may also be characterized by For certain random processes of interest, s does not belong to H −1 0 (Ω) almost surely.E.g., white noise is in H −s 0 (Ω) almost surely if and only if s > d/2.Nevertheless, the solution formula may still make sense if G v,A (x, •) is sufficiently smooth on Ω 0 .This is always the case if the support of s and v, A are disjoint or if s ∈ H −1 0 (Ω) almost surely, which is typically true if s is spatially correlated.Otherwise, we have to impose smoothness conditions on v and A such that G v,A is sufficiently smooth and G has suitable mapping properties.
Assume we have observations Tr Γ ψ 1 , . . ., Tr Γ ψ N where ψ j solves (2) for independent samples s 1 , . . ., s N of s.As E [s] = 0, we have E Tr Γ ψ j = 0, and we can compute the correlations by This is an unbiased estimator of the covariance C will be the forward operator of our inverse problem.Recall that if belongs to the space of Hilbert-Schmidt operators HS L 2 (Γ) on L 2 (Γ), and Therefore, HS L 2 (Γ) is the natural image space of the forward operator.It is a Hilbert space with an inner product ⟨T, S ⟩ HS = tr(S * T ).Here tr(K) denotes the trace of a linear operator K : H → H in a separable Hilbert space H defined by tr(K) := ∞ j=1 ⟨Ke j , e j ⟩ H for any is an orthonormal basis {e j : j ∈ N} of H.
Let us also consider the case that in addition to sources s in the interior of Ω there are sources s ∂Ω on the boundary ∂Ω.Such sources generate a field ψ(x) = ∂Ω G v,A (x, y)s ∂Ω (y) dy.Its restriction to Ω 0 is given by which is the single layer potential operator for Γ = ∂Ω.It is easy to see that K admits a factorization K = Tr Γ G v,A Tr * ∂Ω in the spaces [35, Thm.1(iii)]).Therefore, in the presence of boundary sources the measured covariance operator is given by Often one assumes that the source process s is spatially uncorrelated and for some source strength S ∈ L ∞ (Ω 0 ), where M S denotes the multiplication operator M S f := S • f .If S is treated as an additional unknown, the forward operator becomes Of course, we could also assume that s ∂Ω is spatially uncorrelated and treat its source strength as a further unknown, but for the sake of notational simplicity, we assume that Cov[s ∂Ω ] ∈ L L 2 (∂Ω) is known.We first study the continuity and Fréchet differentiability of G v,A with respect to the parameters (v, A).We will assume that v and A are known in Ω \ Ω 0 .Let (v ref , A ref ) ∈ B k be some reference solution.Then the set B k of admissible parameters in Assumption 1 satisfies where A is well-defined and continuous, and Fréchet differentiable in the interior of Proof.Again, we only sketch the proof and refer to [31, §5.3] for a more detailed proof of a similar result.Let L v,A : H 1 (Ω) → H −1 0 (Ω) denote the operator associated to the sesquilinear form in the weak formulation (4) A is continuous and affine linear in the parameters.As L ′ v,A (∂v, ∂A) = −2iM ∂A • ∇ + M ∂v , the result follows from the continuity of operator inversion and the formula for its derivative,

Diagonals of operator kernels
The present section serves as a preparation for computing adjoints of the Fréchet derivative of the forward operator defined by (10).A crucial step will be the characterization of adjoints of the mapping S → M S (in a sense to be specified later).
In the discrete setting, M S corresponds to diagonal matrices diag(S ) ∈ C d×d with diagonal S .The adjoint of the mapping with respect to the Frobenius norm is given by where Diag(A) ∈ C d denotes the diagonal of the matrix A ∈ C d×d .We wish to generalize this to an infinite dimensional setting, with the Frobenius norm replaced by the Hilbert-Schmidt norm.Recall that any operator A ∈ HS L 2 (Ω) has a Schwartz kernel A ∈ L 2 (Ω × Ω) such that (Aφ)(x) = Ω A(x, y)φ(y) dy and ∥A∥ HS = ∥A∥ L 2 .It is tempting to define (Diag A)(x) := A(x, x).However, as A is only a L 2 -function and the diagonal {(x, x) : x ∈ Ω} ⊂ Ω × Ω has measure zero, the restriction of A to the diagonal is not well-defined.
To address this problem, we first recall that for Hilbert spaces X, Y and p ∈ [1, ∞) the p-Schatten class S p (X, Y) consists of all compact operator A ∈ L (X, Y) for which the singular values σ j (A) (counted with multiplicity) form a ℓ p sequence.S p (X, Y) is a Banach space equipped with the norm ∥A∥ S p := ( j σ j (A) p ) 1/p .S 2 (X, Y) coincides with HS (X, Y).We write S p (X) := S p (X, X).The elements of S 1 (X) are called trace class operators.For such operators, the trace tr(A) := k ⟨Ae j , e j ⟩ is well-defined for any orthonormal basis {e k } of X, and | tr(A)| ≤ ∥A∥ S 1 .
Let us first recall Mercer's theorem: It states that for a positive definite operator A with continuous kernel A, we have tr A = Ω A(x, x) dx and A(x, x) ≥ 0 for all x.Since not all (positive semidefinite) Hilbert-Schmidt operators are trace class, we cannot expect that x → A(x, x) belongs to L 1 (Ω) for general Hilbert-Schmidt operators.However, with the help of Mercer's theorem, we can show the following result.Diag : for all operators A ∈ S 1 (L 2 (Ω)) with continuous kernel A. Moreover, Eq. ( 12) is shown in [36,Thm. 3.5] where it is also shown that Diag(A) can be constructed by local averaging, but the first part is not explicitly stated.We sketch an alternative, more elementary proof: Proof of Proposition 3. If A is positive semidefinite, it may be factorized as A = B * B with B ∈ HS L 2 (Ω) and ∥A∥ S 1 = ∥B∥ 2 S 2 , e.g., by choosing B = A 1/2 .By density of C(Ω × Ω) in L 2 (Ω × Ω), there exists a sequence (B n ) converging to B in L 2 (Ω × Ω).For the corresponding operators B n it follows that lim n→∞ ∥B n −B∥ HS = 0 and lim n→∞ ∥A n −A∥ S 1 = 0 for A n := B * n B n (see Prop. 4, part ii below).Thus, we have constructed a sequence of positive semidefinite operators with continuous kernels converging to A in S 1 L 2 (Ω) , and the statement follows from the classical Mercer theorem.
We decompose a general A ∈ S 1 L 2 (Ω) as linear combination of trace class operators: We start with A = Re(A) + i Im(A) where Re(A) := 1 2 (A + A * ) and Im(A) := 1 2i (A − A * ).

can be written as a linear combination of positive semi-definite trace class operators:
we can apply the first proven special case to all P j to obtain the result.□ To speak of an adjoint of the operator M : S → M S , we have to treat M S in some space with a dual pairing.We will use Hilbert-Schmidt spaces between suitable Sobolev spaces.
of Hilbert-Schmidt spaces with dual pairing, given by ⟨A, B⟩ HS := tr(B * A) for A ∈ HS (V, V ′ ) and B ∈ HS (V ′ , V).
We give some preliminary results on p-Schatten class embeddings.
Proposition 4. Let Ω ⊂ R d be a bounded Lipschitz domain and let X, Y, Z be Hilbert spaces.
Then the following holds true: (i) The Sobolev embedding: j : Then, BA ∈ S r (X, Z) and we have the bound: and we have the bounds: Proof.Part (i) follows from Theorem 1 of [37].Part (ii) and (iii) follow from Lemma 16.7 of [38].
Let A ∈ S p (X, Y) and B ∈ S q (X ′ , Y ′ ).By part (ii) and the boundedness of the trace in S 1 ([38,Lemma 16.23]), we get: By the Hahn-Banach theorem, we have the sequence S p (X, Y) is a uniformly convex Banach space ([39, Chapter 5]) and therefore reflexive by Milman-Pettis theorem ( [40]).Hence, S p (X, Y) ′′ = S p (X, Y) and the assertion follows.□ Using this proposition, we can prove that multiplication operators are Hilbert-Schmidt in suitable Sobolev spaces: Lemma 5. Let Ω ⊂ R d be a bounded Lipschitz domain, S ∈ L ∞ (Ω), and s > d/4, s−1/2 N 0 .(In particular, for d ∈ {2, 3} we may choose s = 1.)Then M S ∈ HS H s (Ω), H −s 0 (Ω) , and the following mapping is continuous Then, we consider M S in the function spaces: By Proposition 4, part i, the embedding j is an element of the Schatten class The continuity of M follows from the continuity of the mapping: S → MS .□ Lemma 6.Under the assumptions of Lemma 5, the adjoint operator M * : HS Proof.Let f ∈ L ∞ (Ω), j : H s (Ω) → L 2 (Ω), and A ∈ HS H −s 0 (Ω), H s (Ω) .It follows from Proposition 4 that Ã := jA j * ∈ S 1 L 2 (Ω) .We identify Ã and A, i.e. a more precise formulation of ( 14) is M * (A) = Diag( Ã) for all A. By the density result established at the beginning of the proof of Proposition 3, it suffices to establish the relation for operators A with continuous kernel A. Choosing an orthonormal basis {e k : k ∈ N} of L 2 (Ω), we obtain Since Ω is bounded and hence A(y, In Lemma 2 we consider M ∂v and M ∂A in the following function spaces: The multiplication operator M ∂v has been discussed in Lemma 5 and Lemma 6.Although for M ∂A we have less regularity, the following analogs still hold true: and the following map is continuous: Proof.In this proof, j will denote the embedding H 1 (Ω 0 ) → L 2 (Ω 0 ) and recall from Proposition 4, part i that j ∈ S 4 H 1 (Ω 0 ), L 2 (Ω 0 ) and hence . As in Lemma 6, the assertion now follows.□

Fréchet derivative and adjoint of the forward operator
A characterization of the adjoint of S → Tr Γ GM S G * Tr * Γ was given in [10].There, a characterization of the adjoint of S → M S in a functional analytic framework was circumvented, resulting in a rather technical formulation of the result.
With the results of the previous section, the proof of the following central results is now mostly straightforward.Theorem 8. Assumptions 1 and 2 hold true for some wave number k ∈ C and d ∈ {2, 3}.Let Then the following holds true: (i) The forward operator (10) is well-defined and continuous as a mapping and it is Fréchet differentiable on the interior of B. The derivative We consider the factors defining C 1 [v, A, S ] in the following function spaces: Here M S : Lemma 5, and all other operators are bounded.By part (iii) of Proposition 4, it follows that Similarly, we consider the factors defining C 2 [v, A, S ] in the following function spaces: By part (i) of Proposition 4, every embedding is an element of S 8 .By part ii, iii of Proposition 4, it follows that Together with Lemma 2, it follows that C is continuous.Fréchet differentiability and the formula for the derivative follow from Lemma 2 and the chain rule.
Part (ii): If X 1 , . . ., X 4 are Hilbert spaces, A ∈ L (X 1 , X 2 ) and B ∈ L (X 3 , X 4 ), then a straightforward computation shows that the adjoint of the linear mapping HS (X 2 , X 3 ) → HS (X 1 , X 4 ), T → BT A is given by the mapping HS (X 1 , X 4 ) → HS (X 2 , X 3 ), S → B * SA * and that Re ∈ L HS X j is a self-adjoint projection operator.Note that the Fréchet derivative and its adjoint in Theorem 8 can be reformulated as These propagators H α , H β have a physical interpretation in helioseismology that will be discussed in Section 7.1.

On the algorithmic realization of iterative regularization methods
For notational simplicity, we will use q = (v, A, S ) throughout this section.Then we formally have to solve the operator equation Tr Γ ψ n ⊗ Tr Γ ψ n .

Avoiding the computation of Corr
Computing Corr from the primary data Tr Γ ψ n in a preprocessing step drastically increases the dimensionality of the data.In helioseismology, the data set with the best resolution consists of Doppler images of size 4096×4096.This leads to approximately 10 14 independent two-point correlations, at each frequency.Hence, the surface cross-correlation is a noisy five-dimensional data set of immense size, which is infeasible to use in inversions directly.Moreover, these two-point correlations are extremely noisy.In traditional approaches such as time-distance helioseismology, one usually reduces the cross-correlation in an a priori step to a smaller number of physically interpretable quantities with an acceptable signal-to-noise ratio.However, this leads to a significant loss of information, see [42].
To use the full information content of Corr without the need to ever compute these correlations, we exploit the fact that the adjoint of the Fréchet derivative of the forward operator is of the form C ′ [q] * D = Diag(H q α * Re(D)H q β ), see eq. ( 16).As Re(Corr) = Corr, pulling the sum outside yields We will show in Section 7.1 that traditional helioseismic holography can be interpreted as the application of C ′ [q] * to Corr.Since 1 N N n=1 Diag(. . . ) can be interpreted as computing diagonal correlations of the back-propagated signals, eq. ( 17) may be paraphrased as first back-propagating signals and then correlating them, rather than vice versa.

Iterative regularization methods without image space vectors
For ill-posed inverse problems, the adjoint of the linearized forward operator is typically a bad approximation of the inverse.To obtain a quantitative imaging method, we can improve the approximation in (17) by implementing an iterative regularization method.We will focus on the Iteratively regularized Gauss-Newton Method (IRGNM) with inner Conjugate Gradient iterations, but the discussion below also applies to other commonly used methods such as Landweber iteration or the Newton-CG method.IRGNM is defined by Here q 0 defines the initial guess.Since the image space Y of the forward operator is high dimensional, direct evaluations of C[q] and C ′ [q] must be avoided.However, IRGNM with inner CG iterations as well as other iterative regularization methods only require the operations q → C ′ [q] * C[q] and We will refer to the integral kernels K of C ′ [q] * C ′ [q] as sensitivity kernels for the normal equation.In Section 7.2 they will be described for various settings of interest in terms of forward-backward operators In our numerical tests in helioseismology reported in Section 8, the bottleneck concerning computation time is the evaluation of the Green function involved in the definitions of the propagators H α and H β .To accelerate these computations, we use separable reference Green's functions G 0 := G q 0 discussed in Appendix A and Appendix B and the corresponding integral operator G 0 as well as the algebraic identity This identity, with L q := G −1 q and L 0 := G −1 0 , is equivalent to G 0 − G q = G 0 (L q − L 0 )G q .The operators L q , L 0 : H 1 (Ω) → H −1 0 (Ω) represent the corresponding sesquilinear forms in eq. ( 4) in the proof of Proposition 1 and involve both the differential operator and the boundary condition.As both the boundary operator B and the leading order differential operator are independent of the parameters q, they cancel out, and This approach is efficient since the operator G 0 can be solved with one-dimensional code and the operator the calculation of Id +G 0 (L q − L 0 ) can typically be restricted to a supported area of L q −L 0 .Usually, we compute a pivoted LU-decomposition of Id +G 0 (L q −L 0 ) and solve for a list of input sources G(•, x), x ∈ Γ.Furthermore, we can use low-rank approximations for G 0 based on the expansions in Appendix A for solar-like medium and Appendix B for uniform medium.

Noise and likelihood modelling
In this section, we study the noise model in order to step forward to the full likelihood modeling and to create realistic noise.The main noise term is realization noise.Recall that the wavefield ψ is a realization of a Gaussian random process with covariance operator C[q].
The covariance matrix of Gaussian correlation data can be computed by Isserlis theorem [43] and is given by fourth-order correlations (e.g.[44], [45], [29]): The third term vanishes as E ψ(r 1 )ψ(r 2 ) = Ω Ω G(r r r 1 , z z z 1 )G(r r r 2 , z z z 2 )E [s(z 1 )s(z 2 )] dz 1 dz 2 and E [s(z 1 )s(z 2 )] = 0. Hence, we observe Cov(C(r 1 , r 2 ), C(r 3 , r 4 )) = E ψ(r 1 )ψ(r 2 )ψ(r 3 )ψ(r 4 ) − E ψ(r 1 )ψ(r 2 ) E ψ(r 3 )ψ(r 4 ) = E ψ(r 1 )ψ(r 3 ) E ψ(r 2 )ψ(r 4 ) = C(r 1 , r 3 )C(r 4 , r 2 ).Therefore, we can define the data covariance operator by Hence, if we choose a quadratic log-likelihood approximation, we are formally lead to replace However, with this replacement, the iteration (18) would in general not be well defined and numerically unstable since C 4 [q n ] is not boundedly invertible.Note that the operator and this cannot be bounded in infinite dimensions since C[q n ] is Hilbert-Schmidt.Hence, the mentioned replacement is not well defined and numerically unstable.Therefore, we choose a bounded, self-adjoint, positive semi-definite approximation e.g., by a truncated eigenvalue decomposition or by Lavrentiev regularization The parameter β may model the presence of measurement and modelling errors in addition to the realization noise modelled by C 4 [q n ] −1/2 .Then the iteration ( 18) is replaced by Note that for numerical efficiency it is very fortunate that the covariance operator C 4 of the correlation data has the separable structure (22).Further note from the second line of the last equation that we only need Γ n , not Γ 1/2 n .

Forward problems in local helioseismology
In this section, we discuss applications of the model problem considered in the previous sections to helioseismology.

Acoustic oscillations in the Sun
Ω 0 will denote the interior of the Sun (typically Ω 0 = B(0, R ⊙ ) with R ⊙ = 696 Mm), whereas Ω may also include parts of the solar atmosphere.The measurement region we consider is an open subset Γ of the visible surface ∂Ω 0 , accounting for the fact that in typical helioseismic applications, measurements are only available on the near side of the solar surface.Given that solar oscillations near the solar surface are primarily oriented in the radial direction [46], there is also a lack of Doppler information near the poles.This phenomenon results in leakage, causing challenges such as incomplete decoupling of normal modes of oscillation (e.g.[47,48]).In the subsequent analysis, we will exclusively work in the frequency domain.
The propagation of acoustic waves in a heterogeneous medium like the Sun can be described by the differential equation where we have ignored gravitational effects and have assumed an adiabatic approximation ( [49]).The random source term F describes the stochastic excitation of waves by turbulent motions and ζ is the Lagrangian wave displacement vector.As usual, we denote with ρ the density, c the sound speed, γ the damping, and u the flow field.If we furthermore neglect second order terms in γ, u, equation ( 24) can be converted into a Helmholtz-like equation ( [29], inspired by [50]) where The potential V is defined by The frequency ω c is recognized as the acoustic cutoff frequency.This cutoff frequency arises due to the abrupt decline in density near the solar surface and results in the trapping of acoustic modes with frequencies below the acoustic cutoff frequency.Modes with frequencies surpassing the acoustic cutoff frequency can propagate through the solar atmosphere.The conditions on c, ρ, γ, u are summarized in Assumption 3. Assumption 3. Suppose that for some B 0 ∈ L H 1/2 (∂Ω), B : H 1/2 (∂Ω) → H −1/2 (∂Ω) satisfies the conditions (3c)-(3e) (27e) For the flow field, we incorporate a mass conservation constraint (equation (27d)).Additionally, we assume that the flow field does not intersect the computational boundary (equation (27c)).Various boundary conditions, in particular radiation boundary conditions and learned infinite elements, and their efficacy are extensively discussed in [51,52,53].It is notable that the most popular choices of boundary conditions in helioseismology, such as radiation boundary conditions, Sommerfeld boundary conditions, or free boundary conditions, are incorporated in Assumption (27e).
We define the operator P that transforms the parameters in the wave equation ( 25) into the form of equation ( 2) by where In Figure 1, we present the acoustic sound speed, the density, and the scalar potential v as obtained from the Solar Model S and smoothly extended to the atmosphere.Modeling the forward problem for the Sun remains challenging due to the substantial density gradients near the surface, leading to strong variations of the scalar potential v near the solar surface.Lemma 9.The operator P, defined in equation (28), is well-defined in the sense that for all parameters (c, ρ, γ, u, S ) ∈ A c min ,ρ min ,k we have P(c, ρ, γ, u, S ) ∈ B k , and this map is continuous.
Proof.By equations (29a), (29b), we have v ∈ L ∞ (Ω, C), A ∈ W ∞ (div, Ω 0 ), and the mapping is continuous.The conditions (3c)-(3e) are obviously satisfied, and (3b) is satisfied by Assumption 27d.For condition (3a), we note that , where we have used (27d).Therefore, Because of Assumption (27b), c, ρ, ∇ρ, γ, u, ∇u are fixed at ∂Ω 0 and in the exterior.Therefore, the space of parameter perturbations is where For c min , ρ min > 0 and k ∈ C, the operator P is Fréchet differentiable in the interior of A c min ,ρ min ,k with Fréchet derivative P ′ : X P → X G given by For arguments (ṽ, Ã, S ) ∈ L 1 (Ω 0 ; C) × L 1 (Ω 0 ; R d ) × L 1 (Ω 0 ; R) ⊂ X ′ , the values of the adjoint Proof.We rephrase the potential v in the form: where The operator P is Fréchet differentiable with Fréchet derivative (30) since the terms ∂ q v, ∂ q A are well-defined for q ∈ {c, ρ, γ, u}.The claim follows with the mapping properties of (∂ q v) * , (∂ q A) * .□ In analogy to equation ( 16), we can write the Fréchet derivative in the form: The operators L q play the role of local correlation operators.The propagators and local correlation operators in the flow-free case can be read in Table 1.
Table 1.Distributional kernel of back-propagator and local correlation operator for the different parameters.The functions g 0 ρ , g 1 ρ , g 2 ρ are defined in equation (31).The coordinates are chosen such that x ∈ Γ and y ∈ Ω.Here, we assume that Cov[s ∂Ω ] = M B with B ∈ L ∞ (∂Ω) and use the notation S Ω := S + Bδ ∂Ω .

Quantity q
Propagator H αq Propagator Note that despite the fact that the adjoint with respect to the standard L 2 dual pairings takes values in negative Sobolev spaces, it is usually not necessary to deal with such functions (or distributions) numerically in iterative regularization methods.For instance, in Landweber iteration in Banach spaces, the application of the adjoint is followed by the application of a duality mapping which takes values in positive Banach spaces.For Hilbert space methods, one would choose a L 2 -based Sobolev space W s,2 with sufficiently large s and compute the adjoint with respect to the W s,2 inner product, which amounts to an evaluation of the adjoint of the embedding W s,2 → L 2 .

Source model in helioseismology
It remains to discuss the seismic source model in helioseismology.It has been shown in several settings that the cross-correlation is roughly linked to the imaginary component of the outgoing Green's function ( [30]).In helioseismology, this relation takes the form ( [49]) where Π(ω) is the source power spectrum.This relation leads to a power spectrum in good agreement with the observations [49].As outlined in [49], Equation ( 34) holds true for an outgoing radiation condition and random sources that are appropriately excited across the volume in proportion to the damping rate Moreover, there are surface integrals that persist for frequencies above the acoustic cutoff frequency, and these are dependent on the chosen boundary condition.The relationship between source power and damping rate emerges from the idea of equipartition among distinct acoustic modes ( [55]).This choice of covariance couples the source strength with wave attenuation and sound speed.Nevertheless, we consider the source strength as an additional individual parameter.This source model is included in the discussion of the previous sections.In helioseismology, the relation ( 34) is the standard choice to reduce the computational costs of the operator evaluation.Furthermore, it allows us to evaluate the back-propagator in Table 1 efficiently.

Iterative helioseismic holography
In this section, we discuss the application of the approach outlined in Section 5 to local helioseismology.We first show that it can be interpreted as an extension of conventional helioseismic holography.For this reason, we will refer to this approach as iterative helioseismic holography.We also discuss relations to other methods in local helioseismology.
In a second subsection, we will describe sensitivity kernels for the normal equation as introduced in (19) for the following three scenarios: (i) Inversion for the source strength, (ii) Inversion for scalar parameter q ∈ {ρ, c, γ}, (iii) Inversion for mass-conserved flow field u.

Relations to conventional helioseismic holography and other methods
Conventional helioseismic holography is based on the Huygens principle in the sense that the observed wavefield is described as a superposition of seismic point sources on the wavefront.This principle allows holography to propagate the correlations of acoustic waves at the solar surface forward in time ("ingression" using H β ) or backward in time ("egression" using H α ) to a pre-defined target location in the interior in order to image anomalies in the background medium (e.g.[56]).There exists a close connection to seismic migration in terrestrial seismology, which re-locates seismic events on the earth's surface in time and space, based on the wave equation (e.g.[57], [58]).Furthermore, similar back-propagators are used in conventional beamforming in aeroacoustics ( [30], [10]).
The Lindsey-Braun holographic image (see [27]) is constructed by the wave propagators where Γ 1 , Γ 2 ⊂ Γ are called pupils.In Lindsey-Braun holography the information is extracted from the so-called egression-ingression correlation for parameters q ∈ {c, ρ, u, γ} and the egression power for seismic sources where ⊗ the standard tensor product, and C v,A,S is from equation (7).The comparison of equations ( 16) and (36) shows that the adjoint of the Fréchet derivative of the covariance operator is linked to traditional helioseismic holography.Denoting potential additional dependence of I α,β on the unknown parameters q through H α and H β by superscripts, in terms of conventional holography the Newton step ( 18) is a regularized solution to where K q α,β are the sensitivity kernel of traditional holography, see (20).We will discuss the sensitivity kernels in more detail in Section 7.2.
In traditional helioseismic holography, one has freedom in the choice of the pupils and back-propagators.For example, the pupils can be chosen such that the hologram intensity becomes sensitive to specific flow components [59].As a further example, Porter-Bojarski holograms, introduced to the field of helioseismology in [60,61], make use of the normal derivative at the surface in addition to the Dirichlet data.In contrast, the backward propagators in iterative helioseismic holography are determined by the wave equation, and the image is improved by iteration.
While many techniques in helioseismology including traditional helioseismic holography are limited to linear scenarios, iterative holography naturally allows to tackle nonlinear problems.
Among the commonly used imaging techniques in local helioseismology, holography is the only method that uses the complete cross-correlation data.As already discussed in Section 5.1, these data are used only in an implicit manner without the (usually infeasible) requirement of computing or storing the cross-correlation data explicitly.Whereas traditional helioseismic holography only provides feature maps [16], iterative helioseismic holography additionally allows to retrieve quantitative information.

Kernels and resolution
In the following, we compute the sensitivity kernels for the normal equation as defined in eq. ( 19) which are set up explicitly in our current implementation.It is important to note that sensitivity kernels are typically 3D × 3D operators and should be avoided in computations.Nevertheless, in the spherically symmetric case or two-dimensional medium, computation becomes feasible.Therefore, for the purpose of this paper, we can compute the sensitivity kernels in each iteration and do not study more sophisticated approaches.These kernels are infinitely smooth for smooth coefficients, but they are well localized.It turns out that the width of these kernels is of the order of the classical resolution limit of half a wavelength.This provides an upper bound on the achievable resolution.For simplicity, we will assume a spherically symmetric background without a flow field.
(i) Inversion for source strength: It follows from Theorem 8 that where the multiplication operator is defined in Lemma 5, so with the sensitivity kernel The real part comes from the fact that the source strength has to be a real parameter, it is the adjoint of the embedding of a vector space of real-valued functions into the corresponding vector space of complex-valued functions.The source forward-backward kernel takes the form: where D n is the integral kernel of Γ n in (23).The sensitivity kernel becomes |K α S ,α S | 2 and is therefore non-negative.Furthermore, there are almost no sidelobes after averaging over frequency.
(ii) Inversion for scalar parameters q ∈ {ρ, c, γ}: The operators ∂ q v and ∂ q A for q ∈ {ρ, c, γ} are computed in equations ( 31) and (32).For a flow-free background medium, we have ∂ q A = 0 for all scalar parameters q.It follows from Theorem 8 that In particular, for q ∈ {c, γ}, we have ∂ q v = M(g 0 q ), and the kernel takes the form (iii) Inversion for mass-conserved flow field u: The flow field sensitivity kernel takes the form for i, j ∈ {r, θ, ϕ} where F α u ,α u = Γ Γ H α u (z, x)D n (z, z ′ )H α u (z ′ , y) dz dz ′ and analogue the further kernels.It is remarkable that we are encountering a gradient in the local correlation (∂ u L u ) * .In Chapter 4 of [59], enhancements were accomplished by calculating the difference between two holograms which are spaced apart by half the local wavelength.This can be understood as an approximation to the gradient in the target direction and is therefore naturally incorporated in our framework.
In Figures 2 and 3, we present the sound speed kernel for a uniform and a solar-like radially stratified medium for four different target positions.The kernels are computed for spherical harmonic degrees 0 ≤ l < 100 and averaged over 100 evenly spaced frequencies between 2.75-3.25 mHz.Since there are no strong ghost images on the backside, we show only half of the geometry.The sound speed kernels are very sharp near the target location.Therefore, we can expect the holograms to catch the main features of the image.It is important  (38) in a spherically stratified solar-like background medium and spherical harmonics degrees 0 ≤ l < 100 for four different target positions.We have averaged the sensitivity kernels over 100 frequencies in the frequency regime 2.75-3.25 mHz and normalized with K(x, x) at the target location x.For better comparisons, we have multiplied the sensitivity kernels with the sound speed.
to highlight that the kernels maintain their sharpness even in deep regions within the interior.In addition, this result holds true for a radial stratification similar to that of the Sun.We observe similar behavior for the sensitivity kernels for wave damping, density, source strength, and the components of the flow field.Similar to the sensitivity kernels for the source strength, it is important to note that there are only small visible sidelobes in the sensitivity kernels for sound speed perturbations.This is an additional advantage compared to traditional techniques used in helioseismology.
Figure 4 provides a comparison of the width of the sensitivity kernels for the normal equation and the local half wavelength λ/2.Note that both are of similar size in all cases, .The sound speed sensitivity kernels in a spherically stratified background medium and the l range is 0 ≤ l < 100.In the top panels we present the kernels for a uniform medium with c 0 = 200 km/s (as in Fig. 2), and in the second line the kernels for a solar-like medium (as in Fig. 3).In the first column, we show the kernels in the radial direction, and in the second column the kernels in the angular direction.We have averaged the sensitivity kernels over 100 frequencies in the frequency regime 2.75-3.25 mHz.Furthermore, we compare the width of the sensitivity kernels to the classical resolution limit of λ/2.For better comparisons, we have multiplied the sensitivity kernels with the sound speed.and similar results hold true in angular direction.Therefore, we can expect a resolution of (at least) λ/2.However, in the case of a solar-like stratification, the sensitivity kernels are increasing close to the solar surface.
In helioseismology, and particularly in helioseismic holography, a common issue is the indistinguishability of various sources of perturbations, which complicates the interpretation of seismic data.The design of a holographic back-propagator holds the promise of separating different perturbations.In Figure 5, we present the sensitivity kernels for a perturbation in sound speed, a perturbation in damping, and the cross-kernel in a uniform two-dimensional medium.We show the kernels in a region around the target location.Note that the sound speed kernel is on one scale bigger than the damping kernel and the cross-kernel.Furthermore, the cross-kernel exhibits a different shape with positive and negative maxima around the target location.Therefore, we expect that iterative holography can separate different perturbations in the background medium.

Inversions
In this section, we analyze the performance of iterative holography.The geometry is meshed with a resolution of 10 internal points per local wavelength.Furthermore, we impose a Sommerfeld boundary condition throughout the inversions.
Throughout the following inversions, we employ a L 2 -term as the penalty term and introduce a non-negativity constraint for both sound speed and source strength.The regularization parameter is determined by a power law: α n = α 0 • 0.9 n , where α 0 represents the maximal eigenvalue of the first iteration.The stopping criterion for the inversions is a version of the discrepancy principle for the normal equation, with the noise level determined by the trace of the covariance operator C 4 .In more advanced inversions, stopping rules may be investigated in the hologram space.We set a limit of at most 50 inner conjugate gradient steps per Newton step.Furthermore, we opt for a spatial resolution of 7 grid points per local wavelength.

Holographic image for source perturbation
We have performed a numerical test for a uniform, flow-free two-dimensional medium with source region [0.5, 0.7] 2 and 100 uniformly sampled receivers located on ∂B(0, 1) (see Figure 6).We choose a constant sound speed c = 350 km/s, which corresponds to the solar sound speed at ≈ 0.38R ⊙ .The frequency is fixed to be ω/2π = 3 mHz, which corresponds to the solar 5-minute oscillations.In the case of uniform medium, the differential equation simplifies to a Helmholtz equation, such that the Green's function is analytically known (see Appendix B).Note that Lindsey-Braun holography (H α = H β = G) provides sharp feature maps in the case of small wave damping.For stronger wave damping, the quality of these feature maps deteriorates rapidly.Even without damping, the feature maps are not quantitative at all.

Source strength inversion
Due to its linear nature, inversion for source strength is the simplest case.Therefore, it is in general possible to work with a much finer grid than in the case of parameter identification problems.We add a strong perturbation in the source region [0, 0.5 R ⊙ ] 2 .The inversion results at 3 mHz are shown in the first row of Figure 7 for 10000 realizations.Note that even very deep source terms can be inverted using only one frequency.The reconstructions exhibit a remarkable quality, strongly improving the results by traditional Lindsey-Braun holography (see Figure 6).

Parameter identification
We add a perturbation in the quadratic region [0.5R ⊙ , 0.7R ⊙ ] 2 .Furthermore, we choose 100 evenly spaced frequencies in the frequency range of 2.75 − 3.25 mHz and assume 1000 realizations for each frequency.Note that in helioseismology we have many more frequencies available.
The inversions are shown in the second row of Figure 7.The Newton iteration was stopped after 15 iterations.The resolution of the reconstruction is again below the classical limit of half a wavelength.We observed qualitatively similar results in the inversions for wave damping and density.
The total number of Dopplergrams is given by N ω × N obs , where N ω is the number of frequencies and N obs the number of realizations for each frequency.Note that the total size of Doppler data is fixed by the observation time.We observe that a larger number of frequencies leads to better reconstructions.On the other hand, the computational costs scale roughly linearly with the number of frequencies.This becomes particularly important for large-scale forward problems like for the Sun.Therefore, the choice of N ω often is a trade-off between quality of reconstructions and computation time.

Flow fields
The inversion is performed in a solar-like three-dimensional medium.The example flow field is computed by u = curl ψ, where ψ is a stream function.This guarantees conservation of mass and axisymmetry of the flow field.The stream function is chosen similar to models of meridional circulation profiles in the Sun [62].In the inversion process, we guarantee conservation of mass through Lagrange multipliers, as discussed in Appendix C. The inversion for a symmetric flow field is presented in Figure 8.We inverted with 100 evenly spaced frequencies between 2.75-3.25 mHz and assumed 1500 realizations for each frequency.Since meridional flows are a small perturbation, the iteration is stopped after one iteration.Because of the symmetry, we show only one half-space of the flow field.Besides the strength of the flow field at larger depths, there is no difference visible in the eye-norm.

Conclusions
We have developed a theoretical framework for quantitative passive imaging problems in helioseismology.It shows that traditional holography can be interpreted as an adjoint imaging method.Holographic back-propagation can be seen as part of the adjoint of the Fréchet derivative of the forward operator mapping physical parameters to the covariance operator of the observations.In contrast to traditional holography, the backward propagators are uniquely determined by the wave equation, and the holograms can be improved by iteration rather than clever choices of back-propagators.Iterative helioseismic holography surpasses traditional helioseismic techniques by the quantitative nature of its imaging capabilities and its ability to solve nonlinear problems.
We have demonstrated the performance of iterated holography in inversions for the right hand side of wave equation (source strength), parameters of the zeroth order term (sound speed, absorption) and of the first order term (flows).In all three cases, we have achieved reconstructions with a resolution of slightly less than half of the local wave-length by the iteratively regularized Gauss-Newton method, even for strong realization noise.This is well below the spatial resolution of traditional time-distance helioseismology (see [15]).
Inversions in other more challenging solar setups and for real solar oscillation data are planned as future work and will be presented elsewhere.
In view of the huge size of solar oscillation data, the main bottleneck that prevents the immediate application of iterative holography to interesting large-scale problems in helioseismology is computational complexity.The results of this paper encourage further algorithmic research on iterative regularization methods tailored to passive imaging problems, e.g., by more efficient treatments of sensitivity kernels and Green's functions.
An interesting feature of correlations of Gaussian fields is the structure of the realization noise as described in Section 5.3.A thorough mathematical treatment will require further investigation concerning appropriate stopping rules, consistency, and convergence rates as the sample size tends to infinity.
where γ 0 /2π = 4.29 µHz and ω 0 /2π = 3 mHz.We extend the computational boundary by 500 km above the solar surface (compare with the density scale height of H = 105 km) and apply the radiation boundary condition "Atmo Non Local" (see [52]), assuming an exponential decay of density and constant sound speed in the solar atmosphere.In a spherical symmetric background, we can decompose the Green's function into spherical harmonics: where the Y lm are spherical harmonics.The functions G l (r 1 , r 2 ) satisfy a one-dimensional differential equation and are computed with NGsolve [65], [66].The computation of the Green's function is usually expensive, as the stiffness matrix has to be inverted.The two-step algorithm of [67] allows us to obtain the full modal Green's function from two computations only.Furthermore, this expansion allows us to use a low-rank approximation for the Green's function.

Appendix B. Green's function in uniform medium
We perform numerical toy examples in uniform flow-free two-dimensional and threedimensional background mediums and consider a Sommerfeld boundary condition.The differential equation ( 25) reduces to the Helmholtz equation where k is constant.In this setting, the Green's function is well known (e.g.[31]): for |x| ≥ |y|.Here, J n , h 1 n , j n denote the Bessel function, spherical Hankel function, and spherical Bessel function.Moreover, θ x,y denotes the angular distance between x and y.Furthermore, this basis transformation allows a natural implementation of the singularity.We use this expansion in order to use low-rank approximations for the Green's function.

Proposition 3 .
Let Ω ⊂ R d be open and non-empty.Then there exists a unique bounded linear operator

Figure 1 .
Figure 1.The left panel shows the sound speed and density obtained from the Solar Model S [54] in the solar core, the convection zone (CZ), and the radiation zone (RZ).The right panel shows the potential close to the surface for ω/2π = 3 mHz.
admissible parameters c, ρ, γ, u, S containing some reference parameters c ref , ρ ref , γ ref , u ref , S ref such that the following holds true:

Figure 2 .
Figure2.The sound speed sensitivity kernel K(x, •) in the r − θ-plane as defined in(38) for a three-dimensional uniform medium with c 0 = 200 km/s and the l-range is 0 ≤ l < 100 for four different target positions x.We have averaged the sensitivity kernels over 100 frequencies in the frequency regime 2.75-3.25 mHz and normalized with K(x, x) at the target location x.

Figure 3 .
Figure3.The sound speed sensitivity kernel K(x, •) in the r − θ-plane as defined in(38) in a spherically stratified solar-like background medium and spherical harmonics degrees 0 ≤ l < 100 for four different target positions.We have averaged the sensitivity kernels over 100 frequencies in the frequency regime 2.75-3.25 mHz and normalized with K(x, x) at the target location x.For better comparisons, we have multiplied the sensitivity kernels with the sound speed.

Figure 4
Figure 4.The sound speed sensitivity kernels in a spherically stratified background medium and the l range is 0 ≤ l < 100.In the top panels we present the kernels for a uniform medium with c 0 = 200 km/s (as in Fig.2), and in the second line the kernels for a solar-like medium (as in Fig.3).In the first column, we show the kernels in the radial direction, and in the second column the kernels in the angular direction.We have averaged the sensitivity kernels over 100 frequencies in the frequency regime 2.75-3.25 mHz.Furthermore, we compare the width of the sensitivity kernels to the classical resolution limit of λ/2.For better comparisons, we have multiplied the sensitivity kernels with the sound speed.

Figure 5 .
Figure 5. Matrix-valued sensitivity kernel for joint inversion for sound speed c and damping γ.The left two panels exhibit the diagonal entries, and the right panel the cross-kernel in a uniform two-dimensional medium in a rectangular box of [0.2 R ⊙ , 0.4 R ⊙ ] 2 .The target location is indexed by a red cross.The kernels are normalized by the maximal value of the sound-speed kernel.

Figure 6 .
Figure 6.Lindsey-Braun holographic image intensities in uniform two-dimensional medium (Helmholtz equation) for different degrees of damping with 100 equidistant receivers on ∂B(0, 1).The wave number is such that k 2 = ω 2 (1 + iγ)/c 2 with constant sound speed c = 350 km/s and ω/2π = 3 mHz corresponding to a wavelength of ≈ 0.17 R ⊙ .Note the different scalings of the color maps illustrating the non-quantitative nature of Lindsey-Braun holography.

Figure 7 .
Figure 7. Inversions in two-dimensional uniform background with sound speed c = 350 km/s.In the first row, we present inversions for the source strength at 3 mHz and 10000 realizations.In the second row, we present inversions for the sound speed for 100 frequencies evenly spaced in the band 2.75-3.25 mHz and 1000 realizations.The black cross indicates λ/2.

Figure 8 .
Figure 8. Inversion for the flow field in solar-like three-dimensional background medium in the r − θ-plane.The inversion is performed with 100 evenly spaced frequencies between 2.75-3.25 mHz and 1500 realizations.Due to the symmetry, we show only one half-space of the flow field.