Focusing and Green’s function retrieval in three-dimensional inverse scattering revisited: A single-sided Marchenko integral for the full wave ﬁeld

The Marchenko integral, key to inverse scattering problems across many disciplines, is a long-standing equation that relates single-sided reﬂection data and Green’s functions for virtual source locations inside of an inaccessible, one-dimensional volume. The concept was later expanded to two and three dimensions, yielding important advances in imaging complex media, particularly in the context of geophysics. However, this expansion is based on a set of coupled Marchenko equations which requires up and down decomposition of the wave ﬁelds at both the level of the measurement surface and the level of the virtual source of the desired Green’s function. The underlying theory implies that the recently developed Marchenko relations, while enabling novel applications, carry intrinsic limitations. For example, this scheme cannot incorporate evanescent or refracted waves, and in turn practical implementations must discard data to meet such requirements. We present a derivation that circumvents these limitations, thereby yielding a Marchenko integral akin to those in recent advances, but that is more general than previously assumed. We set up a wave equation based framework to describe the physical concept of focusing functions by introducing homogeneous Green’s functions of the second kind. Based on this, we derive integral representations for both closed and open boundary volumes. Owing to our perspective on the integral formalism, we present an inverse scattering approach for retrieving Green’s functions from single-sided reﬂection data—with the same practical applicability of recent methods but without any limitations due to one-way decomposition. Finally, we illustrate the capability of the Marchenko method to obtain the full wave ﬁeld, including evanescent and refracted waves, within an unknown scattering medium by means of a numerical example.


I. INTRODUCTION
Imaging the interior of an object that is only accessible at its boundary is a key problem in many fields, such as seismology [1,2], helioseismology [3], quantum mechanics [4], medical imaging [5][6][7][8][9][10], and nondestructive testing [11,12].Imaging methods rely on acoustic, seismic [13], or electromagnetic [14] waves to probe the interior of objects with sources and receivers located on their boundaries.The objective of inverse scattering theory [15] is the retrieval of the physical characteristics of the medium from measuring its remotely observed scattering response.
The Marchenko integral is an elementary equation in onedimensional inverse scattering theory [16][17][18].While the medium parameters, e.g., the scattering potential, can be directly inferred from the Marchenko equation in one dimension [19], Broggini et al. [20,21] studied the Marchenko integral's capability to produce Green's functions for virtual sources inside of an inaccessible medium.In contrast to popular interferometric methods for Green's function re-the theoretical advancement, is the fact that those representations can be reliably translated into practical approaches to retrieve the Green's functions of real unknown media from single-sided reflection data [30].While there are certain Marchenko-based primary estimation methods [31,32] that require no a priori knowledge of the medium, solving the three-dimensional Marchenko scheme for in-volume Green's functions requires knowledge of a background mediumparameter model, i.e., a wave-speed model that allows for approximating the direct arrivals of the intended Green's functions-but that does not contain information on the unknown scatterers within the medium.Furthermore, the solution is tied to certain causality arguments that generally hold for one dimension but become less general in higher dimensions [26].Hence, complicated models with spatially varying, strong-contrast medium perturbations can pose an issue for the Marchenko method (e.g., Refs.[33,34]).Nonetheless, the Marchenko integral proved to be a valuable extension of existing Green's function retrieval methods [35,36] and is becoming widely used for geophysical applications (e.g., Refs.[30,[37][38][39]).
In this paper, we present a more general, alternative strategy to deriving Marchenko-type integral relations for Green's function retrieval from remote, single-sided scattering data.We start by introducing homogeneous Green's functions of the second kind, which are an extension of the concept of conventional homogeneous Green's functions.Using reciprocity, we obtain integral representations of these fields for both closed and open boundary systems.The open boundary representation is similar to the previously introduced Marchenko equations [26], but it is obtained without the need for (i) defining an auxiliary, truncated medium wave state or (ii) imposing wave-field decomposition within the medium at the location of the desired Green's function.Hence, the Marchenko integral we propose here can be used to obtain the full Green's function-with wave components propagating in all directions, i.e., including the medium's evanescent and refracted field response.Finally, relying on constraints used by previous approaches, we present a practical scheme to solve our single-sided Marchenko equation for an unknown scattering medium's Green's function.We illustrate our findings with a numerical example.

II. INTEGRAL REPRESENTATIONS FOR FOCUSING AND GREEN'S FUNCTIONS
In the following section, we present a derivation for integrals that relate so-called focusing functions and Green's functions.While previous derivations [26] do not include evanescent and refracted waves, our approach is applicable to the full wave field.We start by introducing a partial differential equation (PDE) for focusing functions.Then we use reciprocity to obtain integral representations.

