Classical acoustic wave-field representations consist of volume and boundary integrals, of which the integrands contain specific combinations of Green's functions, source distributions, and wave fields. Using a unified matrix-vector wave equation for different wave phenomena, these representations can be reformulated in terms of Green's matrices, source vectors, and wave-field vectors. The matrix-vector formalism also allows the formulation of representations in which propagator matrices replace the Green's matrices. These propagator matrices, in turn, can be expressed in terms of Marchenko-type focusing functions. An advantage of the representations with propagator matrices and focusing functions is that the boundary integrals in these representations are limited to a single open boundary. This makes these representations a suitable basis for developing advanced inverse scattering, imaging and monitoring methods for wave fields acquired on a single boundary.

## I. INTRODUCTION

The aim of this paper is to give a systematic treatment of different types of wave-field representations (with Green's functions, propagator matrices, and Marchenko-type focusing functions), discuss their mutual relations, and indicate some new applications.

Representations with Green's functions. A Green's function is the response of a medium to an impulsive point source. It is named after George Green, who, in a privately published essay,

^{1}introduced the use of impulse responses in field representations.^{2}Wave-field representations with Green's functions have been formulated, among others, for optics,^{3}acoustics,^{4,5}elastodynamics,^{6–9}and electromagnetics.^{10–12}They find numerous applications in forward modeling problems,^{13–15}inverse source problems,^{12,16}inverse scattering problems,^{5,17–19}imaging,^{20–25}time-reversal acoustics,^{26}and Green's function retrieval from ambient noise.^{27–29}Representations with propagator matrices. In elastodynamic wave theory, a matrix formalism has been introduced to describe the propagation of waves in laterally invariant layered media.

^{30,31}This formalism was refined by Gilbert and Backus,^{32}who coined the name propagator matrix. In essence, a propagator matrix “propagates” a wave field (represented as a vectorial quantity) from one plane in space to another. Using perturbation theory, the connection between wave-field representations with Green's functions and the propagator matrix formalism was discussed.^{33}The propagator matrix has been extended for laterally varying isotropic and anisotropic layered media.^{34,35}Propagation invariants for laterally varying layered media have been introduced, and the propagator matrix concept has been proposed for the modeling of reflection and transmission responses.^{36–39}The propagator matrix has also been used in a seismic imaging method, which accounts for multiple scattering in a model-driven way.^{40}Representations with Marchenko-type focusing functions. Building on a one-dimensional (1D) acoustic autofocusing method, it has been shown that the wave field inside a laterally invariant layered medium can be retrieved with the Marchenko method from the single-sided reflection response at the surface of the medium.

^{41–43}This concept was extended to a three-dimensional (3D) Marchenko wave-field retrieval method for laterally varying media.^{44}Central in the 3D Marchenko method are wave-field representations containing focusing functions. These representations have found applications in imaging methods^{45–47}and inverse source problems,^{48}which account for multiple scattering in a data-driven way.

The setup of this paper is as follows. In Sec. II, we briefly review the matrix-vector wave equation for laterally varying media,^{34,35} generalized for different wave phenomena, and we briefly discuss the concept of the Green's matrix, propagator matrix, and Marchenko-type focusing function. The symmetry properties of the matrix-vector wave equation allow the formulation of unified matrix-vector wave-field reciprocity theorems,^{49,50} which are reviewed in Sec. III. These reciprocity theorems form the basis for a systematic treatment of the different types of wave-field representations mentioned above. Traditionally, a wave-field representation is obtained by replacing one of the states in a reciprocity theorem by a Green's state. In Sec. IV, we follow this approach for the matrix-vector reciprocity theorems. By replacing one of the wave-field vectors by the Green's matrix (and the source vector by a unit source matrix), we obtain wave-field representations with Green's matrices. Analogous to this, in Sec. V, we replace one of the wave-field vectors in the reciprocity theorems by the propagator matrix and, thus, obtain wave-field representations with propagator matrices. In Sec. VI, we discuss a mixed form, obtained by replacing one of the wave-field vectors by the Green's matrix and the other by the propagator matrix. In Sec. VII, we discuss the relation between the propagator matrix and Marchenko-type focusing functions and use this relation to derive wave-field representations with focusing functions. We end with conclusions in Sec. VIII.

## II. THE UNIFIED MATRIX-VECTOR WAVE EQUATION, GREEN'S MATRIX, PROPAGATOR MATRIX, AND FOCUSING FUNCTION

### A. The matrix-vector wave equation

We use a unified matrix-vector wave equation as the basis for the derivations in this paper. In the space-frequency domain, it has the following form:^{32–36}

with

Here, $q(x,\omega )$ is a space- and frequency-dependent $N\xd71$ wave-field vector, where $x$ denotes the Cartesian coordinate vector $(x1,x2,x3)$ (with the positive *x*_{3}-axis pointing downward) and *ω* is the angular frequency. The $N/2\xd71$ sub-vectors, $q1(x,\omega )$ and $q2(x,\omega )$, contain wave-field quantities, which are specified for different wave phenomena in Table I. Operator $\u22023$ stands for the differential operator $\u2202/\u2202x3$. Matrix $A(x,\omega )$ is an *N *×* N* operator matrix; it contains the space- and frequency-dependent anisotropic medium parameters and horizontal differential operators $\u22021$ and $\u22022$. The definitions of this operator matrix for different wave phenomena can be found in many of the references mentioned in Sec. I. The $N\xd71$ vector $d(x,\omega )$ contains the space- and frequency-dependent source functions. A comprehensive overview of the operator matrices and source vectors for the wave phenomena considered in Table I (with some minor modifications) is given in Ref. 51.

. | N
. | $q1$ . | $q2$ . |
---|---|---|---|

Acoustic | 2 | p | v_{3} |

Electromagnetic | 4 | $E0=(E1E2)$ | $H0=(H2\u2212H1)$ |

Elastodynamic | 6 | $v=(v1v2v3)$ | $\u2212\tau 3=\u2212(\tau 13\tau 23\tau 33)$ |

Poroelastodynamic | 8 | $(vs\varphi (v3f\u2212v3s))$ | $(\u2212\tau 3bpf)$ |

Piezoelectric | 10 | $(vH0)$ | $(\u2212\tau 3E0)$ |

Seismoelectric | 12 | $(vs\varphi (v3f\u2212v3s)H0)$ | $(\u2212\tau 3bpfE0)$ |

. | N
. | $q1$ . | $q2$ . |
---|---|---|---|

Acoustic | 2 | p | v_{3} |

Electromagnetic | 4 | $E0=(E1E2)$ | $H0=(H2\u2212H1)$ |

Elastodynamic | 6 | $v=(v1v2v3)$ | $\u2212\tau 3=\u2212(\tau 13\tau 23\tau 33)$ |

Poroelastodynamic | 8 | $(vs\varphi (v3f\u2212v3s))$ | $(\u2212\tau 3bpf)$ |

Piezoelectric | 10 | $(vH0)$ | $(\u2212\tau 3E0)$ |

Seismoelectric | 12 | $(vs\varphi (v3f\u2212v3s)H0)$ | $(\u2212\tau 3bpfE0)$ |

For the wave phenomena considered in Table I, the operator matrix $A$ obeys the symmetry properties

with

where $O$ and $I$ are the $N/2\xd7N/2$ zero and identity matrices, respectively. The superscript *t* denotes transposition (meaning that the matrix is transposed and the operators in the matrix are also transposed with $\u22021t=\u2212\u22021$ and $\u22022t=\u2212\u22022$), The “$*$” denotes complex conjugation, and “$\u2020$” denotes transposition and complex conjugation. In general, the medium parameters in $A$ are complex valued and frequency dependent, accounting for losses. The bar above a quantity means that this quantity is defined in the adjoint medium. Hence, if $A$ is defined in a lossy medium, then $A\xaf$ is defined in an effectual medium and vice versa^{52} (a wave propagating through an effectual medium gains energy). For lossless media, the bar can be dropped. For the wave phenomena considered in Table I, the power-flux density *j* in the *x*_{3}-direction is related to the sub-vectors $q1$ and $q2$ according to

As a special case, we consider acoustic waves in an inhomogeneous stationary medium with complex-valued and frequency-dependent compressibility $\kappa (x,\omega )$ and mass density $\rho kl(x,\omega )$. The latter is defined as a tensor to account for effective anisotropy, for example, as a result of fine layering at the micro scale.^{53} The mass density tensor is symmetric, that is, $\rho kl(x,\omega )=\rho lk(x,\omega )$. We introduce the inverse of the mass density tensor, the specific volume tensor $\u03d1kl(x,\omega )$, via $\u03d1kl\rho lm=\delta km$. Einstein's summation convention holds for repeated subscripts (unless otherwise noted); Latin subscripts run from 1 to 3 and Greek subscripts run from 1 to 2. For the acoustic situation, the vectors and matrix in Eq. (1) are given by

with

where *i* is the imaginary unit and $q(x,\omega )$ and $fl(x,\omega )$ are sources in terms of the volume-injection rate density and external force density, respectively.

Finally, for an isotropic medium, using $\u03d1kl=\rho \u22121\delta kl$ (where *ρ* is the scalar mass density), we obtain for the source vector and operator matrix^{54–58}

### B. The Green's matrix

In the space-time domain, a Green's function is the response to an impulsive point source with the impulse defined as $\delta (t)$. The Fourier transform of $\delta (t)$ equals one; hence, in the space-frequency domain, the Green's function is the response to a point source with unit amplitude for all frequencies. We introduce the *N *×* N* Green's matrix $G(x,xA,\omega )$ for an unbounded, arbitrary inhomogeneous, anisotropic medium as the solution of

where $xA=(x1,A,x2,A,x3,A)$ defines the position of the point source, and $I$ is an *N *×* N* identity matrix. Here, $I$ has a size that is different from that in Eq. (6). For simplicity, we use one notation for differently sized identity matrices (the size always follows from the context). Also, for the zero matrix $O$, we use a single notation for differently sized matrices. Similar to the operator matrix $A$, the Green's matrix is partitioned as

Equation (12) does not have a unique solution. To specify a unique solution, we demand that the time-domain Green's function $G(x,xA,t)$ is causal, hence,

This condition implies that $G$ is outward propagating for $|x\u2212xA|\u2192\u221e$.

The simplest representation involving the Green's matrix is obtained when $q$ and $G$ reside in the same medium throughout space and both are outward propagating for $|x\u2212xA|\u2192\u221e$. Whereas $q(x,\omega )$ is the response to a source distribution $d(x,\omega )$ [Eq. (1)], $G(x,xA,\omega )$ is the response to a point source $I\delta (x\u2212xA)$ for an arbitray source position **x**_{A} [Eq. (12)]. Because Eq. (1) and (12) are linear, a representation for $q(x,\omega )$ follows by applying the superposition principle, according to

where $\mathbb{R}$ is the set of real numbers. This representation is a special case of the more general representations with the Green's matrices, which are derived in a more formal way in Sec. IV.

Next, we discuss the 2 × 2 acoustic Green's matrix as a special case of the *N *×* N* Green's matrix. For this situation, $G$ is partitioned as

Here, the first superscript (*p* or *v*) refers to the observed wave-field quantity at $x$ (the acoustic pressure or vertical component of the particle velocity), whereas the second superscript (*f* or *q*) refers to the source type at $xA$ (the vertical component of the force or volume-injection rate). The unit of a specific element of the Green's matrix is the ratio of the units of the observed wave-field quantity at $x$ and the source quantity at $xA$. For example, $[Gp,f(x,xA,\omega )]=[p]/[f]=m\u22122$. For an isotropic medium, all elements can be expressed in terms of the upper-right element as

Here, $\u22023,A$ stands for the differentiation with respect to the source coordinate $x3,A$. Equations (17) and (19) follow directly from Eqs. (11), (12), and (16). Equation (18) follows from Eq. (17) and a source-receiver reciprocity relation, which is derived in Sec. IV A.

We illustrate $Gp,q$, decomposed into plane waves, for a horizontally layered lossless isotropic medium. To this end, we first define the spatial Fourier transform of a space- and frequency-dependent function $u(x,\omega )$ along the horizontal coordinates $xH=(x1,x2)$ according to

where $s=(s1,s2)$, and *s*_{1} and *s*_{2} are the horizontal slownesses. This transform decomposes function $u(x,\omega )$ at a given depth level *x*_{3} into monochromatic plane wave components. Next, we define the inverse temporal Fourier transform per slowness value $s$ as

where $\u211c$ denotes the real part and *τ* is the intercept time.^{59} We apply these transforms to the Green's function $Gp,q(x,xA,\omega )$ and source function $\delta (x\u2212xA)$, choosing $xA=(0,0,x3,0)$ and setting $s2=0$ (the field is cylindrically symmetric in the considered horizontally layered isotropic medium). We, thus, obtain $Gp,q(s1,x3,x3,0,\tau )$, which is the plane wave response (as a function of *x*_{3} and *τ*) to a source function $\delta (x3\u2212x3,0)\delta (\tau )$.

An example of a horizontally layered medium is shown in Fig. 1(a). The propagation velocities in the layers are indicated by *c _{n}* (with $c=1/\kappa \rho $). The depth level of the source is chosen as $x3,0=0$ m; see Fig. 1(b). The vertical green line in Fig. 1 is the line

*τ*= 0, left of which the field is zero due to the causality condition [similar to that in Eq. (14)]. Figure 1(b) further shows the numerically modelled Green's function $Gp,q(s1,x3,x3,0,\tau )$ as a function of

*x*

_{3}and

*τ*for a single horizontal slowness $s1=1/3500$ s/m (the red and blue arrows indicate the downgoing and upgoing waves, respectively). This is the wave field that would be measured by a series of acoustic pressure receivers, vertically above and below the volume-injection rate source at $x3,0=0$ m. Each trace shows $Gp,q(s1,x3,x3,0,\tau )$ for a specific depth

*x*

_{3}as a function of

*τ*(actually, at each depth, the Green's function has been convolved with a time-symmetric wavelet with a central frequency of 50 Hz to get a nicer display).

The horizontal slowness *s*_{1} is related to the propagation angle *α _{n}* in layer

*n*via $sin\u2009\u2009\alpha n=s1cn$, hence, in layer 1 (with $c1=1600$ m/s), the propagation angle is $\alpha 1=27.2\xb0$. In the thin layer, with $c3=3600$ m/s, we obtain $sin\u2009\u2009\alpha 3=1.03$, meaning that

*α*

_{3}is complex valued. This implies that the wave is evanescent in this layer. Because the layer is thin, the wave tunnels through the layer and continues with a lower amplitude as a downgoing wave in the lower half-space.

### C. The homogeneous Green's matrix

For a lossless medium, a homogeneous Green's function is the superposition of a Green's function and its complex conjugate (or, in the time domain, its time-reversed version). The superposition is chosen such that the source terms of the two functions cancel each other, hence, a homogeneous Green's function obeys a wave equation without a source term.^{19,20} Here, we extend this concept for the matrix-vector wave equation for a medium with losses.

Let $G(x,xA,\omega )$ be, again, the outward propagating solution of Eq. (12) for an unbounded, arbitrary inhomogeneous, anisotropic medium. We introduce the Green's matrix of the adjoint medium, $G\xaf(x,xA,\omega )$, as the outward propagating solution of

Pre- and post-multiplying all terms by $J$ and, subsequently, using Eq. (5) and $JJ=I$ gives

Subtracting the complex conjugate of all terms in Eq. (23) from the corresponding terms in Eq. (12) yields

with

Because $Gh(x,xA,\omega )$ obeys a wave equation without a source term, we call it the homogeneous Green's matrix. According to Eqs. (6), (13), and (25), it is partitioned as

### D. The propagator matrix

We introduce the *N *×* N* propagator matrix $W(x,xA,\omega )$ for an arbitrary inhomogeneous anisotropic medium as the solution of the unified matrix-vector wave equation (1) but without the source vector $d$. Hence,

Similar to the operator matrix $A$ and Green's matrix $G$, the propagator matrix is partitioned as

Equation (27) does not have a unique solution. To specify a unique solution, we impose the boundary condition

where $xH,A$ denotes the horizontal coordinates of $xA$, hence, $xH,A=(x1,A,x2,A)$. Because Eq. (27) is first order in $\u22023$, a single boundary condition suffices. Note that $W(x,xA,\omega )$ only depends on the medium parameters between the depth levels $x3,A$ and *x*_{3}. This is different from the Green's matrix $G(x,xA,\omega )$, which, for an arbitrary inhomogeneous medium, depends on the entire medium [this is easily understood from the illustration in Fig. 1(b)].

The simplest representation involving the propagator matrix is obtained when $q$ and $W$ reside in the same medium in the region between $x3,A$ and *x*_{3} and both have no sources in this region. Whereas $q(x,\omega )$ obeys no boundary conditions in this region, $W(x,xA,\omega )$ collapses to $I\delta (xH\u2212xH,A)$ at the depth level $x3,A$ [Eq. (29)]. Applying the superposition principle again yields

where $\u2202DA$ is the horizontal boundary defined as $x3=x3,A$.^{32–36} Note that $W(x,xA,\omega )$ propagates the field vector $q$ from the depth level $x3,A$ to *x*_{3}, hence, the name “propagator matrix.” This representation is a special case of the more general representations with the propagator matrices, which are derived in a more formal way in Sec. V.

When we replace $q(xA,\omega )$ by a Green's matrix $G(xA,xB,\omega )$ with $xB$ outside the region between $x3,A$ and *x*_{3}, we obtain

This is the simplest relation between the Green's matrix and propagator matrix. It is a special case of the more general representations with the Green's matrices and propagator matrices, which are derived in a more formal way in Sec. VI A. Alternatively, we may replace $G(x,xB,\omega )$ by the homogeneous Green's matrix $Gh(x,xB,\omega )$, where $xB$ may be located anywhere because the homogeneous Green's matrix has no source at $xB$, hence,

This relation will be derived in a more formal way in Sec. VI B.

Next, we discuss the 2 × 2 acoustic propagator matrix as a special case of the *N *×* N* propagator matrix. For this situation, $W$ is partitioned as

The first and second superscripts refer to the wave-field quantities at $x$ and **x**_{A}, respectively (where the superscript *p*, again, stands for the acoustic pressure and *v* is the vertical component of the particle velocity). The unit of a specific element of the propagator matrix is the ratio of the units of the wave-field quantities at $x$ and $xA$. For example, $[Wp,v(x,xA,\omega )]=[p]/[v]=kg\u2009s\u22121\u2009\u2009m\u22122$. For an isotropic medium, the elements can be expressed in terms of the upper-right element as

Equations (34) and (36) follow directly from Eqs. (11), (27), and (33). Equation (35) follows from Eq. (34) and a source-receiver reciprocity relation, which is derived in Sec. V A.

We illustrate the elements $Wp,p$ and $Wp,v$, decomposed into plane waves, for the horizontally layered medium of Fig. 1(a). We choose, again, $xA=(0,0,x3,0)$ and set $s2=0$. Usually, the propagator matrix is considered in the frequency domain, but to facilitate the comparison with the Green's function in Fig. 1(b), we consider the time-domain functions $Wp,p(s1,x3,x3,0,\tau )$ and $Wp,v(s1,x3,x3,0,\tau )$ [obtained via the transforms of Eqs. (20) and (21)], once more, for a single horizontal slowness $s1=1/3500$ s/m; see Figs. 2(a) and 3(a). At $x3=x3,0$ (the horizontal green lines in Figs. 2(a) and 3(a)), the boundary conditions for these functions are $Wp,p(s1,x3,0,x3,0,\tau )=\delta (\tau )$ and $Wp,v(s1,x3,0,x3,0,\tau )=0$, respectively [this follows from applying the transforms of Eqs. (20) and (21) to Eq. (29), with $xH,A=(0,0)$]. Following these functions along the depth coordinate, starting at $x3=x3,0$, we observe an acausal upgoing wave and a causal downgoing wave until we reach the first interface. Here, both events split into upgoing and downgoing waves below the interface. The waves tunnel through the high-velocity layer, split again, and continue with higher amplitudes in the next layer. This illustrates that the evanescent waves may lead to unstable behaviour of the propagator matrix and should be handled with care. Note that $Wp,p(s1,x3,x3,0,\tau )$ and $Wp,v(s1,x3,x3,0,\tau )$ are symmetric and asymmetric in time, respectively. This is best seen in Figs. 2(b) and 3(b), which show the last traces of both elements. The same holds for elements $Wv,v$ and $Wv,p$, which are not shown.

### E. The Marchenko-type focusing function

For an arbitrary inhomogeneous anisotropic medium, the time-domain version of boundary condition (29) reads

This boundary condition has a form similar to the focusing condition for the focusing functions appearing in the multidimensional Marchenko method.^{44} In Sec. VII A, we discuss the general relations between the propagator matrix and Marchenko-type focusing functions. Here, we present a short preview of these relations by considering the acoustic propagator matrix in the horizontally layered lossless isotropic medium of Fig. 1(a). We combine the elements $Wp,p$ and $Wp,v$, decomposed into plane waves (Figs. 2 and 3), as^{60}

where the vertical slowness $s3,0$ is defined as

where *c*_{0} and *ρ*_{0} are the propagation velocity and mass density of the upper half-space $x3\u2264x3,0$, respectively. Due to the symmetric form of $Wp,p$ and the asymmetric form of $Wp,v$, half of the events double in amplitude and the other half of the events cancel. The result is shown in Fig. 4. For the interpretation of this focusing function, we start at the bottom of Fig. 4(a). The blue arrows indicate the upgoing waves, which are tuned such that after the interaction with the tunneling waves in the thin layer and downgoing wave just above the thin layer (indicated by a red arrow), they continue as a single upgoing wave, which finally focuses at the depth level $x3,0$ as a temporal delta function $\delta (\tau )$ and continues as an upgoing wave into the homogeneous upper half-space.

Note that, although we interpret the focusing function here in terms of upgoing and downgoing waves, it is derived from the propagator matrix (which does not rely on up/down decomposition), vertical slowness $s3,0$, and the mass density *ρ*_{0} of the upper half-space. Moreover, the focusing function is defined in the actual medium rather than in a truncated version of the actual medium as is usually the case for the Marchenko-type focusing functions.^{44} Hence, the only assumption is that a real-valued vertical slowness $s3,0$ exists in the upper half-space. Other than that, the focusing function $Fp(s1,x3,x3,0,\tau )$, as defined in Eq. (38), does not require a truncated medium, does not rely on up/down decomposition inside the medium, and does (at least in principle) not break down when the waves become evanescent inside the medium. These properties also hold for the more general version of the Marchenko-type focusing functions defined in Sec. VII.

## III. UNIFIED MATRIX-VECTOR WAVE-FIELD RECIPROCITY THEOREMS

As the basis for the derivation of the general representations with Green's matrices (Sec. IV), propagator matrices (Sec. V), or a combination thereof (Sec. VI), we introduce unified matrix-vector wave-field reciprocity theorems. In general, a wave-field reciprocity theorem interrelates two wave states (sources, wave fields, and medium parameters) in the same spatial domain.^{12} The reciprocity theorems have been formulated for acoustic,^{4} electromagnetic,^{61} elastodynamic,^{62,63} poroelastodynamic,^{64} piezoelectric,^{65,66} and seismoelectric waves.^{67} The matrix-vector equation, discussed in Sec. II, allows a unified formulation of the reciprocity theorems for these different wave phenomena,^{49,50} which extends the theory of the propagation invariants.^{36–39} Using the symmetry properties of operator matrix $A$, formulated in Eqs. (3) and (4), the matrix-vector reciprocity theorems can be derived^{49,50}

and

Here, $D$ denotes a domain enclosed by two infinite horizontal boundaries $\u2202D0$ and $\u2202DM$ at the depth levels $x3,0$ and $x3,M$ with outward pointing normals $n3=\u22121$ and $n3=1$, respectively; see Fig. 5. The subscripts *A* and *B* refer to two independent states. These theorems hold for lossless media and media with losses.^{51} Equation (40) is a convolution-type reciprocity theorem as products like $qAtNqB$ in the frequency domain correspond to convolutions in the time domain (similar to those in Ref. 63). Equation (41) is a correlation-type reciprocity theorem because products such as $qA\u2020KqB$ in the frequency domain correspond to correlations in the time domain (as in Ref. 18).

A special case is obtained when the sources, wave fields, and medium parameters are identical in both states. We may then drop the subscripts *A* and *B*, and Eq. (41) simplifies to

Because $14q\u2020Kq=14(q1\u2020q2+q2\u2020q1)=j$ [see Eq. (7)], Eq. (42) formulates the unified power balance. The term on the left-hand side is the power generated by the sources in $D$. The first term on the right-hand side is the power flux through the boundary $\u2202D0\u222a\u2202DM$ (i.e., the power leaving the domain $D$), and the second term on the right-hand side is the dissipated power in $D$.

## IV. REPRESENTATIONS WITH GREEN'S MATRICES

### A. Symmetry property of the Green's matrix

Before we derive wave-field representations, we first derive a symmetry property of the Green's matrix. To this end, we replace both of the wave-field vectors **q**_{A} and **q**_{B} in the reciprocity theorem (40) by the Green's matrices $G(x,xA,\omega )$ and $G(x,xB,\omega )$, respectively. Accordingly, we replace the source vectors **d**_{A} and **d**_{B} by $I\delta (x\u2212xA)$ and $I\delta (x\u2212xB)$, respectively, where **x**_{A} and **x**_{B} denote the source positions; see Fig. 5. Furthermore, we replace $D$ by $\mathbb{R}3$ such that the boundary integral in Eq. (40) vanishes (the Sommerfeld radiation condition). Both of the Green's matrices are defined in the same medium, hence, $AA=AB$. This implies that the second integral on the right-hand side of Eq. (40) also vanishes. From the remaining integral, we, thus, obtain the following symmetry property of the Green's matrix:

The Green's matrix on the left-hand side is the response to a source at **x**_{A}, observed by a receiver at **x**_{B}. Similarly, the Green's matrix on the right-hand side is the response to a source at **x**_{B}, observed by a receiver at **x**_{A}. Hence, Eq. (43) is a unified source-receiver reciprocity relation.

Using this relation and $JN=\u2212NJ$, for the homogeneous Green's matrix defined in Eq. (25), we find

### B. Representations of the convolution type with the Green's matrix

We derive a representation of the convolution type for the actual wave-field vector $q(x,\omega )$, emitted by the actual source distribution $d(x,\omega )$ in the actual medium; the operator matrix in the actual medium is defined as $A(x,\omega )$. We let state *B* in the reciprocity theorem (40) be this actual state, hence, we drop the subscript *B* from **q**_{B}, **d**_{B}, and $AB$. For state *A*, we choose the Green's state. Therefore, we replace $qA(x,\omega )$ in the reciprocity theorem (40) by $G(x,xA,\omega )$ and $dA(x,\omega )$ by $I\delta (x\u2212xA)$. We keep the subscript *A* in $AA(x,\omega )$ to account for the fact that, in general, this operator matrix is defined in a medium that may be different from the actual medium. We, thus, obtain

where $\chi D(x)$ is the characteristic function for the domain $D$, which is defined as

Using the symmetry property of the Green's matrix, formulated by Eq. (43), we obtain

This is the unified wave-field representation of the convolution type with the Green's matrix. The left-hand side is the wave-field vector $q$ at a specific point **x**_{A}, multiplied with the characteristic function. According to the right-hand side, this field consists of a contribution from the source distribution in $D$ (the first integral), a contribution from the wave field on the boundary of $D$ (the second integral), and a contribution caused by the contrasts between the operator matrices in the Green's and the actual state (the third integral).

Note that Eq. (15) follows as a special case of Eq. (47) if we choose the same medium parameters in both states and replace $D$ by $\mathbb{R}3$ [except that the roles of $x$ and $xA$ are interchanged because we used the symmetry property of Eq. (43) in the derivation of Eq. (47)].

Next, we consider another special case. Again, we choose the same medium parameters in both states, but this time we replace $D$ by the entire half-space below $\u2202D0$, which we choose to be source-free in the actual state. Assuming **x**_{A} lies in the lower half-space, we obtain

Using Eqs. (8), (16), and (18) for the isotropic acoustic situation, for the upper element of $q(xA,\omega )$, Eq. (48) yields

which is the well-known acoustic Kirchhoff-Helmholtz integral. Hence, the representation of Eq. (48) is the unified Kirchhoff-Helmholtz integral. It can be used for forward wave-field extrapolation of the wave-field vector $q(x,\omega )$ from $\u2202D0$ to any point **x**_{A} below $\u2202D0$. Note that the Green's matrix $G(xA,x,\omega )$ depends on the medium parameters of the entire half-space below $\u2202D0$. In Sec. V B, we derive a relation similar to Eq. (48) but with the Green's matrix replaced by the propagator matrix, which depends only on the medium parameters between $x3,0$ and $x3,A$.

Finally, we replace $q(x,\omega )$ by $G(x,xB,\omega )$ in Eq. (48). Because we assumed that the lower half-space is source-free for $q$, we choose **x**_{B} in the upper half-space. We, thus, obtain

This expression shows that the Green's matrix can be composed from a cascade of Green's matrices, assuming a specific order of the depth levels at which the Green's sources and receivers are situated. In Sec. V B, we derive a similar relation for propagator matrices for an arbitrary order of depth levels.

### C. Representations of the correlation type with the Green's matrix

We derive a representation of the correlation type for the actual wave-field vector $q(x,\omega )$, emitted by the actual source distribution $d(x,\omega )$ in the actual medium. As in Sec. IV B, we let the state *B* be this actual state, and state *A* be the Green's state. Hence, making the same substitutions as in Sec. IV B, this time in the reciprocity theorem (41), using Eq. (43) and $N\u22121K=\u2212J$, yields

This is the unified wave-field representation of the correlation type with the Green's matrix.

As a special case, we derive a representation for the homogeneous Green's matrix $Gh$. To this end, for state *B*, we choose the Green's matrix in the actual medium, hence, we replace $q(x,\omega )$ by $G(x,xB,\omega )$ and $d(x,\omega )$ by $I\delta (x\u2212xB)$. For state *A*, we replace the Green's matrix by that in the adjoint of the actual medium, therefore, we replace $G(xA,x,\omega )$ by $G\xaf(xA,x,\omega )$ and $AA$ by $A\xaf$. With this choice, the contrast operator $\u2009\u2009\u2009A\xafA\u2212A=A\xaf\xaf\u2212A$ vanishes. Making these substitutions in Eq. (51), taking **x**_{A} and **x**_{B} both inside of $D$, and using Eq. (25), we obtain

This is the unified homogeneous Green's matrix representation. It finds applications in optical, acoustic, and seismic holography,^{20,25} inverse source problems,^{16} inverse scattering methods,^{19} and Green's function retrieval by cross correlation.^{27–29,70,71} A disadvantage is that the integral is taken along the two boundaries $\u2202D0$ and $\u2202DM$, whereas in many practical situations, the measurements are only available on a single boundary. Using the propagator matrix, in Sec, VI B, we present a single-sided unified homogeneous Green's function representation.

Finally, using Eqs. (16)–(18) for the isotropic acoustic situation, for the upper-right element of $Gh(xA,xB,\omega )$, Eq. (52) yields

where

## V. REPRESENTATIONS WITH PROPAGATOR MATRICES

We follow an approach similar to that used in Sec. IV to derive wave-field representations. However, this time, we replace one of the states in the reciprocity theorems by a propagator state (instead of a Green's state). Hence, this leads to the wave-field representations with the propagator matrices.

### A. Symmetry properties of the propagator matrix

We start with deriving the symmetry properties of the propagator matrix. To this end, we replace the wave-field vectors **q**_{A} and **q**_{B} in the reciprocity theorem (40) by the propagator matrices $W(x,xA,\omega )$ and $W(x,xB,\omega )$, respectively. We define horizontal boundaries $\u2202DA$ and $\u2202DB$, containing the points **x**_{A} and **x**_{B}, respectively; see Fig. 5. We replace $\u2202D0$ and $\u2202DM$ in Eq. (40) by these boundaries (and $D$ by the region enclosed by these boundaries). Note that the boundary condition (29) implies $W(x,xA,\omega )=I\delta (xH\u2212xH,A)$ for $x$ at $\u2202DA$ and $W(x,xB,\omega )=I\delta (xH\u2212xH,B)$ for $x$ at $\u2202DB$. The propagators obey the wave equation (27) without sources, hence, we set **d**_{A} and **d**_{B} to $O$. This implies that the integral on the left-hand side of Eq. (40) vanishes. Both propagator matrices are defined in the same medium, hence, $AA=AB$. This implies that the second integral on the right-hand side vanishes. Evaluating the remaining integral along the boundary $\u2202DA\u222a\u2202DB$ yields (irrespective of the arrangement of $\u2202DA$ and $\u2202DB$)

This is the first unified symmetry relation for the propagator matrix.

A second symmetry relation can be derived from the reciprocity theorem (41). To this end, we replace **q**_{A} and **q**_{B} by $W\xaf(x,xA,\omega )$ and $W(x,xB,\omega )$, respectively. We replace $\u2202D0\u222a\u2202DM$ in Eq. (41) by $\u2202DA\u222a\u2202DB$ (and $D$ by the enclosed region). We set **d**_{A} and **d**_{B} to $O$, therefore, the integral on the left-hand side vanishes. The propagator matrices are defined in mutually adjoint media, hence, $AA=A\xafB$. This implies that the second integral on the right-hand side of Eq. (41) vanishes. From the remaining integral, we, thus, obtain

Note that in this last equation, **x**_{B} and **x**_{A} appear in the same order on the left- and right-hand sides.

### B. Representations of the convolution type with the propagator matrix

We derive a representation of the convolution type for the actual wave-field vector $q(x,\omega )$, emitted by the actual source distribution $d(x,\omega )$ in the actual medium; the operator matrix in the actual medium is defined as $A(x,\omega )$. We let state *B* in reciprocity theorem (40) be this actual state, hence, we drop the subscript *B* from **q**_{B}, **d**_{B}, and $AB$. For state *A*, we choose the propagator state. Therefore, we replace $qA(x,\omega )$ in the reciprocity theorem (40) by $W(x,xA,\omega )$ and $dA(x,\omega )$ by $O$. We keep the subscript *A* in $AA(x,\omega )$ to account for the fact that, in general, this operator matrix is defined in a medium that may be different from the actual medium. We replace $\u2202D0\u222a\u2202DM$ in the reciprocity theorem (40) by $\u2202D0\u222a\u2202DA$, where $\u2202DA$ is the boundary containing **x**_{A}. Here and in the following, $\u2202DA$ is below $\u2202D0$, thus, $x3,A>x3,0$. The domain enclosed by these boundaries is called $DA$; see Fig. 6. Unlike in the classical, decomposition-based derivations of the Marchenko method,^{44,72} the domain $DA$ does not define a truncated medium; in general, the medium below $\u2202DA$ is inhomogeneous. Applying the mentioned substitutions in Eq. (40), using the boundary condition (29) and symmetry relation (55) (with **x**_{B} replaced by $x$), yields

This is the unified wave-field representation of the convolution type with the propagator matrix. Note the analogy with the representation of the convolution type with the Green's matrix, Eq. (47). An important difference is that the boundary integral in Eq. (58) is single sided.

We consider a special case by choosing the same medium parameters in $DA$ in both states and choosing $DA$ to be source free. We obtain

This is the special case that was already presented in Eq. (30) [except that the roles of $x$ and $xA$ are interchanged because we used the symmetry property of Eq. (55) in the derivation of Eq. (58)].

Note the analogy of Eq. (59) with Eq. (48), which contains a Green's matrix instead of the propagator matrix. The propagator matrix in Eq. (59) depends only on the medium parameters between $x3,0$ and $x3,A$. Equation (59) is used in Sec. VII as the basis for deriving the Marchenko-type representations.

Next, we again consider Eq. (58) in which we replace $q(x,\omega )$ by $W(x,xB,\omega )$ and, hence, $d(x,\omega )$ by $O$ and $A$ by $AA$. This implies that the first and third integrals on the right-hand side vanish. Moreover, we replace $\u2202D0$ by $\u2202DC$ at an arbitrarily chosen depth level $x3,C$; see Fig. 6. This yields

This expression shows that the propagator matrix can be composed from a cascade of propagator matrices.^{32–36} It is similar to Eq. (50) with the cascade of the Green's matrices, but unlike in Eq. (50), where $x3,A>x3,0>x3,B$, the arrangement of $x3,A,\u2009\u2009x3,B$, and $x3,C$ in Eq. (60) is arbitrary (because there are no sources in both states).

Finally, we replace **x**_{B} in Eq. (60) by $xA\u2032=(xH,A\u2032,x3,A)$. Applying the boundary condition (29) to the left-hand side, we obtain

This equation defines $W(xA,x,\omega )$ as the inverse of $W(x,xA,\omega )$.^{32–36}

### C. Representations of the correlation type with the propagator matrix

We aim to derive a representation of the correlation type for the actual wave-field vector $q(x,\omega )$, emitted by the actual source distribution $d(x,\omega )$ in the actual medium. Similar to Sec. V B, we let state *B* be this actual state. For state *A*, we choose the adjoint propagator state. Hence, we replace $qA(x,\omega ),\u2009AA(x,\omega )$, and $dA(x,\omega )$ by $W\xaf(x,xA,\omega ),A\xafA(x,\omega )$, and $O$, respectively. Furthermore, we again replace $\u2202D0\u222a\u2202DM$ by $\u2202D0\u222a\u2202DA$. Applying these substitutions in the reciprocity theorem (41), using the boundary condition (29) and symmetry relation (56) (with **x**_{B} replaced by $x$), yields, once more, Eq. (58). Therefore, due to the additional symmetry property of the propagator matrix [Eq. (56)], the correlation-type reciprocity theorem does not lead to a new representation.

## VI. REPRESENTATIONS WITH GREEN'S MATRICES AND PROPAGATOR MATRICES

We follow an approach, which is similar to that in Secs. IV and V, to derive wave-field representations. However, this time, we replace one of the states in the reciprocity theorem of the convolution type by a Green's state and the other state by a propagator state. This leads to wave-field representations with Green's matrices and propagator matrices.

### A. Single-sided Green's matrix representation

In the reciprocity theorem (40), we choose the propagator state for state *A* and the Green's state for state *B*. Hence, in state *A*, we replace $qA(x,\omega )$ and $dA(x,\omega )$ by $W(x,xA,\omega )$ and $O$, respectively. In state *B*, we replace $qB(x,\omega )$ and $dB(x,\omega )$ by $G(x,xB,\omega )$ and $I\delta (x\u2212xB)$, respectively. The medium parameters in states *A* and *B* may be different, thus, we keep the subscripts in $AA(x,\omega )$ and $AB(x,\omega )$. Last but not least, we replace $\u2202D0\u222a\u2202DM$ by $\u2202D0\u222a\u2202DA$ and $D$ by the enclosed region $DA$; see Fig. 6. With these substitutions, using boundary condition (29) and symmetry relation (55), the reciprocity theorem (40) yields

where $\chi DA(x)$ is the characteristic function for domain $DA$, defined similarly as $\chi D(x)$ in Eq. (46). Equation (62) is a representation for $G\u2212W$ (when $xB$ is inside $DA$) in terms of integrals containing $G$ and $W$. When $G$ and $W$ are defined in the same medium, this representation simplifies to

Finally, when $xB$ (the position of the source of the Green's matrix) lies outside of $DA$ (i.e., above $\u2202D0$ or below $\u2202DA$), then the latter expression simplifies to

This is a single-sided representation of the Green's matrix, which uses the propagator matrix. This special case was already presented in Eq. (31) (except that the roles of $x$ and $xA$ are interchanged as we used the symmetry property of Eq. (55) in the derivation of Eq. (64)].

Note the analogy of Eq. (64) with Eq. (50), which contains a Green's matrix instead of the propagator matrix. The propagator matrix in Eq. (64) depends only on the medium parameters between $x3,0$ and $x3,A$. Moreover, unlike the representation of Eq. (50), which only holds for $\u2009x3,B<x3,0$, Eq. (64) also holds for $x3,B>x3,A$.

### B. Single-sided homogeneous Green's matrix representation

We follow the same procedure as in Sec. VI A, except that in state *B*, we choose the homogeneous Green's matrix; hence, we replace $qB(x,\omega )$ and $dB(x,\omega )$ in the reciprocity theorem (40) by $Gh(x,xB,\omega )$ and $O$, respectively. We, thus, obtain

When $Gh$ and $W$ are defined in the same medium, this representation simplifies to

This is a single-sided representation of the homogeneous Green's matrix, which uses the propagator matrix. This special case was already presented in Eq. (32) [except that the roles of $x$ and $xA$ are interchanged as we used the symmetry property of Eq. (55) in the derivation of Eq. (66)]. Note that unlike the Green's matrix representation of Eq. (64), the homogeneous Green's matrix representation of Eq. (66) has no restrictions for the position of **x**_{B} (as there is no source at **x**_{B}).

Equation (66) is the single-sided counterpart of the unified homogeneous Green's matrix representation of Eq. (52). Whereas in Eq. (52) the integral is taken along two boundaries $\u2202D0$ and $\u2202DM$, in Eq. (66) the integral is taken only along $\u2202D0$. This is an important advantage for the practical situations in which the measurements are often available only on a single boundary. In Sec. VI C, we indicate how Eq. (66) can be used in a process called source and receiver redatuming.

Finally, using Eqs. (16), (17), (26), (33), and (35) for the isotropic acoustic situation, assuming real-valued $\rho (x)$, Eq. (66) yields for the upper-right element of $Gh(xA,xB,\omega )$

with $Ghp,q(xA,xB,\omega )$ as defined in Eq. (54).

### C. Source and receiver redatuming

Redatuming is the process of moving sources and/or receivers from the acquisition boundary to positions inside of the medium.^{22,57,73,74} Traditionally, this is done with the Kirchhoff-Helmholtz integrals with the Green's functions defined in a smooth background medium, thus, ignoring multiple scattering. Here, we derive a unified redatuming method, which accounts for multiple scattering, using the single-sided homogeneous Green's matrix representation of Sec. VI B as the starting point.

Let **x**_{S} and **x**_{R} denote the source and receiver coordinates, respectively, at the acquisition boundary $\u2202D0$. Then the Green's matrix $G(xR,xS,\omega )$ represents the medium's response, measured with the sources and receivers at the surface; see Fig. 7(a). Assuming that the medium is lossless, the homogeneous Green's matrix $Gh(xR,xS,\omega )$ is obtained from this Green's matrix via Eq. (25).

First, we discuss a method for source redatuming. We rename some variables in Eq. (66) ($xB\u2192xR,\u2009xA\u2192xB,x\u2192xS$) and transpose all matrices, which gives

The latter expression describes the redatuming of the sources from **x**_{S} at the acquisition boundary to a virtual-source position **x**_{B} inside of the medium; see Fig. 7(b).

Next, we discuss receiver redatuming. We replace $x$ by **x**_{R} in Eq. (66), which gives

This expression describes the redatuming of the receivers from **x**_{R} at the acquisition boundary to a virtual-receiver position **x**_{A} inside of the medium; see Fig. 7(c). The Green's matrix under the integral is the output of the source redatuming method, which is described by Eq. (69).

Combining Eqs. (69) and (70) gives the following expression for the combined source and receiver redatuming:

This expression redatums the actual sources and receivers from the acquisition boundary $\u2202D0$ to the virtual sources and receivers inside of the medium. Because it takes multiple scattering into account, it generalizes the classical source and receiver redatuming.^{22,57,73,74} Equation (71) resembles expressions for source-receiver interferometry^{75,76} but with the integrals along a closed boundary replaced by the integrals along the open acquisition boundary $\u2202D0$ and two of the Green's functions replaced by the propagator matrices. When the medium is known, these matrices can be numerically modelled. Alternatively, the propagator matrices can be expressed in terms of focusing functions which, at least for the acoustic situation, can be retrieved via the Marchenko method from the reflection response at the acquisition boundary. In Sec. VII, we introduce the relations between the propagator matrix and focusing functions (Sec. VII A) and derive Marchenko-type representations, which relate the focusing functions to the reflection response at the acquisition boundary (Sec. VII B).

## VII. REPRESENTATIONS WITH FOCUSING FUNCTIONS

### A. Relation between the propagator matrix and focusing functions

In Sec. II E, we indicated that there is a relation between the propagator matrix and the Marchenko-type focusing functions. Here, we derive this relation for the propagator matrix of the unified matrix-vector wave equation, assuming that the medium is lossless. Our starting point is Eq. (59), which is repeated here for convenience,

with the boundary condition

for $x$ at $\u2202D0$. From here on, we assume that the half-space above $\u2202D0$ (including $\u2202D0$) is homogeneous and isotropic, the medium below $\u2202D0$ is arbitrary inhomogeneous, anisotropic, and source-free, and $xA$ is chosen at or below $\u2202D0$ (hence, $x3,A\u2265x3,0$).

Without loss of generality, we can decompose $q(x,\omega )$ in the upper half-space into downgoing and upgoing plane waves. To this end, we defined the spatial Fourier transform $u\u0303(s,x3,\omega )$ of a space- and frequency-dependent function $u(x,\omega )$ in Eq. (20). For a function of two space variables, $u(xA,x,\omega )$, we define the spatial Fourier transform along the horizontal components of the second space variable as

Using these definitions and Parseval's theorem, we rewrite Eq. (72) as

with the boundary condition

In the upper half-space, we relate the wave-field vector $q\u0303$ to the downgoing and upgoing wave-field vectors $p\u0303+$ and $p\u0303\u2212$, respectively, via

(Refs. 55–57 and 77). Next, we renormalize the downgoing and upgoing wave-field vectors. To this end, we define the downgoing and upgoing parts of $q\u03031$ as

and rewrite Eq. (77) as

with

In Appendix A, we give explicit expressions for $D\u03031\xb1$ for the acoustic, electromagnetic, and elastodynamic situation. The substitution of Eq. (79) into Eq. (75) gives

with

Let $F\u03031(xA,s,x3,0,\omega )$ be the upper-right $N/2\xd7N/2$ matrix of the block matrix $Y\u03031(xA,s,x3,0,\omega )$. According to Eqs. (28), (82), and the definition of $D\u03031(s)$ in Eq. (79), it is defined as

This is a generalization of the definition of the acoustic focusing function for a horizontally layered medium, which is defined in Eq. (38). From Eqs. (28) and (76), it follows that $F\u03031(xA,s,x3,0,\omega )$ obeys the boundary condition

or, applying an inverse spatial and temporal Fourier transform,

for $x$ at $\u2202D0$. Hence, $F1(xA,x,t)$ is indeed a focusing function (where $xA$ is a variable and $x$ is the focal point at $\u2202D0$).

To analyse the upper-left $N/2\xd7N/2$ matrix of the block matrix $Y\u03031(xA,s,x3,0,\omega )$, we first establish some symmetry properties. The symmetry of $W\u0303$ follows from the Fourier transform of Eq. (57) for the lossless situation, hence,

For the sub-matrices of $W\u0303$, this implies that

(no summation convention) where $J11=\u2212J22=I$. The symmetry of $D\u03031\xb1$ follows from its explicit definitions in Appendix A. In all of the cases, the following symmetry holds:

When we ignore the evanescent waves at and above $\u2202D0$, then $D\u03031\xb1(s)$ is real valued; see Appendix A. Hence, given that $F\u03031(xA,s,x3,0,\omega )$ is the upper-right $N/2\xd7N/2$ matrix of $Y\u03031(xA,s,x3,0,\omega )$, we derive from Eqs. (28), (82), (87), and (88) and the structure of $D\u03031$ indicated in Eq. (79) that the upper-left $N/2\xd7N/2$ matrix is equal to $F\u03031*(xA,\u2212s,x3,0,\omega )$. Using this in Eq. (81), for the upper $N/2\xd71$ vector of $q(xA,\omega )$, we obtain

By applying Parseval's theorem again, we obtain

This relation has previously been derived via another route for the situations of acoustic waves ($q1=p$) and elastodynamic waves ($q1=v$),^{78} but without explicitly defining the focusing function $F1(xA,x,\omega )$. Here, we have an explicit expression for the Fourier transform of this focusing function [Eq. (83)]. In Sec. VII B, we use Eq. (90) as the basis for deriving the Marchenko-type Green's matrix representations.

Next, we derive a representation similar to Eq. (90) for $q2(xA,\omega )$. To this end, we define the downgoing and upgoing parts of $q\u03032$ as

with

with

Let $F\u03032(xA,s,x3,0,\omega )$ be the lower-right $N/2\xd7N/2$ matrix of the block matrix $Y\u03032(xA,s,x3,0,\omega )$. According to Eqs. (28), (95), and the definition of $D\u03032(s)$ in Eq. (92), it is defined as

$F\u03032$ is a focusing function with focusing properties similar to those of $F\u03031$, expressed in Eqs. (84) and (85).

To analyse the lower-left $N/2\xd7N/2$ matrix of the block matrix $Y\u03032(xA,s,x3,0,\omega )$, from Eqs. (88) and (93), we first derive

When we ignore the evanescent waves at and above $\u2202D0$, then $D\u03032\xb1(s)$ is real valued. Hence, given that $F\u03032(xA,s,x3,0,\omega )$ is the lower-right $N/2\xd7N/2$ matrix of $Y\u03032(xA,s,x3,0,\omega )$, we derive from Eqs. (28), (87), (95), and (97) and the structure of $D\u03032$ indicated in Eq. (92) that the lower-left $N/2\xd7N/2$ matrix is equal to $F\u03032*(xA,\u2212s,x3,0,\omega )$. Using this in Eq. (94), we obtain for the lower $N/2\xd71$ vector of $q(xA,\omega )$

Applying Parseval's theorem again, we obtain

From Eqs. (78), (80), and (91), we obtain $q\u03032\xb1=D\u03031\xb1q\u03031\xb1$. Substituting this into Eq. (98), we obtain a representation for $q2(xA,\omega )$ in terms of $q\u03031+$ and $q\u03031\u2212$. Combining this with Eq. (89) into a single equation, we obtain Eq. (81) with [using Eq. (88)]

Once $Y\u03031$ is known, the propagator matrix $W\u0303$ follows from inverting Eq. (82), according to

Equations (100) and (101), with $D\u03031(s)$ defined in Eq. (79), express the propagator matrix explicitly in terms of the focusing functions.

Using this in Eq. (101), together with

we obtain an explicit expression for $W\u0303(xA,s,x3,0,\omega )$ in terms of the acoustic focusing functions. Transforming this to the space-frequency domain, replacing $s3,0$ by $(1/\omega )H1$ (where $H1$ is the square-root of the Helmholtz operator $\omega 2/c02+\u2202\alpha \u2202\alpha $ in the homogeneous upper half-space^{54,56,57}) yields^{60}

for $x$ at $\u2202D0$, where $\u2111$ denotes the imaginary part.

### B. Marchenko-type Green's matrix representations with focusing functions

We use Eqs. (90) and (99) as the starting point for deriving two Marchenko-type Green's matrix representations. We follow a similar procedure as in Ref. 78 (Appendix B), generalized for the different wave phenomena considered in this paper. We replace the $N/2\xd71$ vector $q1(x,\omega )$ by a modified version $\Gamma 12(x,xS,\omega )$ of the $N/2\xd7N/2$ Green's matrix $G12(x,xS,\omega )$. Here, $G12(x,xS,\omega )$ stands for the $q1$-type field observed at $x$ in response to a unit $d2$-type source at **x**_{S}. We choose $xS=(xH,S,x3,S)$ in the upper half-space at a vanishing distance *ϵ* above $\u2202D0$, hence, $x3,S=x3,0\u2212\u03f5$. Our aim is to modify this Green's matrix such that for $x$ at $\u2202D0$ (i.e., just below the source), its downgoing part simplifies to

First, we derive the properties of the downgoing part of $G12(x,xS,\omega )$ for $x$ at $\u2202D0$. To this end, consider the inverse of Eq. (79), which reads

for $\u2009x3\u2264x3,0$ with

The upper-right $N/2\xd7N/2$ matrix $(\Delta \u03031)\u22121$ transforms $q\u03032$ into a downgoing field vector $q\u03031+$. In a similar way, this matrix transforms a unit $d2$-type source in a homogeneous half-space into the downgoing part of the Green's matrix $G\u030312$ just below this source, according to

The explicit expressions for $(\Delta \u03031)\u22121$ are given in Appendix A. In Eq. (110), the source is located at $(0,x3,S)$. Next, we consider $G12(x,xS,\omega )$ for a laterally shifted source position $(xH,S,x3,S)$. Applying a spatial Fourier transform along the horizontal source coordinate $xH,S$, using Eq. (74) with $xH$ replaced by $xH,S$, yields $G\u030312(x,s,x3,S,\omega )$. For the downgoing part just below the source, we obtain a phase-shifted version of the Green's matrix of Eq. (110), according to

This suggests that the modified Green's matrix be defined as

such that

Transforming this back to the space domain yields, indeed, Eq. (107). Next, we define the reflection response $R12(x,xS,\omega )$ of the medium below $\u2202D0$ as the upgoing part of the modified Green's matrix $\Gamma 12(x,xS,\omega )$ for $x$ at $\u2202D0$, hence,

Substituting $q1(xA,\omega )=\Gamma 12(xA,xS,\omega )$ and $q1\xb1(x,\omega )=\Gamma 12\xb1(x,xS,\omega )$ into Eq. (90), using Eqs. (107) and (114), gives

for $x3,A\u2265x3,0$. This is the first Marchenko-type representation. In a similar way, we derive a second Marchenko-type representation for a modified version $\Gamma 22(x,xS,\omega )$ of the Green's matrix $G22(x,xS,\omega )$. We derive the properties of the downgoing part of $G22(x,xS,\omega )$ for $x$ at $\u2202D0$. Consider the inverse of Eq. (92), which reads

for $x3\u2264x3,0$ with

The upper-right $N/2\xd7N/2$ matrix $\u2212(\Delta \u03032)\u22121D\u03032\u2212$ transforms $q\u03032$ into a downgoing field vector $q\u03032+$. In a similar way, this matrix transforms a unit $d2$-type source in a homogeneous half-space into the downgoing part of the Green's matrix $G\u030322$ just below this source, according to

The explicit expressions for $\u2212(\Delta \u03032)\u22121D\u03032\u2212$ are given in Appendix A. Steps similar to those below Eq. (110) lead to

This suggests that the modified Green's matrix be defined as

such that

We define the reflection response $R22(x,xS,\omega )$ as

for $x$ at $\u2202D0$. Substituting $q2(xA,\omega )=\Gamma 22(xA,xS,\omega )$ and $q2\xb1(x,\omega )=\Gamma 22\xb1(x,xS,\omega )$ into Eq. (99), using Eqs. (121) and (122), gives

for $x3,A\u2265x3,0$. This is the second Marchenko-type representation.

or in the space-frequency domain [using Eq. (18)],

For this situation, $R12(x,xS,\omega )$ and $R22(x,xS,\omega )$ in the representations of Eqs. (115) and (123) are the upgoing parts of $2Gp,f(x,xS,\omega )$ and $2Gv,q(x,xS,\omega )$, respectively, for $x$ at $\u2202D0$. Note that according to Eq. (43), $Gv,q(x,xS,\omega )=\u2212Gp,f(xS,x,\omega )$.

In their general form, the representations of Eqs. (115) and (123) are generalizations of the previously derived representations for the 3D Marchenko method for acoustic^{72,79,80} and elastodynamic wave fields^{81,82} (but note that the subscripts 1 and 2 of the focusing functions have a different meaning than in those papers; here, they refer to the wave-field components $q1$ and $q2$). In most of the previous work on the 3D Marchenko method, one of the underlying assumptions is that the wave field can be decomposed into downgoing and upgoing waves in the interior of the medium. Only recently, several authors proposed to avoid the decomposition inside the medium.^{78,83,84} The representations discussed in this section expand on this. In these representations, decomposition into downgoing and upgoing waves and negligence of the evanescent waves occurs only in the upper half-space. However, inside the medium, no wave-field decomposition takes place. Moreover, the evanescent waves inside the medium (for example, in high-velocity layers) are accounted for by the representations of Eqs. (115) and (123). These representations, transformed to the time domain, form the basis for the development of Marchenko schemes, aiming at resolving the focusing functions $F1$ and $F2$ from the reflection responses $R12$ and $R22$, respectively, at the acquisition boundary $\u2202D0$. Such schemes have been successfully developed for precritical acoustic data.^{45–47,85–88} For more complex situations, research on retrieving the focusing functions from the reflection responses is ongoing (for example, Ref. 89). Once the focusing functions are found, they can be used to define the propagator matrix via Eqs. (100) and (101), which can, subsequently, be used in the representations of Secs. V and VI.

## VIII. CONCLUSIONS

We have discussed different types of wave-field representation in a systematic way. Classical wave-field representations contain Green's functions. By starting with a unified matrix-vector wave equation, we have formulated wave-field representations with Green's matrices, analogous to the classical representations. For example, the classical Kirchhoff-Helmholtz integral follows as a special case of the unified representation with the Green's matrix. Another special case is the classical homogeneous Green's function representation.

Using the same matrix-vector formalism, we formulated wave-field representations with propagator matrices. Unlike a Green's matrix, a propagator matrix depends only on the medium parameters between the two depth levels for which this matrix is defined. The representations with the propagator matrices have a similar form as those with the Green's matrices. An important difference is that the boundary integrals in the representations with the propagator matrices are single sided.

We also formulated representations containing a mix of Green's matrices and propagator matrices. A special case is the single-sided homogeneous Green's function representation as the counterpart of the classical closed-boundary homogeneous Green's function representation.

We have shown that the propagator matrix is related to Marchenko-type focusing functions. We have used this relation to reformulate the representations with the propagator matrix into representations with focusing functions. For the acoustic situation, these focusing functions can be retrieved from the single-sided reflection response of the medium by applying the Marchenko method. For more complex situations, research on retrieving these focusing functions from the reflection response is ongoing. Once the focusing functions are known, they can be used to construct the propagator matrix. Subsequently, the propagator matrix can be used in the representations to obtain the wave field inside of the medium which, in turn, can be used, for example, for imaging or monitoring. Unlike earlier imaging methods, which use the propagator matrix, this is a data-driven approach (because the propagator matrix is retrieved from the reflection response) and, hence, it does not require a detailed model of the medium.

## ACKNOWLEDGMENTS

The author thanks Roel Snieder, Sjoerd de Ridder, and Evert Slob for fruitful discussions and two anonymous reviewers for their positive evaluations of the original manuscript and useful suggestions for further improvements. This research is funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant No. 742703).

### APPENDIX A: EXPLICIT EXPRESSIONS FOR SOME MATRICES IN THE HOMOGENEOUS ISOTROPIC UPPER HALF-SPACE

We derive the explicit expressions for $D\u03031\xb1,\u2009\u2009(\Delta \u03031)\u22121$, and $\u2212(\Delta \u03032)\u22121D\u03032\u2212$ in the homogeneous isotropic upper half-space $x3\u2264x3,0$ for the acoustic, electromagnetic, and elastodynamic situations.

##### 1. Unified expressions in the horizontal slowness domain

For the lossless homogeneous isotropic upper half-space $x3\u2264x3,0$, we transform the unified matrix-vector wave equation (1) to the horizontal slowness domain using the spatial Fourier transform of Eq. (20). This yields

where

The symmetry property, defined by Eq. (3), transforms to

We decompose matrix $A\u0303$ as^{55,57,77}

where

With a proper scaling of the columns of matrix $L\u0303$, these matrices obey the following symmetry properties:

The substitution of Eq. (A4) into Eq. (A1) and premultiplying all terms with $L\u0303\u22121$ yields

where

where $p\u0303+$ and $p\u0303\u2212$ are wave-field vectors containing the downgoing (+) and upgoing (-) waves, respectively. Following Eqs. (80), (93), (109), and (117), we define

##### 2. Acoustic wave equation

For the acoustic wave equation, the 1 × 1 sub-matrices of $A\u0303$ are

where the radial slowness *s _{r}* is defined as

Here, *κ* and *ρ* are the compressibility and mass density of the homogeneous upper half-space, respectively. For convenience, here and in the remainder of this appendix, we do not use the subscripts 0 to indicate the parameters of the upper half-space. The 1 × 1 sub-matrices of $L\u0303$ and $\Lambda \u0303$ are

where the vertical slowness *s*_{3} is defined as

where the propagation velocity *c* is defined as $c=1/\kappa \rho $. The two expressions in Eq. (A15) represent the propagating and evanescent waves. On substitution of Eq. (A14) into Eqs. (A10) and (A11), we obtain

##### 3. Electromagnetic wave equation

For the electromagnetic wave equation, the 2 × 2 sub-matrices of $A\u0303$ are

Here, *ε* and *μ* are the permittivity and permeability of the upper half-space, respectively. The 2 × 2 sub-matrices of $L\u0303$ and $\Lambda \u0303$ are

where *s _{r}* and

*s*

_{3}are defined in Eqs. (A13) and (A15), respectively, but this time with the propagation velocity

*c*defined as $c=1/\epsilon \mu $. On substitution of Eqs. (A20) and (A21) into Eqs. (A10) and (A11), we obtain

##### 4. Elastodynamic wave equation

For the elastodynamic wave equation, the 3 × 3 sub-matrices of $A\u0303$ are

where

where *λ* and *μ* are the Lamé parameters and *ρ* is the mass density of the upper half-space. The 3 × 3 sub-matrices of $L\u0303$ and $\Lambda \u0303$ are

where *s _{r}* is defined as in Eq. (A13), and the vertical slownesses $s3P$ and $s3S$ are defined as

Here, *c _{P}* and

*c*are the

_{S}*P*- and

*S*-wave velocities of the upper half-space, respectively, defined as $cP=(\lambda +2\mu )/\rho $ and $cS=\mu /\rho $. On substitution of Eqs. (A32) and (A33) into Eqs. (A10) and (A11), we obtain

where $S=s3Ps3S+sr2$.

## References

**44**(4), 356–374 (1852); and

**47**(3), 161–221 (1854).]