A. The homogeneous Green's function of the second kind
The acoustic wave equation in the frequency domain is given by with the wave operator where u(x, ω) is the acoustic pressure field at location x = (x 1 , x 2 , x 3 ) and frequency ω.The source function is of volume injection rate density and denoted by s(x, ω); i.e., there are no force sources.The medium is defined by density ρ(x) and propagation velocity c(x).The variable i denotes the imaginary unit.Although we consider acoustics in the frame of this paper, the following derivations should be applicable for other PDEs as well (e.g., Refs.[40,41]).A Green's function g(x, ω; x f ) is a wave field that obeys It is the medium's response to a filtered impulse source at x f .Adding the complex conjugate of Eq. ( 3) to Eq. ( 3) delivers the so-called homogeneous Green's function, i.e., a sourcefree superposition of Green's functions, according to In the context of this paper, we refer to this equation as the homogeneous Green's function of the first kind.The star denotes complex conjugation in the frequency domain, which is similar to time reversal in the time domain.The Green's function g(x, ω; x f ) is a causal wave field, i.e., the wave propagates after the source triggering at time t = 0. Furthermore, it is an outgoing wave field, i.e., with respect to a certain volume of interest that contains the source location x f .Thus, g * (x, ω; x f ) in Eq. ( 4) is an incoming wave field for times t < 0. This wave field focuses at the source location and then, according to Eq. ( 4), keeps on propagating as the forward Green's function, i.e., as an outgoing field for times t > 0. At time t = 0, when focusing occurs, both the incoming and outgoing fields coalesce to equal impulsive sources of opposite polarity, thus satisfying Eq. ( 4).As a next step, we rewrite the time-reversed Green's function; i.e., we replace it by where These equations hold for real-valued q(x, ω; x f ), such that q(x, ω; x f ) − q * (x, ω; x f ) = 0. Apart from this restriction, the source field q(x, ω; x f ) is arbitrary up to this point.Its argument x f does not actually represent a necessary dependence but is kept for the sake of consistency.This reformulation allows for replacing the Green's function by a superposition of so-called focusing fields f (x, ω; x f ).We note here that although we refer to f (x, ω; x f ) as focusing fields, these are more general than in previous definitions (e.g., Ref. [26]), since they are tied to the arbitrary source fields q(x, ω; x f ).The Green's and focusing functions are different wave-field realizations but obey the same source function.When modeling numerically, we conventionally propagate the source energy forward in time from t = 0.This representation is different in that the source produces causal and anticausal fields, suggesting the corresponding wave excitation may also be active at t < 0 relative to reference origin time that is associated with the Green's function.Take, for instance, q(x, ω; x f ) = 0. We can then write One way of visualizing the physics behind this choice is to imagine an incoming field, in this case g * (x, ω; x f ) with half-amplitude, that interacts with the source at zero time.The source energy is larger then the field energy, i.e., the incoming field gets absorbed and additionally a new field, in this case −g(x, ω; x f ) with half-amplitude, is created.Note that the focusing function f (x, ω; x f ) is, in contrast to the Green's function, not unique-because q(x, ω; x f ) can be chosen arbitrarily, so long as it satisfies Eq. ( 6).
We may now write the homogeneous Green's function of the second kind according to Mathematically, this is equivalent to Eq. ( 4), but it gives an additional physical insight.Let f (x, ω; x f ) be an incoming field.This field focuses at x f and, afterward, keeps on propagating as the forward Green's function.On the way to the focal point, however, it also produces a scattered field which is not related to the Green's function, namely the outgoing field f * (x, ω; x f ).The choice which of the two, f (x, ω; x f ) or − f * (x, ω; x f ), is the incoming field is, because of the fields' time symmetry, indeed arbitrary.We call the field that satisfies Eq. ( 8) the homogeneous Green's function of the second kind, for it is a source-free field that yields the causal Green's function.The involved fields, however, are not necessarily the same as those in Eq. (4)-this distinction is essential to our approach.
While we did not specify the source q(x, ω; x f ) in more detail, this subsection introduced the general idea of the homogeneous Green's function of the second kind.In the next subsection, we discuss the interferometric representation of the partial differential Eq. (8).

B. Integral representations
In this subsection, we derive an integral representation for the wave fields in Eq. ( 8).The derivation is similar to derivations for multidimensional convolution and deconvolution [42][43][44].We start our derivation from the general form of Rayleigh's reciprocity theorem for acoustic waves [45,46], i.e., where the superscripts A and B mark two different wave states.In addition to the previously introduced pressure field u = u(x, ω), we also require the particle velocity field v =  x, ω) and s = s(x, ω) denote sources of volume force density and volume injection rate density, respectively.The compressibility is given by κ = κ (x) = 1/[ρ(x)c 2 (x)] and density by ρ = ρ(x).We are investigating a volume V , bounded by the smooth surface ∂V .
We consider two states with identical compressibility and density in a lossless volume V , and thus the last integral in Eq. ( 9) vanishes.Furthermore, we choose not to have any force sources within the volume, i.e., f A = f B = 0.This yields v = 1/(ρiω)∇u (e.g., Ref. [45]).Inside the volume, we use Note that x r lies in V , while the sources for u B are outside of the volume.The variable x shows the dependency of u B on the source distribution and appearance of s B outside V .This wave-state configuration is sketched in Fig. 1.Inserting these definitions into Eq.(9) delivers where we show dependencies on space but omit those on frequency for brevity.Let the medium outside V be reflection free.Then the Green's function is purely outgoing, while p(x; x) is both incoming and outgoing with respect to V .Hence, we can write p(x; x) = p in (x; x) + p out (x; x), (15) where the superscripts mark incoming and outgoing fields, respectively.Here we are neglecting waves that travel along ∂V .Using a far-field approximation (e.g., Refs.[22,47]) or pseudodifferential operator theory [48], the latter circumventing the need for an approximation, one finds that the terms g∇ p out and p out ∇g are identical, thus canceling each other.On the other hand, the terms g∇ p in and p in ∇g deliver the same outcome but with opposite sign, as p in is an incoming field and g an outgoing field.Therefore, Eq. ( 14) becomes If p in (x; x) happens to be a time-reversed field, i.e., p in * (x; x), this equation still holds when neglecting evanescent waves on the boundary ∂V [49].As a next step, we can insert the fields from the wave Eq. ( 8) into Eq.( 16).This is possible because the respective overall field is source free; i.e., just like our previously defined field p(x; x), it has no sources in V .Furthermore, we replace the variable x by x f in order to emphasize the dependence on the focusing location x f .The effective sources, however, are at x and inject the incoming field.This incoming field is then given by f (x; x f ), and thus we get This is the integral form for a closed boundary.We note here that this is the most general representation relating focusing and Green's functions, potentially having applications of its own, which will be the subject of further research.Ultimately, one is not limited to the case of a closed boundary.This is important because one of our main goals is to retrieve a medium Green's function response from remote, single-sided wave data-without access to enclosing boundaries.To that end, let the volume be bounded by a horizontal interface ∂V 0 , e.g., at x 3 = 0, and a half-sphere ∂V 1 .Setting the radius of the half-sphere to infinity and considering only the subset of focusing functions for which f (x ∈ ∂V 1 ; x f ) = 0 the contribution of the respective surface integral vanishes.Incoming and outgoing fields are now down-and up-going fields, respectively.Let x r and x be immediately below and above ∂V 0 , respectively; i.e., both are very close to the surface but the receiver is still in V and the source outside V .As the medium is reflection free outside V and the down-going field p down (x; x) is as such not propagating along the horizontal boundary ∂V 0 , it follows that the field p(x r ; x) on the left-hand side of Eq. ( 16) becomes purely up-going, i.e., where we assume a downward-pointing x 3 axis.This is also sketched in Fig. 2. Inserting Eq. ( 8) into Eq.( 18) again, we now get with This relation has exactly the same form as the single-sided Green's function representation shown by, e.g., Wapenaar et al. [26].However, given our PDE-based derivation, the wave fields in our equation are significantly more general than previously understood.Because we rely on up and down decomposition only on the surface ∂V 0 , the integral fully accounts for evanescent waves inside the medium V .Furthermore, the concept of a truncated model space becomes unnecessary, so refracted and diving waves are included in the representation.Additional notes regarding a comparison with the traditional Marchenko scheme can be found in the Appendix.
We want to stress that integral Eq. ( 17) is entirely general with respect to the focusing function f (x s ; x f ).If we consider the special case of q(x; x f ) = 0, for instance, Eq. ( 17) represents injecting the half-amplitude, time-reversed Green's function governing the half-amplitude homogeneous Green's function.In fact, it is well known that this time-reversal homogeneous Green's function retrieval works for closed boundaries, but not for open boundaries, as depicted in Eq. ( 19) [35].This representation generally produces nonphysical artifacts.In the next section, we present a way of solving a particular form of Eq. ( 19), i.e., the Marchenko integral, for the Green's and focusing function based on an estimate of the first arrival of the Green's function.We also discuss the approach in the context of the aforementioned artifacts arising from the open boundary representation.

III. SOLVING THE MARCHENKO INTEGRAL FOR THE GREEN'S FUNCTION
In this section, we want to focus on solving the open boundary Eq. ( 19) for the medium's unknown Green's function, based on having an estimate of its first arrival as a priori information.The physical arguments closely follow those of Wapenaar et al. [26]; however, we are directly solving the more general uncoupled Marchenko integral, i.e., without up and down decomposition.For the sake of brevity, we neglect all arguments in the following equations and instead use the discrete matrix-operator form [50] of Eq. ( 19), i.e., where we assume any numerical integration details, such as scaling, to be included within the discrete kernel of the R operator.In this paper, we use the composite rectangle rule for numerical integration.We want to stress that the focusing function f which we investigate in the following is just one of many possible realizations of f (x; x f ) in Eq. ( 8), i.e., for a particular source q(x; x f ).The realization that we are interested in is defined by certain causality arguments, more specifically, that f and g are separable in time-space domain such that f is preceding g.These conditions are numerically imposed by introducing the windowing function [26,29].This function acts as a mask; i.e., is zero for times |t| > t d (x r ; x f ) + and one elsewhere.In this context, t d (x r ; x f ) is the first arrival time for a source at x f and a receiver at x r and being greater zero accounts for the limited bandwidth of the wavelet [51].We postulate that this is the realization of f (x; x f ) that is most compact in the time-space domain.This focusing function has to extend up to the time-reversed first arrival of the Green's function in order to be able to produce a focus at x f but can be assumed to be zero before it.In fact, this bound has to be symmetric in time, since we cannot record anything before having injected energy; thus the up-going, time-reversed focusing function must be zero before the time-reversed first arrival of the Green's function, too.These separability assumptions hold in the reflection regime; however, they break in complex media with significant diffracted energy and laterally localized velocity perturbations that induce caustics on the first arrivals.
When presenting the open boundary integral representation above, we alluded to the fact that artifacts may arise in the retrieval of the Green's function.Such artifacts are well known to arise from open-boundary systems when retrieving Green's functions by means of time reversal [35].It is absolutely key that we solve for the injected field instead of simply injecting, e.g., the time-reversed Green's function as done in the context of time reversal.This important distinction is what makes our Green's function retrieval a full wave-field inverse scattering approach, as opposed to a direct application of conventional time-reversal principles.Hence, we find a focusing field that only produces its own time-reversed copy with opposite sign as well as the desired causal-only Green's function when injected from a single-sided measurement surface, but no noise or artifacts by apparently missing boundary data.This specific focusing function is therefore not only defined by the temporal preconditioning, i.e., the window function , but also by the general form of the integral representation itself, i.e., f for an open boundary is different from f for a closed boundary.The window function is designed such that where g = g d + g m , i.e., the sum of the first arrival g d and the coda g m .Note that we use the subscript d which was originally proposed in the Marchenko context to denote the direct arrival in a truncated medium [26], but we are in fact referring to the first arrival in the actual medium, which not only propagates upward toward the boundary but also in all other directions.This important difference means that, here, g d and also the respective travel times t d (x r ; x f ) include diving and refracted waves from the medium below x f .Because the window operator is a filter in the time domain, it acts as a convolutional operator in the frequency domain [26].Applying this windowing to Eq. ( 21) gives This represents the three-dimensional Marchenko integral.It follows from Eq. ( 19) when choosing the most compact focusing function in time along with consequent causality arguments.Time-reversing this relation yields noting that is symmetric in time.Applying the window operator = I − to Eq. ( 21), we also get Rearranging Eq. ( 24) for f * and inserting it into Eq.( 25), we obtain Inserting this into Eq.( 26) and adding g d , we get This expression delivers the Green's function for a virtual source inside of an inaccessible volume measured at receivers on the surface from an estimate of the first arrival of the Green's function g d , the windowing function , and the single-sided reflection data R.While R can be obtained from a measurement, one needs to have an estimate of the physical medium to approximate g d .Generally, we assume that a smooth, kinematically correct version of the actual model is sufficient.Similarly, we require such a model to build the windowing operator , i.e., to find the first arrival times, e.g., via ray tracing.Rather than solving for either g or f , we can also solve for such that both f and g follow from the respective filtering in the time domain, i.e., g = g d + b (30) and If the operator norm of R * R is smaller than one, we can use a Neumann series to estimate the inverse in Eqs. ( 27), (28), and (29), i.e., Regarding the accordingly gained infinite series for g, i.e., one can indeed derive the same result from the coupled Marchenko equations [26].However, the previous result assumes up and down decomposition inside the volume, which means that laterally propagating waves and those which are evanescent in the vicinity of x f are excluded.Furthermore, the previous g d is defined in a truncated medium and thus excludes refracted or diving waves that are often present in real media.We derived the solution without the need of a truncated FIG. 3. Velocity model.The star denotes the location of an exemplary source for which we investigate its wave propagation through the medium.The dark rectangle marks the volume of interest, which is considered for this propagation study.The two black triangles refer to virtual source locations that are examined in more detail in the Appendix.medium, circumventing up and down decomposition at the focusing level.As a result of this, our newly adjusted Marchenko method is, in principle, able to retrieve the full wave-field response of the medium.We did not prove that the focusing function always exists under the above conditions, i.e., separated from the Green's function and for an open boundary.We suspect, however, that it does as long as the model complexity is moderate.If it exists, the inverse in Eqs. ( 27), (28), and ( 29) exists and we can use the described method to estimate the Green's function.There are, of course, still limitations with regard to, e.g., finite apertures [52], spatial sampling of the integrands [53,54], model complexity [33], band limitation [55], and the accessible background model information.In the next section, we show a numerical example in support of our theoretical findings and illustrate some remaining issues.

IV. NUMERICAL EXAMPLE: FULL WAVE-FIELD RETRIEVAL IN A 2D HETEROGENEOUS MEDIUM
In this section, we investigate a numerical example in the geophysical context.We want to show that the Marchenko integral can be used to retrieve the full wave field, including evanescent and refracted waves.Here, we focus on the single-sided integral representation in Eq. ( 19) and solve it as suggested in Eq. (33).Hence, we rely on wave-field observations at the horizontal upper boundary of a heterogeneous, scattering half-space.In real-life applications, such measurements are typically limited to a finite aperture, i.e., sources and receivers only cover a certain extension of ∂V 0 .We use the two-dimensional model in Fig. 3 with constant mass density ρ = 2000 kg/m 3 .We rename x 1 and x 3 for two dimensions, i.e., x = (x, z).We use 501 equally spaced receivers on the surface and record 501 shots for sources on the same grid.Note that the model is more complex than conventional velocity media used for numerical studies of the Marchenko integral (e.g., Refs.[26,36]); i.e., here we consider a shallow high-velocity layer and relatively large velocity contrasts.
In order to visualize the estimated wave propagation through the volume, we solve the Marchenko integral for the Green's functions of all points within the dark rectangle in Fig. 3 according to Eq. ( 33)-regularly sampled from x = 1192 m to x = 2808 m and from z = 0 m to z = 1600 m at every 4 m both in x and z directions.Then we make use of source-receiver reciprocity and choose the same receiver location x r for all these Green's functions to get the Green's function for a source at the respective surface location x r measured at all x f in the volume under investigation.
We start by using the correct first arrival wave field g d and window operator , which in this particular case we obtain from modeling in the correct velocity model.This is done solely as a proof of concept, i.e., to show that the Marchenko integral performs sufficiently well.Figure 4 presents five snapshots of the wave propagation for a source in the top left corner of the target volume.In the left column, we show g d .As we are using the correct medium, these wave fields actually show the correct Green's function; however, they are masked by as we only use the first arrivals.In the right column, we show the correct and estimated snapshots.The two wave fields are superimposed and visualized in an alternating fashion to facilitate their comparison.These snapshots show a very good match between true and estimated fields.In particular, the estimated wave fields include refracted and evanescent waves, observable, e.g., in the region of the fourth layer; see arrows in Figs.4(d), 4(f), and 4(h).While one can only see a single event in Fig. 4(d), Figs. 4(f), and 4(h) show a separate refracted (vertical arrow) and evanescent (diagonal arrow) wave.Looking more closely into the accuracy of the estimated field, Fig. 5 depicts the absolute error of the true and the estimated wave field at 0.6 s.Note that the color bar is clipped at the same amplitudes as in Fig. 4. Hence, Fig. 5 shows that the fields match well almost everywhere within the medium, but it also reveals regions where the Marchenko field is reconstructed slightly worse.These errors mainly manifest as poorly matched amplitudes rather than unwanted wave-field artifacts.The first of these poorly matched regions, marked by I in Fig. 5, refers to the nearly horizontally traveling wave near the surface.Given the up and down decomposition at the surface ∂V 0 in our derivation, it might be challenging for the Marchenko integral to incorporate the energy of this particular field that arrives at the receivers under such a near-horizontal propagation angle.This would probably improve for wider bandwidth and/or larger aperture data.Furthermore, the misfit is right at the edge of the window function , i.e., right after the first arrival of the Green's function.This makes it hard for the method to retrieve accurate amplitudes.The second poorly matched area, marked by II, refers to a rather steep, deeper event.As this up-going wave front is yet to travel through the fourth, high-velocity layer, its slope, i.e., in terms of horizontal wave number, can be assumed to increase even further on its way up.We surmise that this amplitude mismatch is mainly caused by the limited aperture of the numerical experiment, i.e., missing sources and receivers, especially at x > 4 km.We want to stress that the Neumann expansion was truncated at a constant term for this experiment.Including additional terms may also lead to potentially further improved amplitudes for I and II.Overall, our numerical example supports the claim that our uncoupled Marchenko integral can be used to reproduce the full nonlinear scattering of the Green's function, including evanescent and refracted waves within the medium.There are still limits in the accuracy of the here retrieved wave fields mainly related to the band limitation of the data, the windowing operator, and the limited extension of the measurement surface. in Fig. 6, with increasing degrees of inaccuracy relative to the actual medium.Such models may be obtained by tomographic inversion methods, for instance [56].
While the models are potentially too smooth to produce reflections, they can be used to approximate the first arrival g d and the window operator .Then, the Marchenko integral can be utilized to approximate all orders of scattering inside the medium.Generally speaking, the estimated Green's function combines the kinematic information of the background model in the first arrival with the reflection measurement, seeking to  find a consistent wave field that matches both.Figure 7 shows the estimated wave fields for a propagation time of 0.3 s.At first glance, all three estimated fields show comparable results.While the reflections are not known a priori, we retrieve them through the Marchenko scheme, using the information from the smooth background models.It is, however, easily observed that the estimated wave fields also contain a significant amount of noise.There are artifacts, mostly in form of apparently steep, coherent events.While they are lower in amplitude than the desired signal, they are not negligible.Furthermore, we observe that these artifacts become more pronounced for smoother background models; i.e., the field in Fig. 7(b) is better than that in Fig. 7(f).Figure 8 presents the snapshots at a later propagation time of 0.6 s.Again, all estimated fields appear to be of similar, good quality at first sight.In fact, the wave fields seem to be better for higher propagation times; i.e., there are fewer visible artifacts compared to the earlier time counterparts in Fig. 7. Evanescent and refracted waves are still retrieved, albeit with lower accuracy.Upon closer inspection, however, there are still evident biases.Especially for the smoothest model, Fig. 6(c), the interlaced snapshot, Fig. 8(f), reveals significant phase shifts between true and estimated arrivals.This is a result of an inaccurate g d estimate, which is observable by comparing the first arrivals in the bottom right area of Fig. 8(f) between true and estimated fields.This misfit affects all later reflections and produces significant phase and amplitude errors at all orders of scattering.Both local and coherent artifacts can be inspected in greater detail in the enlarged regions in Figs.8(d) and 8(f).While the wave field is still rather good for the former, it reveals a different shape and several pointlike structures in the latter.In fact, these artifacts also relate to the stability of the Marchenko integral.It is not evident that the inverse in Eqs. ( 27), (28), and (29) should always exist, in particular when considering complex media and/or inaccurate background models, the latter reflecting upon the quality of the window operator .Empirically, we find this inverse to be generally stable for moderately heterogeneous media, e.g., media where velocity increases rather monotonically with depth.In these cases, even strong smoothing of the true model delivers an appropriate and a stable inverse.Furthermore, the Neumann expansion in Eq. ( 33) appears to deliver a convergent series then.For more complex settings, as, e.g., the model in Fig. 3, the inverse can become unstable, in particular as the background models used to obtain become smoother.We also observed, in addition, that using a truncated model as required by the original Marchenko approach [26] leads to even greater instability (see the Appendix).Finally, we point out that we use a finite Neumann series to solve the Marchenko equation in this paper, where the order of the last term can be thought of as playing a regularization role.In fact, the leading-order solution already yields an accurate first guess, while remaining stable even for relatively inaccurate background models, i.e., For the results shown in this paper, we use the tenth-order truncated series.This order is chosen for its accuracy seems adequate and the computational cost reasonable, while still delivering sensibly regularized results for the investigated smoothed models.These numerical results support our hypothesis that the Marchenko integral is generally valid for the full wave field.The quality of the reconstructed wave fields is shown to depend on the quality of the required a priori model in terms of producing an accurate first arrival estimate.While we can use smooth medium-parameter estimates, they can introduce artifacts and phase shifts in the resulting wave fields.We found these effects to be strong for highly complex media and addressing a better practical scheme for Green's function retrieval in such cases is the topic of current research.However, with the current practical scheme, the Marchenko integral can generally be used to obtain a rather reliable approximation of the entire wave field within a volume of interest for a wide range of medium configurations.

V. DISCUSSION
We present a derivation for the Marchenko integral which proves that the equation is more general than previously assumed.Introducing the concept of the homogeneous Green's function of the second kind is the key point in this derivation.While previous versions of Marchenko representations target very particular choices for the focusing function, e.g., defined by means of the transmission operator of an auxiliary, truncated medium [26], we present a generalization of focusing functions, which encompasses previous choices but accommodates other approaches.Using conventional reciprocity theorems, one can easily obtain integral representations that relate the Green's and focusing functions to observed reflection data.These data can be obtained either on a closed or an open boundary, and for either case our respective Marchenkolike formalism is well defined.Furthermore, we present a general strategy for solving the uncoupled Marchenko integral to infer the medium's response from open boundary observations, i.e., for a single-sided reflection experiment.The physical arguments that lead this solution are generally equivalent to those of Wapenaar et al. [26].To obtain the Green's function from Eq. ( 19), we rely on the special realization of a focusing function that is most compact in time space.This focusing function has the benefit of being separable from the Green's function in time, thus allowing for a solution of the Marchenko integral based on an estimate of the first arrival g d .Additionally, the method circumvents artifacts that are conventionally introduced by open boundary integral representations, delivering, in principle, an unbiased Green's function estimate.Kiraz et al. [57] recently presented a heuristic, iterative scheme to solve the closed boundary integral Eq. ( 17) using the very same, time-compact, physical realization of a focusing function.In this case, injecting f (x; x f ) into the medium delivers a wave field that, when adding its complex conjugate, equals the homogeneous Green's function of the first kind.Our generalized framework for focusing functions might be useful for future studies, e.g., directly involving the partial differential Eq. ( 8).In particular, there might be possibilities of including focusing functions in other inverse scattering approaches without explicitly relying on the Marchenko equation, such as full waveform inversion [2] or the contrast source method [14].
Our numerical studies show that the Marchenko integral can be used to obtain an accurate approximation of the full wave-field Green's function from only an estimate of its first arrival and single-sided reflection data.The necessary a priori estimate of the first arrival can be based on a reference model, i.e., a smooth estimate of the actual model.We show how, under this theory, the corresponding adaptions to the existing Marchenko workflow produce reliable wave fields for relatively complex models-with the key addition of retrieving evanescent and refracted fields within an unknown scattering medium.Furthermore, we present the estimated wave fields for an entire volume, allowing for a more thorough, spatially dependent analysis of propagation effects.Finally, we discuss the impact of the accuracy of the background model used to approximate g d .For complex media, a poor estimate of the first arrival from an inaccurate reference model can produce significant local artifacts and phase shifts in the recovered fields, but it still allows for a good approximation of the global, scattered field.When in doubt about the quality of the reference medium in achieving sufficiently accurate first arrival estimates, we suggest using a first-order truncated Neumann series, generally allowing for a stable yet reasonably accurate estimate of the Green's function.Alternatively, one may solve Eqs. ( 27), (28), and (29) directly using a numerical solver such as LSQR [58,59].
As mentioned earlier, there are variations of the conventional Marchenko scheme that can be used for primary estimation; i.e., multiple reflections can be filtered out [31,32].Such primary estimation schemes can be applied without the need for a parameter, e.g., wave-speed, model, making them rather attractive for processing wave-field data.It is yet to be investigated how our findings can be linked to these methods.While we only consider acoustic waves in this paper, some studies already investigated the conventional Marchenko method for elastic waves [60][61][62].This extension is valuable for certain data applications, e.g., seismic imaging or medical elastography, and it will be a topic of future research to connect our insights with these studies.
We note that the coupled Marchenko system [26] remains very useful.In the example of geophysical applications, there is reason to directly estimate the decomposed up-and downgoing Green's functions at a certain level in the medium.This allows for so-called redatuming and target-oriented imaging, independent of model perturbations above the redatuming level (e.g., Refs.[63,64]).The implicit model truncation, however, that is inherent to the coupled Marchenko approach might be an issue in complex media.In these cases, it might be beneficial to solve the uncoupled Marchenko integral first (see Sec. III), and then decompose the estimated Green's functions afterward if desired or, alternatively, adjust the coupled Marchenko equations to match our scheme.There are also applications that can make use of the full Marchenkoestimated wave field inside a volume of interest.Such a wave field might, for instance, be used to estimate the scattering potential, i.e., the perturbations of the medium that are missing in the reference model and induce the scattered wave field [65,66].It will be a topic of future research to see what other ways there are for the Marchenko approach to add value to related inverse scattering and imaging schemes and, in this regard, for it to be applied not only in geophysical but also in, e.g., medical applications and nondestructive testing.

VI. CONCLUSION
We introduce the homogeneous Green's function of the second kind, delivering a framework for focusing functions that is substantially more general than in previous Marchenkorelated applications.Based on the resulting partial differential equation, we can construct integral representations for the involved wave fields.The single-sided representation is identical in its form to the previously derived three-dimensional Marchenko relations; however, our derivation imposes significantly fewer limitations on the retrievable wave fields, while accommodating also for closed boundary representations.As such, we find that the Marchenko equation is therefore more general than previously assumed and can indeed be used to obtain the full wave-field response from an unknown scattering medium including evanescent and refracted waves.For practical Green's function retrieval, we present a direct solution of the uncoupled Marchenko integral that follows certain causality assumptions-in a manner analogous to the current approach for the coupled Marchenko representation.It is only in this step that we make use of a particular realization of the focusing function that is purposefully chosen to be separated from the Green's function in time.Lastly, we show numerical examples that illustrate the Marchenko integral's capability of estimating the entire, full field Green's function from an estimate of the first arrival and single-sided reflection data.This paves the way for more complex data applications and a potentially broader usage in related imaging sciences.FIG. 10.Similar to Fig. 9, but for the virtual source location x f marked by the upward-pointing triangle in Fig. 3. a priori information and not an inverse transmission response.As already pointed out several times throughout the paper, we never consider a truncated medium wave state.Therefore, the required a priori information is the first arrival of the Green's function in the actual medium, including diving and/or refracted waves-no approximations involved.
Our newly adjusted Marchenko scheme can incorporate diving and refracted as well as evanescent waves.While we show results to prove these points for our method above, we want to compare them with the conventional approach here to directly highlight differences.Figures 9 and 10 show time-space domain focusing functions f (x r ; x f ) and Green's functions g(x r ; x f ) for the two virtual source locations marked in Fig. 3, both using our approach and the conventional approach.We use the correct first arrival and direct arrival Green's functions, respectively, where the latter is obtained by modeling in the accordingly truncated, correct medium.For both approaches, we use the same, unfiltered reflection data, i.e., including also, e.g., refracted arrivals.The arrows in Figs.9(b) and 10(b) denote refracted waves that are reconstructed by our scheme but not by the conventional one.The conventional method can not obtain any Green's function contributions before the direct arrival in the truncated medium, neglecting refracted waves and, in general, diving waves, as they obey the same physics.Regarding evanescent waves, we note that the traditional theory excludes them [26] while our new theory includes them.
Last but not least, Figs. 9 and 10 illustrate the instability of the conventional approach.Comparing the focusing functions in Figs.9(a) and 10(a) with those in Figs.9(c) and 10(c), one clearly sees largely increased amplitudes in the latter.This indicates a divergence-related energy growth of the Neumann series for the conventional approach.

FIG. 1 .
FIG. 1. Sketch of the wave-state setup for a closed boundary.The rays indicate involved Green's functions, i.e., from x to x s , from x s to x r , and from x to x r .

FIG. 2 .
FIG. 2. Illustration of the wave-state setup for an open boundary.Rays indicate involved Green's functions.

FIG. 4 .
FIG. 4. Snapshots showing the wave propagation for the source and area marked in Fig. 3.The left column, panels (a) to (i), shows the first arrival.The right column, panels (b) to (j), shows the true wave fields and the Marchenko solutions in an interlaced manner.The bottom colors black and gray specify columns showing the true and the estimated wave fields, respectively, separated by white lines.The first row, panels (a) and (b), is for 0.3 s, panels (c) and (d) are for 0.4 s, panels (e) and (f) for 0.5 s, panels (g) and (h) for 0.6 s, and panels (i) and (j) for 0.7 s.The color bars are clipped at 1% of the overall absolute maximum amplitude.In most practical scenarios, one does not typically have access to the correct first arrival.Therefore, we want to analyze the outcomes for three different approximations of g d based on smooth estimates of the correct wave-speed model, displayed

FIG. 6 .
FIG. 6. Smoothed versions of the model in Fig. 3.The smoothing degree increases linearly from panels (a) to (b) to (c).

FIG. 7 .
FIG. 7. Snapshots showing the wave propagation for the source and area marked in Fig. 3. Left and right columns show, as before, the first arrival and the true/estimated wave field, respectively.Black bottom color marks the true, and gray indicates the estimated wavefield columns.The first row [panels (a) and (b)] is obtained for the lightly smoothed velocity model in Fig. 6(a), the second row [panels (c) and (d)] for the moderately smoothed model in Fig. 6(b), and the last row [panels (e) and (f)] for the considerably smoothed model in Fig. 6(c).The color bars are clipped at 1% of the overall absolute maximum amplitude.All snapshots show the same propagation time, i.e., 0.3 s.

FIG. 8 .
FIG. 8. Similar to Fig. 7, but for a propagation time of 0.6 s.The big green windows in panels (d) and (f) show enlarged regions.Their respective locations are marked by the small green windows.