Propagator and transfer matrices, Marchenko focusing functions and their mutual relations

TL;DR

This paper explores the relationships between propagator matrices, transfer matrices, and Marchenko focusing functions in seismic wave field extrapolation, aiming to improve imaging and wave field retrieval in heterogeneous media.

Contribution

It presents a unified framework expressing these operators and functions, avoiding common approximations, and discusses their mutual relations for enhanced seismic imaging techniques.

Findings

Expressed Marchenko focusing functions using propagator matrices.

Established the mutual relations between propagator, transfer matrices, and focusing functions.

Discussed implications for wave field focusing and Green's function retrieval.

Abstract

Many seismic imaging methods use wave field extrapolation operators to redatum sources and receivers from the surface into the subsurface. We discuss wave field extrapolation operators that account for internal multiple reflections, in particular propagator matrices, transfer matrices and Marchenko focusing functions. A propagator matrix is a square matrix that `propagates' a wave-field vector from one depth level to another. It accounts for primaries and multiples and holds for propagating and evanescent waves. A Marchenko focusing function is a wave field that focuses at a designated point in space at zero time. Marchenko focusing functions are useful for retrieving the wave field inside a heterogeneous medium from the reflection response at its surface. By expressing these focusing functions in terms of the propagator matrix, the usual approximations (such as ignoring evanescent waves) are avoided. While a propagator matrix acts on the full wave-field vector, a transfer matrix (according to the definition employed in this paper)`transfers' a decomposed wave-field vector (containing downgoing and upgoing waves) from one depth level to another. It can be expressed in terms of decomposed Marchenko focusing functions. We present propagator matrices, transfer matrices and Marchenko focusing functions in a consistent way and discuss their mutual relations. In the main text we consider the acoustic situation and in the appendices we discuss other wave phenomena. Understanding these mutual connections may lead to new developments of Marchenko theory and its applications in wave field focusing, Green's function retrieval and imaging.

Equations

$$\displaystyle\partial_{3}{\bf q}={{\mbox{\boldmath${\cal A}$}}}\,{\bf q}+{\bf d}$$

$$\mathbf{q}(\mathbf{x},\omega )=\left(\begin{array}{c}p\\ v_3\end{array}\right)(\mathbf{x},\omega ).$$

$$\displaystyle{{\mbox{\boldmath${\cal A}$}}}({\bf x},\omega)$$

$${\mathcal{H}}_2(\mathbf{x},\omega )=k^2(\mathbf{x},\omega )+{\partial }_{\alpha }{\partial }_{\alpha },$$

$$k^2(\mathbf{x},\omega )={\omega }^2\kappa \rho -\frac{3({\partial }_{\alpha }\rho )({\partial }_{\alpha }\rho )}{4{\rho }^2}+\frac{({\partial }_{\alpha }{\partial }_{\alpha }\rho )}{2\rho }$$

$$\mathbf{d}(\mathbf{x},\omega )=\left(\begin{array}{c}{\hat{f}}_3\\ \frac{1}{i\omega }{\partial }_{\alpha }(\frac{1}{\rho }{\hat{f}}_{\alpha })+\hat{q}\end{array}\right)(\mathbf{x},\omega ).$$

$$\displaystyle\partial_{3}{\bf W}({\bf x},{\bf x}_{F})={{\mbox{\boldmath${\cal A}$}}}({\bf x}){\bf W}({\bf x},{\bf x}_{F}),$$

$$\mathbf{W}(\mathbf{x},{\mathbf{x}}_F){|}_{x_3=x_{3,F}}=\mathbf{I}\delta ({\mathbf{x}}_{\mathrm{H}}-{\mathbf{x}}_{\mathrm{H},F}),$$

$$\mathbf{q}(\mathbf{x})$$

$$\mathbf{W}(\mathbf{x},{\mathbf{x}}_F)=\left(\begin{array}{cc}{W}^{p,p}& W^{p,v}\\ W^{v,p}& W^{v,v}\end{array}\right)(\mathbf{x},{\mathbf{x}}_F),$$

$$\mathbf{W}({\mathbf{x}}^{\prime },{\mathbf{x}}_F)$$

$$\mathbf{I}\delta ({\mathbf{x}}_{\mathrm{H}}^{\prime }-{\mathbf{x}}_{\mathrm{H},F})$$

$$\tilde{u}(\mathbf{s},x_3,\omega )={\int }_{{\mathbb{R}}^2}\exp \{-i\omega \mathbf{s}\cdot {\mathbf{x}}_{\mathrm{H}}\}u({\mathbf{x}}_{\mathrm{H}},x_3,\omega ){\mathrm{d}}^2{\mathbf{x}}_{\mathrm{H}},$$

$$\tilde{\mathbf{W}}(\mathbf{s},x_3^{\prime },x_{3,F})=\tilde{\mathbf{W}}(\mathbf{s},x_3^{\prime },x_3)\tilde{\mathbf{W}}(\mathbf{s},x_3,x_{3,F}).$$

$$u(\mathbf{s},x_3,\tau )=\frac{1}{\pi }\Re {\int }_0^{\infty }\tilde{u}(\mathbf{s},x_3,\omega )\exp \{-i\omega \tau \}\mathrm{d}\omega ,$$

$$\mathbf{W}(\mathbf{s},x_3^{\prime },x_{3,F},\tau )=\mathbf{W}(\mathbf{s},x_3^{\prime },x_3,\tau )\ast \mathbf{W}(\mathbf{s},x_3,x_{3,F},\tau ),$$

$$\displaystyle{{{\mbox{\boldmath${\cal A}$}}}}={{{\mbox{\boldmath${\cal L}$}}}}{\bf\Lambda}{{{\mbox{\boldmath${\cal L}$}}}}^{-1},$$

$$\mathbf{\Lambda }=\left(\begin{array}{cc}i{\mathcal{H}}_1& 0\\ 0& -i{\mathcal{H}}_1\end{array}\right),{\mathcal{H}}_1={\rho }^{1/2}{\mathcal{H}}_2^{1/2}{\rho }^{-1/2},$$

$$\displaystyle{{{\mbox{\boldmath${\cal L}$}}}}=\begin{pmatrix}1&1\\ \frac{1}{\omega\rho}{\cal H}_{1}&-\frac{1}{\omega\rho}{\cal H}_{1}\end{pmatrix},\quad{{{\mbox{\boldmath${\cal L}$}}}}^{-1}=\frac{1}{2}\begin{pmatrix}1&\omega{\cal H}_{1}^{-1}\rho\\ 1&-\omega{\cal H}_{1}^{-1}\rho\end{pmatrix}$$

$${\int }_{{\mathbb{R}}^2}\{{\mathcal{H}}_2^{1/2}g({\mathbf{x}}_{\mathrm{H}})\}h({\mathbf{x}}_{\mathrm{H}}){\mathrm{d}}^2{\mathbf{x}}_{\mathrm{H}}={\int }_{{\mathbb{R}}^2}g({\mathbf{x}}_{\mathrm{H}})\{{\mathcal{H}}_2^{1/2}h({\mathbf{x}}_{\mathrm{H}})\}{\mathrm{d}}^2{\mathbf{x}}_{\mathrm{H}}$$

$$\displaystyle{\bf q}({\bf x})={{{\mbox{\boldmath${\cal L}$}}}}({\bf x}){\bf p}({\bf x}),$$

$$\mathbf{p}(\mathbf{x})=\left(\begin{array}{c}{p}^{+}\\ p^{-}\end{array}\right)(\mathbf{x}).$$

$$p(\mathbf{x})=p^{+}(\mathbf{x})+p^{-}(\mathbf{x}),$$

$$\mathbf{q}(\mathbf{x})={\int }_{\partial {\mathbb{D}}_F}\mathbf{Y}(\mathbf{x},{\mathbf{x}}_F)\mathbf{p}({\mathbf{x}}_F){\mathrm{d}}^2{\mathbf{x}}_F,$$

$$\displaystyle{\bf Y}({\bf x},{\bf x}_{F})={\bf W}({\bf x},{\bf x}_{F}){{{{\mbox{\boldmath${\cal L}$}}}}}({\bf x}_{F}),$$

$$\mathbf{Y}(\mathbf{x},{\mathbf{x}}_F)=\left(\begin{array}{cc}{W}^{p,p}& W^{p,v}\\ W^{v,p}& W^{v,v}\end{array}\right)(\mathbf{x},{\mathbf{x}}_F)\left(\begin{array}{cc}1& 1\\ \frac{1}{\omega {\rho }_0}{\mathcal{H}}_1& -\frac{1}{\omega {\rho }_0}{\mathcal{H}}_1\end{array}\right)({\mathbf{x}}_F).$$

$$\mathbf{Y}(\mathbf{x},{\mathbf{x}}_F)=\left(\begin{array}{cc}{Y}^{p,+}& Y^{p,-}\\ Y^{v,+}& Y^{v,-}\end{array}\right)(\mathbf{x},{\mathbf{x}}_F),$$

$$Y^{p,\pm }(\mathbf{x},{\mathbf{x}}_F)=W^{p,p}(\mathbf{x},{\mathbf{x}}_F)\pm \frac{1}{\omega {\rho }_0}{\mathcal{H}}_1({\mathbf{x}}_F)W^{p,v}(\mathbf{x},{\mathbf{x}}_F),$$

$$Y^{v,\pm }(\mathbf{x},{\mathbf{x}}_F)=W^{v,p}(\mathbf{x},{\mathbf{x}}_F)\pm \frac{1}{\omega {\rho }_0}{\mathcal{H}}_1({\mathbf{x}}_F)W^{v,v}(\mathbf{x},{\mathbf{x}}_F).$$

$$Y^{p,-}(\mathbf{x},{\mathbf{x}}_F){|}_{x_3=x_{3,F}}=\delta ({\mathbf{x}}_{\mathrm{H}}-{\mathbf{x}}_{\mathrm{H},F}),$$

Full text

††volume: 999

{summary}

Many seismic imaging methods use wave field extrapolation operators to redatum sources and receivers from the surface into the subsurface. We discuss wave field extrapolation operators that account for internal multiple reflections, in particular propagator matrices, transfer matrices and Marchenko focusing functions.

A propagator matrix is a square matrix that ‘propagates’ a wave-field vector from one depth level to another. It accounts for primaries and multiples and holds for propagating and evanescent waves.

A Marchenko focusing function is a wave field that focuses at a designated point in space at zero time. Marchenko focusing functions are useful for retrieving the wave field inside a heterogeneous medium from the reflection response at its surface. By expressing these focusing functions in terms of the propagator matrix, the usual approximations (such as ignoring evanescent waves) are avoided.

While a propagator matrix acts on the full wave-field vector, a transfer matrix (according to the definition employed in this paper) ‘transfers’ a decomposed wave-field vector (containing downgoing and upgoing waves) from one depth level to another. It can be expressed in terms of decomposed Marchenko focusing functions.

We present propagator matrices, transfer matrices and Marchenko focusing functions in a consistent way and discuss their mutual relations. In the main text we consider the acoustic situation and in the appendices we discuss other wave phenomena. Understanding these mutual connections may lead to new developments of Marchenko theory and its applications in wave field focusing, Green’s function retrieval and imaging.

1 Introduction

In many seismic imaging methods, wave field extrapolation is used to redatum sources and receivers from the surface to a depth level in the subsurface. In most cases the redatuming process is based on one-way wave field extrapolation operators, which only account for primaries. To account for internal multiple reflections in redatuming, more advanced wave field extrapolation operators are required. This paper is not about the redatuming process itself, but about wave field extrapolation operators that account for internal multiples. In particular, we discuss propagator and transfer matrices, Marchenko focusing functions and their mutual relations.

In elastodynamic wave theory, a propagator matrix is a square matrix that ‘propagates’ a wave-field vector from one depth level to another. It was originally introduced in geophysics for horizontally layered media (Thomson, 1950; Haskell, 1953; Gilbert & Backus, 1966) and later extended for laterally varying media (Kennett, 1972). It has been used for modelling surface waves (Woodhouse, 1974) and reflection and transmission responses of heterogeneous media (Haines, 1988; Kennett et al., 1990; Koketsu et al., 1991; Takenaka et al., 1993). It has also been proposed as an operator for accurate seismic imaging schemes, accounting for high propagation angles (Kosloff & Baysal, 1983) and internal multiple reflections (Wapenaar & Berkhout, 1986). The wave-field vector that the propagator matrix acts on contains components of the full wave field (e.g. particle velocity and stress). Here ‘full’ means that the wave field implicitly consists of downgoing and upgoing, propagating and evanescent waves.

A Marchenko focusing function is a wave field that focuses at a designated point in space at zero time, accounting for primaries and multiples. Marchenko focusing functions were originally introduced to retrieve the wave field inside a horizontally layered medium from the reflection response at the boundary of that medium (Rose, 2001, 2002; Broggini & Snieder, 2012; Slob et al., 2014). This has been extended for laterally varying media (Wapenaar et al., 2013), under the assumption that the wave field inside the medium can be decomposed into downgoing and upgoing components and that evanescent waves can be neglected. It has recently been shown that the propagator matrix can be expressed in terms of Marchenko focusing functions and vice versa (Wapenaar & de Ridder, 2022). Via this relation, the usual assumptions underlying the focusing functions (such as ignoring evanescent waves) are circumvented.

In this paper we define a transfer matrix as a square matrix that ‘transfers’ decomposed wave-field vectors (explicitly containing downgoing and upgoing waves) from one depth level to another (Born & Wolf, 1965; Katsidis & Siapkas, 2002). It is different from the propagator matrix, which acts on full wave-field vectors (but please note that in the literature there is not a clear distinction between the use of the terminologies ‘propagator matrix’ and ‘transfer matrix’). It has recently been shown that the transfer matrix can be expressed in terms of decomposed Marchenko focusing functions (Dukalski et al., 2022a, b), an insight that is expected to be useful in further analysis of the minimum-phase property of elastodynamic focusing functions (Reinicke et al., 2023) and beyond.

The aim of this paper is to present propagator matrices, transfer matrices and Marchenko focusing functions in a consistent way and to discuss their mutual relations. We aim to set up the theory as general as possible, accounting for lateral and vertical variations of the medium parameters, accounting for evanescent waves and taking dissipation into account. Whereas in the main text we consider acoustic waves, in the appendices we generalize the theory for other wave phenomena. The numerical examples, which are meant as illustrations of the different quantities and their relations, are restricted to oblique acoustic plane waves in a lossless horizontally layered medium.

We hope that this consistent treatment will contribute to the understanding of the mutual connections and provide insight in the assumptions and approximations that underlie Marchenko-type wave field retrieval schemes and how to cope with them (Slob, 2016; Dukalski et al., 2019, 2022b; Reinicke et al., 2020, 2023; Elison et al., 2020; Diekmann & Vasconcelos, 2021; Wapenaar et al., 2021; Kiraz et al., 2023). Moreover, we hope to stimulate new research directions.

The setup of this paper is as follows. In section 2 we discuss the $2\times 2$ propagator matrix for acoustic wave fields and its relation with acoustic Marchenko focusing functions. The advantage of concentrating on the acoustic situation is that all expressions are relatively simple and yet contain all essential aspects. In section 3 we discuss the $2\times 2$ transfer matrix for acoustic wave fields and its relation with decomposed acoustic Marchenko focusing functions. In section 4 we present some conclusions.

Appendices A and B are generalisations of sections 2 and 3 for other wave phenomena. Here, the propagator and transfer matrices are $N\times N$ matrices, with $N$ ranging from 2 for acoustic waves to 12 for seismoelectric waves; the Marchenko focusing functions are $\frac{N}{2}\times\frac{N}{2}$ matrices. We derive their mutual relations by exploiting general symmetry properties, which are derived in Appendix C. The appendices not only cover classical waves, but also quantum mechanical waves obeying the Schrödinger equation ($N=2$) and the Dirac equation ($N=4$).

2 Acoustic propagator matrix and focusing functions

2.1 Acoustic matrix-vector wave equation

Our starting point is the following matrix-vector wave equation in the space-frequency domain

[TABLE]

(Woodhouse, 1974; Corones, 1975; Ursin, 1983; Kosloff & Baysal, 1983; Fishman & McCoy, 1984; Wapenaar & Berkhout, 1986; de Hoop, 1996). In Appendix A we discuss this equation for a range of wave phenomena. Here we consider acoustic waves. For this situation, ${\bf q}$ is a vector containing the wave field components $p$ (acoustic pressure) and $v_{3}$ (vertical component of the particle velocity), both as a function of the space coordinate vector ${\bf x}=(x_{1},x_{2},x_{3})$ (with positive $x_{3}$ denoting depth) and the angular frequency $\omega$, hence,

[TABLE]

Operator $\partial_{3}$ stands for the partial differential operator $\partial/\partial x_{3}$. The space- and frequency-dependent operator matrix ${\cal A}$ is defined as

[TABLE]

where $\kappa({\bf x},\omega)$ is the compressibility, $\rho({\bf x},\omega)$ the mass density and $i$ the imaginary unit. Operator $\partial_{\alpha}$ stands for the partial differential operator $\partial/\partial x_{\alpha}$. Greek subscripts take on the values 1 and 2 and Einstein’s summation convention applies to repeated subscripts, unless otherwise noted. In general the medium may be dissipative, meaning that $\kappa$ and $\rho$ may be frequency-dependent and complex-valued, with (for positive $\omega$) $\Im(\kappa)\geq 0$ and $\Im(\rho)\geq 0$, where $\Im$ denotes the imaginary part. For later convenience we rewrite the operator matrix as follows

[TABLE]

Here ${\cal H}_{2}({\bf x},\omega)$ is the Helmholtz operator, defined as

[TABLE]

with wavenumber $k({\bf x},\omega)$ defined via

[TABLE]

(Brekhovskikh, 1960; Wapenaar et al., 2001). Finally, vector ${\bf d}$ in equation (1) contains source terms, according to

[TABLE]

Here $\hat{f}_{\alpha}({\bf x},\omega)$ and $\hat{f}_{3}({\bf x},\omega)$ are the horizontal and vertical components, respectively, of the external force density (the hats are used to distinguish external force components from focusing functions), and $\hat{q}({\bf x},\omega)$ is the volume injection-rate density (where $\hat{q}$ is to be distinguished from the wave field vector ${\bf q}$). From here onward we simplify the notation by not explicitly mentioning the frequency-dependency in the argument lists.

2.2 Acoustic propagator matrix

We define a boundary $\partial\mathbb{D}_{F}$ at depth level $x_{3}=x_{3,F}$. We define a coordinate vector ${\bf x}_{F}$ at this boundary as ${\bf x}_{F}=(x_{1,F},x_{2,F},x_{3,F})$ (with fixed $x_{3,F}$). We introduce the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ as a solution of wave equation (1) for the source-free situation, according to

[TABLE]

with boundary condition

[TABLE]

where ${\bf I}$ is the identity matrix and ${\bf x}_{\rm H}$ and ${\bf x}_{{\rm H},F}$ denote the horizontal coordinates of ${\bf x}$ and ${\bf x}_{F}$, respectively, hence ${\bf x}_{\rm H}=(x_{1},x_{2})$ and ${\bf x}_{{\rm H},F}=(x_{1,F},x_{2,F})$. Since equations (1) and (8) are both linear, Huygens’ superposition principle can be applied to get a representation for ${\bf q}({\bf x})$ in terms of ${\bf W}({\bf x},{\bf x}_{F})$. For a given depth level $x_{3}$, assuming there are no sources for ${\bf q}({\bf x})$ between $x_{3,F}$ and $x_{3}$, we obtain

[TABLE]

(Gilbert & Backus, 1966; Kennett, 1972; Woodhouse, 1974). Note that equation (10) expresses the ‘propagation’ of ${\bf q}$ from depth level $x_{3,F}$ to depth level $x_{3}$, which is why ${\bf W}({\bf x},{\bf x}_{F})$ is called the propagator matrix. It is partitioned as follows

[TABLE]

where $W^{p,p}$, $W^{p,v}$, $W^{v,p}$ and $W^{v,v}$ are the scalar components of the propagator matrix. For each of these components, the second superscript refers to the wave field component ($p$ or $v_{3}$) it acts on at ${\bf x}_{F}$, whereas the first superscript refers to the wave field component it contributes to at ${\bf x}$. Equation (10) is illustrated in the upper-left frame of Figure 1. The solid line at $x_{3,F}$ denotes the boundary $\partial\mathbb{D}_{F}$ (not necessarily a physical boundary). The medium below $\partial\mathbb{D}_{F}$ may be inhomogeneous and dissipative. The dashed line at $x_{3}$ indicates an arbitrary depth level inside the inhomogeneous medium.

By applying equation (10) recursively, it follows that ${\bf W}$ obeys the following recursive expression

[TABLE]

where ${{{\partial\mathbb{D}}}}$ is a horizontal boundary at a constant depth level $x_{3}$. By taking $x_{3}^{\prime}=x_{3,F}$, we obtain from equations (9) and (12)

[TABLE]

from which it follows that ${\bf W}({\bf x}_{F},{\bf x})$ is the inverse of ${\bf W}({\bf x},{\bf x}_{F})$.

The propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ accounts for primaries and multiples between $x_{3,F}$ and $x_{3}$ and it holds for propagating and evanescent waves (for example, Woodhouse (1974) uses the elastodynamic version of the propagator matrix to analyse surface waves). Evanescent field components may lead to instability and should be handled with care (Kennett & Kerry, 1979). Since the underlying wave equation is based on the explicit Helmholtz operator ${\cal H}_{2}$ (rather than on its square-root, appearing in one-way wave equations), Kosloff & Baysal (1983) argue that the numerical evaluation of equation (10) converges much faster and for higher propagation angles than schemes based on one-way wave equations. They exploit this property in wide-angle imaging of seismic reflection responses. They use filters to eliminate evanescent and downgoing waves, so they do not exploit the fact that the propagator matrix can handle multiply reflected and evanescent waves. Wapenaar & Berkhout (1986) propose a seismic imaging scheme based on the propagator matrix that handles internal multiple reflections. Since their scheme is very sensitive to the chosen background model it has not found broad applications. In section 2.3 we show that the propagator matrix can be expressed in terms of Marchenko focusing functions. For a lossless medium, these focusing functions can be derived from seismic reflection data and a smooth background model (section 2.4). Hence, this leads to a propagator matrix that can be used for seismic imaging, which properly handles internal multiple reflections without being highly sensitive to the background model.

We conclude this section with a numerical illustration of the propagator matrix for the horizontally layered lossless medium of Figure 2(a). In each layer the propagation velocity $c=1/\sqrt{\kappa\rho}$ is shown (in m/s). We define the spatial Fourier transformation of a function $u({\bf x},\omega)$ along the horizontal coordinate ${\bf x}_{\rm H}$ for constant $x_{3}$ as

[TABLE]

where ${{\bf s}}=(s_{1},s_{2})$ is the horizontal slowness vector and ${\mathbb{R}}$ is the set of real numbers. For a horizontally layered medium, this transformation decomposes $u({\bf x},\omega)$ into independent plane waves, with propagation angle $\theta$ (with respect to the vertical axis) obeying $\sin\theta=c|{\bf s}|$. We apply this transformation to the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$, choosing ${\bf x}_{F}=(0,0,x_{3,F})$. This yields the transformed propagator matrix $\tilde{\bf W}({\bf s},x_{3},x_{3,F})$, with boundary condition $\tilde{\bf W}({\bf s},x_{3,F},x_{3,F})={\bf I}$. Analogous to equation (12) it obeys the recursive expression

[TABLE]

Next, we define the inverse temporal Fourier transformation for constant ${\bf s}$ and $x_{3}$ as

[TABLE]

where $\Re$ denotes the real part and $\tau$ is the intercept time (Stoffa, 1989). Applying this transformation to $\tilde{\bf W}({\bf s},x_{3},x_{3,F})$ we obtain ${\bf W}({\bf s},x_{3},x_{3,F},\tau)$, with boundary condition ${\bf W}({\bf s},x_{3,F},x_{3,F},\tau)={\bf I}\delta(\tau)$ and with ${\bf W}({\bf s},x_{3},x_{3,F},\tau)$ obeying the recursive expression

[TABLE]

where the inline asterisk denotes temporal convolution. Although the numerical modeling is most efficiently done in the slowness-frequency domain (using equation (15)), the results are more conveniently interpreted when displayed in the slowness intercept-time domain. Setting $s_{2}=0$, the components $W^{p,p}(s_{1},x_{3},x_{3,F},\tau)$ and $W^{p,v}(s_{1},x_{3},x_{3,F},\tau)$, with boundary conditions $W^{p,p}(s_{1},x_{3,F},x_{3,F},\tau)=\delta(\tau)$ and $W^{p,v}(s_{1},x_{3,F},x_{3,F},\tau)=0$, are shown in Figures 2(b) and 2(c) for fixed $s_{1}=1/3500$ s/m, as a function of intercept time $\tau$ and depth $x_{3}$. To get a smooth display, at each depth the components are convolved with a Ricker wavelet with a central frequency of 50 Hz. The upper traces at $x_{3}=x_{3,F}=0$ m represent the aforementioned boundary conditions. Note that $W^{p,p}$ and $W^{p,v}$ are, for each depth $x_{3}$, even and odd functions, respectively, of intercept time $\tau$. The recursive character, described by equation (17), is manifest in Figures 2(b) and 2(c). The propagation velocity in the layer between $x_{3}=760$ m and $x_{3}=800$ m equals 3600 m/s, which implies that for the chosen horizontal slowness $s_{1}=1/3500$ s/m we have $\sin\theta=cs_{1}>1$ (i.e., $\theta$ is complex-valued). Hence, waves become ‘evanescent’ in this layer. The wave field tunnels through this layer and the amplitudes below this layer are higher than above it. In general, evanescent field components of the propagator matrix should be handled with care, because next to exponentially decaying terms they contain exponentially growing terms that may cause numerical inaccuracies (Kennett & Kerry, 1979). In practice this means that beyond a certain horizontal slowness the wave field should be tapered to zero.

2.3 Relation between acoustic propagator matrix and Marchenko focusing functions

From here onward we assume that the medium at and above $\partial\mathbb{D}_{F}$ is homogeneous and may be dissipative, with mass density $\rho_{0}$ and propagation velocity $c_{0}$. The medium below $\partial\mathbb{D}_{F}$ may be inhomogeneous and dissipative, and it is source-free.

In preparation for defining the focusing functions, we decompose operator matrix ${\cal A}$ (in the space-frequency domain) as follows

[TABLE]

with

[TABLE]

(Corones, 1975; Fishman & McCoy, 1984; Wapenaar & Berkhout, 1986; de Hoop, 1996). The square-root operator ${\cal H}_{2}^{1/2}$ is symmetric in the following sense

[TABLE]

(Wapenaar et al., 2001), where $g({\bf x}_{\rm H})$ and $h({\bf x}_{\rm H})$ are test functions in the horizontal plane with ‘sufficient decay at infinity’. Operator ${\cal H}_{1}$, as defined in equation (19) is not symmetric, but operator $\frac{1}{\rho}{\cal H}_{1}$ and its inverse, both appearing in equation (20), are symmetric. We use the operator matrix ${\cal L}$ to express the wave field vector ${\bf q}({\bf x})$ in terms of downgoing and upgoing waves $p^{+}({\bf x})$ and $p^{-}({\bf x})$ via

[TABLE]

with

[TABLE]

Note that these equations imply

[TABLE]

hence, the downgoing and upgoing waves $p^{+}$ and $p^{-}$ are pressure-normalised. An advantage of pressure-normalised (or, more generally, field-normalised) decomposition is that the decomposed quantities simply add up to a field quantity (acoustic pressure in the case of equation (24)). This property does not apply to flux-normalised decomposed wave fields (Frasier, 1970; Kennett et al., 1978; Ursin, 1983). On the other hand, an advantage of flux-normalised decomposition is that the underlying equations obey more simple symmetry properties. For a comprehensive discussion on field-normalised versus flux-normalised decomposition in inhomogeneous media, see de Hoop (1996) and Wapenaar (2020). In this paper we use field-normalised decomposition. In the remainder of section 2 we apply decomposition only at and above $\partial\mathbb{D}_{F}$, where the medium is assumed to be homogeneous. In section 3 we will apply decomposition also inside the inhomogeneous medium.

We use equation (22) at $\partial\mathbb{D}_{F}$ to derive focusing functions and express them in the components of the propagator matrix and vice-versa. Substituting equation (22), with ${\bf x}$ replaced by ${\bf x}_{F}$, into the right-hand side of equation (10) gives

[TABLE]

for $x_{3}\geq x_{3,F}$, with

[TABLE]

or

[TABLE]

The operators $\pm\frac{1}{\omega\rho_{0}}{\cal H}_{1}({\bf x}_{F})$ in equation (27) act, via equation (25), on $p^{\pm}({\bf x}_{F})$. However, since these operators are symmetric (in the sense of equation (21)), we may replace the actions of these operators on $p^{\pm}({\bf x}_{F})$ by actions on the elements $W^{p,v}({\bf x},{\bf x}_{F})$ and $W^{v,v}({\bf x},{\bf x}_{F})$. To be more specific, if we partition ${\bf Y}({\bf x},{\bf x}_{F})$ as follows

[TABLE]

we obtain from equation (27) for the elements of this matrix

[TABLE]

We analyse these expressions one by one. First we consider the element $Y^{p,-}$. From equations (25) and (28) it can be seen that the superscript $p$ refers to the acoustic pressure $p({\bf x})$ contained in ${\bf q}({\bf x})$ and superscript $-$ refers to the upgoing wave field component $p^{-}({\bf x}_{F})$ in ${\bf p}({\bf x}_{F})$. Using equations (9), (11) and (29) we obtain

[TABLE]

which is a focusing condition. Hence, we define

[TABLE]

with $F^{p}({\bf x},{\bf x}_{F})$ denoting a focusing function for the acoustic pressure $p$, which focuses as an upgoing field at ${\bf x}={\bf x}_{F}$ and continues as an upgoing field in the homogeneous upper half-space, see the lower frame of Figure 1. Next, we consider the element $Y^{v,-}$. Superscript $v$ refers to the vertical particle velocity $v_{3}({\bf x})$ contained in ${\bf q}({\bf x})$ and superscript $-$ refers again to the upgoing wavefield component $p^{-}({\bf x}_{F})$ in ${\bf p}({\bf x}_{F})$. Using equations (9), (11) and (30) we obtain

[TABLE]

which is also a focusing condition, but somewhat more complicated than equation (31) because of the mix of the involved wavefield components $v_{3}({\bf x})$ and $p^{-}({\bf x}_{F})$. Hence, we define

[TABLE]

with $F^{v}({\bf x},{\bf x}_{F})$ denoting the particle velocity counterpart of the focusing function $F^{p}({\bf x},{\bf x}_{F})$ (note that the definition of $F^{v}({\bf x},{\bf x}_{F})$ is different from that in Wapenaar (2022), to facilitate the derivations below). The focusing functions $F^{p}({\bf x},{\bf x}_{F})$ and $F^{v}({\bf x},{\bf x}_{F})$, which together form the right column of matrix ${\bf Y}({\bf x},{\bf x}_{F})$, are illustrated in the lower frame of Figure 1. They resemble the focusing function $f_{2}$ introduced in previous work (Wapenaar et al., 2013; Slob et al., 2014), which also focuses at the upper boundary (as opposed to the focusing function $f_{1}$, which focuses inside the medium). However, there are also some notable differences. First, $f_{2}({\bf x},{\bf x}_{F})$ is defined in a truncated version of the actual medium and is obtained from a superposition of downgoing and upgoing components, $f_{2}^{+}({\bf x},{\bf x}_{F})$ and $f_{2}^{-}({\bf x},{\bf x}_{F})$ respectively, at ${\bf x}$ inside the medium (at the lower boundary of the truncated medium). Moreover, representations involving $f_{2}^{+}$ and $f_{2}^{-}$ ignore evanescent waves at $x_{3,F}$ and $x_{3}$. In contrast, $F^{p}({\bf x},{\bf x}_{F})$ and $F^{v}({\bf x},{\bf x}_{F})$ are defined in the actual (i.e., untruncated) medium and represent the full pressure and vertical particle velocity at ${\bf x}$ of a field that focuses at ${\bf x}_{F}$ at the upper boundary. Since they are derived from the propagator matrix, these focusing functions account for evanescent waves (this will be demonstrated below with a numerical example). The only decomposition takes place at the boundary $\partial\mathbb{D}_{F}$, where the medium is homogeneous. This decomposition, formulated by equations (32) and (34), accounts for evanescent waves. Last but not least, $F^{p}$ and $F^{v}$ hold for dissipative media and they are normalised differently from $f_{2}$.

Before we analyse the elements in the left column of matrix ${\bf Y}({\bf x},{\bf x}_{F})$, we introduce an adjoint medium, with parameters $\bar{\kappa}({\bf x})=\kappa^{*}({\bf x})$ and $\bar{\rho}({\bf x})=\rho^{*}({\bf x})$. The bar denotes the adjoint medium and the superscript asterisk denotes complex conjugation. When the original medium is dissipative, the adjoint medium is effectual, with (for positive $\omega$) $\Im(\bar{\kappa})\leq 0$ and $\Im(\bar{\rho})\leq 0$. Waves propagating through an effectual medium gain energy (Bojarski, 1983; de Hoop, 1988). Adjoint media are usually associated to a computational state. The operator matrix {\,\,\,\bar{\mbox{\!\!\!\boldmath{\cal A}}}} and the Helmholtz operator $\bar{\cal H}_{2}$ of the adjoint medium are defined similarly as ${\cal A}$ and ${\cal H}_{2}$ in equations (4) and (5), respectively, but with $\kappa({\bf x})$ and $\rho({\bf x})$ replaced by $\bar{\kappa}({\bf x})$ and $\bar{\rho}({\bf x})$, respectively. Hence, $\bar{\cal H}_{2}={\cal H}_{2}^{*}$. Analogous to equations (8) and (9), we define the propagator matrix $\bar{\bf W}({\bf x},{\bf x}_{F})$ of the adjoint medium as the solution of \partial_{3}\bar{\bf W}({\bf x},{\bf x}_{F})={\,\,\,\bar{\mbox{\!\!\!\boldmath{\cal A}}}}({\bf x})\bar{\bf W}({\bf x},{\bf x}_{F}), with boundary condition $\bar{\bf W}({\bf x},{\bf x}_{F})|_{x_{3}={x_{3,F}}}={\bf I}\delta({{\bf x}_{\rm H}}-{{\bf x}_{{\rm H},F}})$. In Appendix C we derive

[TABLE]

(equation (146)). For the square-root operator we have, similar as for the Helmholtz operator,

[TABLE]

(Wapenaar et al., 2001). Using equations (35) and (36) in equations (29) and (30), we find $\bar{Y}^{p,+}({\bf x},{\bf x}_{F})=Y^{p,-*}({\bf x},{\bf x}_{F})$ and $\bar{Y}^{v,+}({\bf x},{\bf x}_{F})=-Y^{v,-*}({\bf x},{\bf x}_{F})$. Hence, using equations (32) and (34), we find for the elements in the left column of matrix ${\bf Y}({\bf x},{\bf x}_{F})$

[TABLE]

For matrix ${\bf Y}({\bf x},{\bf x}_{F})$ we thus obtain

[TABLE]

Note that $F^{p}$, $F^{v}$, $\bar{F}^{p*}$ and $\bar{F}^{v*}$ are expressed in terms of the components of the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ via equations (32), (34), (37) and (38). Conversely, we can express the components of the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ in terms of the focusing functions $F^{p}$, $F^{v}$, $\bar{F}^{p*}$ and $\bar{F}^{v*}$. Inverting equation (26) yields

[TABLE]

with {{{{\mbox{\boldmath{\cal L}}}}}}^{-1} defined in equation (20). Since operator $\frac{1}{\rho}{\cal H}_{1}$ is symmetric, its inverse ${\cal H}_{1}^{-1}\rho$ is symmetric as well. Hence, in equation (40) these operators can be taken to act on the elements of matrix ${\bf Y}({\bf x},{\bf x}_{F})$. This yields

[TABLE]

Note that up to this point the medium may be dissipative (and its adjoint effectual), and evanescent wave modes are accounted for, inside the medium as well as at the boundary $\partial\mathbb{D}_{F}$. Hence, the expressions in this section are more general than their counterparts in Wapenaar (2022), which were derived for a lossless medium, under the assumption that evanescent waves can be ignored at $\partial\mathbb{D}_{F}$. If we make the same assumptions here, we can omit the bars on $F^{p}$ and $F^{v}$. For this situation equations (41) – (44) simplify to

[TABLE]

We illustrate the focusing function and its relation with the propagator matrix with a numerical example. Applying the transformations of equations (14) and (16) to equation (32) (assuming a laterally invariant medium), taking ${\bf x}_{F}=(0,0,x_{3,F})$ and $s_{2}=0$, we obtain

[TABLE]

with vertical slowness $s_{3,0}=\sqrt{1/c_{0}^{2}-s_{1}^{2}}$ being the spatial Fourier transform of $\frac{1}{\omega}{\cal H}_{1}$ at $x_{3,F}$ for the laterally invariant medium (here we assumed $s_{1}^{2}<1/c_{0}^{2}$). Equation (49) shows how a weighted superposition of the even component $W^{p,p}$ of Figure 2(b) and the odd component $W^{p,v}$ of Figure 2(c) yields the focusing function $F^{p}(s_{1},x_{3},x_{3,F},\tau)$. This focusing function is shown in Figure 3(a) for $s_{1}=1/3500$ m/s. The upper trace at $x_{3}=x_{3,F}=0$ m represents the focusing condition $F^{p}(s_{1},x_{3,F},x_{3,F},\tau)=\delta(\tau)$. At and above $x_{3,F}$ the focusing function is an upgoing field. Note that, similar as in Figure 2, the wave field tunnels through the high-velocity layer between $x_{3}=760$ m and $x_{3}=800$ m, which confirms that this focusing function accounts for evanescent waves inside the medium. The time-reversed focusing function $F^{p}(s_{1},x_{3},x_{3,F},-\tau)$ is shown in Figure 3(b). The focusing function of Figure 3(a) and its time-reversed version of Figure 3(b) can be combined to give components of the propagator matrix. To this end, equations (45) and (46) are transformed to (assuming $s_{1}^{2}<1/c_{0}^{2}$)

[TABLE]

For the acoustic case all components are symmetric in $s_{1}$, i.e., $F^{p}(-s_{1},x_{3},x_{3,F},-\tau)=F^{p}(s_{1},x_{3},x_{3,F},-\tau)$, etc. Hence, equations (50) and (51) show how the even and odd components $W^{p,p}(s_{1},x_{3},x_{3,F},\tau)$ and $W^{p,v}(s_{1},x_{3},x_{3,F},\tau)$ of Figures 2(b) and 2(c) are obtained from the focusing function $F^{p}(s_{1},x_{3},x_{3,F},\tau)$ and its time-reversal $F^{p}(s_{1},x_{3},x_{3,F},-\tau)$ of Figure 3.

2.4 Representations with acoustic Marchenko focusing functions

Substituting the expressions for ${\bf q}({\bf x})$, ${\bf p}({\bf x}_{F})$ and ${\bf Y}({\bf x},{\bf x}_{F})$ (equations (2), (23) and (39)) into equation (25), gives the following representations for the acoustic pressure $p({\bf x})$ and the vertical particle velocity $v_{3}({\bf x})$ inside the inhomogeneous medium

[TABLE]

for $x_{3}\geq x_{3,F}$. These expressions are exact and hold for dissipative media. Equation (52) is a generalisation of equation (17) of Wapenaar & de Ridder (2022) for dissipative media.

We use equations (52) and (53) to derive representations for Green’s functions between the boundary $\partial\mathbb{D}_{F}$ and any position ${\bf x}$ inside the medium. To this end, we define a unit point source of vertical force at ${\bf x}_{S}$ just above $\partial\mathbb{D}_{F}$. For the downgoing field at $\partial\mathbb{D}_{F}$ (i.e., just below the source), we then have $p^{+}({\bf x}_{F})=\frac{1}{2}\delta({\bf x}_{{\rm H},F}-{\bf x}_{{\rm H},S})$, where ${\bf x}_{{\rm H},S}$ denotes the horizontal coordinates of ${\bf x}_{S}$. The upgoing field at $\partial\mathbb{D}_{F}$ is the reflection response to this downgoing source field, hence $p^{-}({\bf x}_{F})=\frac{1}{2}R({\bf x}_{F},{\bf x}_{S})$. The field at ${\bf x}$ inside the medium is the Green’s response to the source at ${\bf x}_{S}$, hence $p({\bf x})=G^{p,f}({\bf x},{\bf x}_{S})$ and $v_{3}({\bf x})=G^{v,f}({\bf x},{\bf x}_{S})$. Here the second superscript ($f$) refers to the vertical force source at ${\bf x}_{S}$, whereas the first superscripts ($p$ and $v$) refer to the observed quantities (pressure and vertical particle velocity) at ${\bf x}$. Substitution of these expressions for $p^{\pm}({\bf x}_{F})$, $p({\bf x})$ and $v_{3}({\bf x})$ into equations (52) and (53) gives

[TABLE]

for $x_{3}\geq x_{3,F}$. Slob (2016) derived similar representations for decomposed wave fields in dissipative media. In the present derivation we only used decomposition at the boundary $\partial\mathbb{D}_{F}$ (similar as Diekmann & Vasconcelos (2021, 2023) and Wapenaar et al. (2021)). This implies that inside the medium the wavefield does not need to be decomposed into downgoing and upgoing waves and that evanescent waves can be present.

When the medium is lossless and evanescent waves are neglected at $\partial\mathbb{D}_{F}$, the bars on $F^{p}$ and $F^{v}$ in representations (54) and (55) can be omitted. Using the Marchenko method, these focusing functions can then be retrieved from the reflection response $R({\bf x}_{F},{\bf x}_{S})$ and a smooth background model (Wapenaar et al., 2013; Elison et al., 2020). Since representations (54) and (55) account for evanescent waves inside the medium, the retrieved focusing functions potentially also account for evanescent waves inside the medium (this is subject of current research (Brackenhoff & Wapenaar, 2023)). Once the focusing functions are found, they can be used to retrieve the Green’s functions $G^{p,f}({\bf x},{\bf x}_{S})$ and $G^{v,f}({\bf x},{\bf x}_{S})$ (from equations (54) and (55)) and all components of the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ (from equations (45) – (48)).

3 Acoustic transfer matrix and decomposed focusing functions

3.1 Acoustic transfer matrix

We introduce the transfer matrix as follows. Given the downgoing and upgoing fields $p^{+}({\bf x}_{F})$ and $p^{-}({\bf x}_{F})$ at the boundary $\partial\mathbb{D}_{F}$, we ‘transfer’ these fields to downgoing and upgoing fields $p^{+}({\bf x})$ and $p^{-}({\bf x})$ at any depth level $x_{3}$ inside the medium using the following expression

[TABLE]

for $x_{3}\geq x_{3,F}$. Vectors ${\bf p}({\bf x}_{F})$ and ${\bf p}({\bf x})$ contain the downgoing and upgoing fields at depth levels $x_{3,F}$ and $x_{3}$ (equation (23)). We call {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) the transfer matrix, which we partition as follows

[TABLE]

For each component of this matrix, the superscripts refer to the propagation direction at ${\bf x}$ and at ${\bf x}_{F}$, respectively. Equation (56) is illustrated in the upper-right frame of Figure 1.

For horizontally layered media, the transfer matrix is usually built up recursively from interface to interface (Born & Wolf, 1965; Katsidis & Siapkas, 2002; Elison, 2020; Dukalski et al., 2022a, b). Here we follow a different approach to derive an expression for {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) for laterally varying media. Substituting equation (22) into equation (10) we obtain equation (56), with

[TABLE]

with {{{{\mbox{\boldmath{\cal L}}}}}}({\bf x}) and its inverse defined in equation (20). Equation (58), which relates the transfer matrix {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) to the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$, is illustrated in the upper half of Figure 1. In the next section we show that the transfer matrix can be expressed in terms of decomposed Marchenko focusing functions.

3.2 Relation between acoustic transfer matrix and decomposed Marchenko focusing functions

From equations (26) and (58) we find

[TABLE]

According to equation (39), the right column of ${\bf Y}({\bf x},{\bf x}_{F})$ contains $F^{p}({\bf x},{\bf x}_{F})$ and $F^{v}({\bf x},{\bf x}_{F})$, i.e., the pressure and vertical particle velocity components at ${\bf x}$ of the focusing function. Hence, analogous to {\bf p}({\bf x})={{{\mbox{\boldmath{\cal L}}}}}^{-1}({\bf x}){\bf q}({\bf x}), we obtain for the right column of {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F})

[TABLE]

with $F^{+}({\bf x},{\bf x}_{F})$ and $F^{-}({\bf x},{\bf x}_{F})$ being the downgoing and upgoing parts at ${\bf x}$ of the focusing function $F^{p}({\bf x},{\bf x}_{F})$.

According to equation (39), the left column of ${\bf Y}({\bf x},{\bf x}_{F})$ contains $\bar{F}^{p*}({\bf x},{\bf x}_{F})$ and $-\bar{F}^{v*}({\bf x},{\bf x}_{F})$. Hence, for the left column of {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) we obtain

[TABLE]

or, using ${\cal H}_{1}=\bar{\cal H}_{1}^{*}$ (equation (36)) and $\rho=\bar{\rho}^{*}$,

[TABLE]

Comparing this with equation (60) we find that this gives a vector with $\bar{F}^{-*}({\bf x},{\bf x}_{F})$ and $\bar{F}^{+*}({\bf x},{\bf x}_{F})$. This is the left column of {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}). Hence, we have obtained

[TABLE]

see the upper-right frame of Figure 1. Hence, the transfer matrix for an inhomogeneous dissipative acoustic medium is expressed in terms of decomposed focusing functions of the medium and its adjoint.

We consider the special case of a horizontally layered medium. Applying the transformations of equations (14) and (16) to equation (63), taking ${\bf x}_{F}=(0,0,x_{3,F})$, we obtain

[TABLE]

Dukalski et al. (2022a, b) used a recursive approach and obtained an expression similar to equation (64). In their derivation they used a path-reversal operator ${\cal P}$, which is equivalent with (i) taking the adjoint medium, (ii) taking the complex conjugate (or in the time domain taking the time-reversal) and (iii) changing the sign of the horizontal slowness. Hence, ${\cal P}\{F^{\pm}({\bf s},x_{3},x_{3,F},\tau)\}$ is equivalent with $\bar{F}^{\pm}(-{\bf s},x_{3},x_{3,F},-\tau)$.

For the lossless medium of Figure 2(a), the decomposed focusing functions $F^{-}(s_{1},x_{3},x_{3,F},\tau)$ and $F^{+}(s_{1},x_{3},x_{3,F},\tau)$ for $s_{1}=1/3500$ m/s and $s_{2}=0$ are shown in Figures 4(a) and 4(b), respectively. For each $x_{3}$, the function $F^{-}(s_{1},x_{3},x_{3,F},\tau)$ can be seen as the intricate field that needs to be emitted upward from $x_{3}$ to arrive as a single upward propagating field at the focal depth $x_{3,F}$ at $\tau=0$. For the same $x_{3}$, the function $F^{+}(s_{1},x_{3},x_{3,F},\tau)$ is the downward reflected response to $F^{-}(s_{1},x_{3},x_{3,F},\tau)$. Figures 4(b) and 4(a) together form the right column of the transformed transfer matrix {{{\mbox{\boldmath{\cal T}}}}}(s_{1},x_{3},x_{3,F},\tau). Their superposition gives the focusing function $F^{p}(s_{1},x_{3},x_{3,F},\tau)$, shown in Figure 3(a).

3.3 Representations with decomposed acoustic Marchenko focusing functions

Substituting the expressions for ${\bf p}({\bf x})$ and {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) (equations (23) and (63)) into equation (56), gives the following representations for the downgoing and upgoing components of the acoustic pressure, $p^{+}({\bf x})$ and $p^{-}({\bf x})$ respectively, inside the inhomogeneous medium

[TABLE]

for $x_{3}\geq x_{3,F}$. These expressions are exact and hold for dissipative media. Making similar substitutions as in section 2.4 we obtain

[TABLE]

for $x_{3}\geq x_{3,F}$. Here $G^{\pm,f}({\bf x},{\bf x}_{S})$ stands for the downgoing ($+$) and upgoing ($-$) part of the Green’s function $G^{p,f}({\bf x},{\bf x}_{S})$.

When the medium is lossless and when evanescent waves are neglected at $\partial\mathbb{D}_{F}$ and at depth level $x_{3}$ inside the medium, the bars on $F^{+}$ and $F^{-}$ in representations (67) and (68) can be omitted. Using the Marchenko method, these decomposed focusing functions can then be retrieved from the reflection response $R({\bf x}_{F},{\bf x}_{S})$ and a smooth background model (Wapenaar et al., 2013; Slob et al., 2014). Once the focusing functions are found, they can be used to retrieve the decomposed Green’s functions $G^{+,f}({\bf x},{\bf x}_{S})$ and $G^{-,f}({\bf x},{\bf x}_{S})$ (from equations (67) and (68)) and all components of the transfer matrix {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) (from equation (63)).

4 Conclusions

We have derived relations between acoustic propagator matrices, transfer matrices and Marchenko focusing functions. In the appendices we generalize the expressions for other wave phenomena. All relations hold for a heterogeneous dissipative medium below a homogeneous upper half-space and account for propagating and evanescent waves. Only for the transfer matrix beyond the acoustic approach (Appendix B) we assume that there are no lateral variations at the depth level inside the medium where decomposition takes place.

The derived relations provide insight in the connections between the propagator matrices, transfer matrices and Marchenko focusing functions and may lead to new modelling algorithms for these quantities. Moreover, several of the derived relations may be useful to develop improved Marchenko-type wave field retrieval and imaging schemes for different wave phenomena, possibly accounting for evanescent waves inside the medium.

Acknowledgements

We thank Ivan Vasconcelos and an anonymous reviewer for their useful comments, which helped us to improve this paper. The research of KW has received funding from the European Union’s Horizon 2020 research and innovation programme: European Research Council (grant agreement 742703).

Data Availability

No data have been used for this study.

Appendix A Unified propagator matrix and focusing functions

In this appendix, we extend the theory of section 2 to unified wave fields.

A.1 Unified matrix-vector wave equation

Consider matrix-vector wave equation (1). We partition the $N\times 1$ wave field vector ${\bf q}$, the $N\times 1$ source vector ${\bf d}$ and the $N\times N$ operator matrix ${\cal A}$ as follows

[TABLE]

This includes the acoustic situation (for $N=2$) discussed in section 2. The vectors and operator matrix for other wave phenomena can be found in various references (Woodhouse, 1974; Ursin, 1983; Van Stralen, 1997; Gelinsky & Shapiro, 1997; Haartsen & Pride, 1997; White & Zhou, 2006; Løseth & Ursin, 2007). A comprehensive overview is given by Wapenaar (2019) for acoustic waves ($N=2$), quantum mechanical waves obeying the Schrödinger equation ($N=2$), electromagnetic waves ($N=4$), elastodynamic waves ($N=6$), poroelastodynamic waves ($N=8$), piezoelectric waves ($N=10$) and seismoelectric waves ($N=12$). For all these wave phenomena, the operator matrix ${\cal A}$ obeys the following symmetry properties

[TABLE]

with

[TABLE]

where ${\bf O}$ and ${\bf I}$ are zero and identity matrices of appropriate size. Superscript $t$ denotes transposition of the matrix and the operators contained in it, with $\partial_{\alpha}^{t}=-\partial_{\alpha}$. Superscript $\dagger$ denotes transposition and complex conjugation. The bar on an operator denotes that it is defined in the adjoint medium. For further details we refer to the aforementioned references.

We show that, with some modifications, equations (69) – (72) also hold for the matrix-vector form of the Dirac equation for the $4\times 1$ spinor $\psi$. The Dirac equation is given by (Sakurai, 1967)

[TABLE]

(summation over $\mu$ from $1$ to $4$), with the Dirac spinor partitioned as

[TABLE]

and with

[TABLE]

with $V({\bf x})$ the space-dependent potential, $h$ Planck’s constant, $c$ the speed of light, $e$ the electron charge and \mbox{\boldmath\sigma}_{1}, \mbox{\boldmath\sigma}_{2} and \mbox{\boldmath\sigma}_{3} the Pauli matrices, defined as

[TABLE]

Assuming a time-dependence $\exp(-iEt/\hbar)$, we replace $\partial_{4}$ by $-(E-eV)/\hbar c$. Equation (74) can be rewritten as follows

[TABLE]

or

[TABLE]

Pre-multiplying all terms by -i\mbox{\boldmath\sigma}_{3} and bringing the $\partial_{3}$-operators to the left-hand side, using \mbox{\boldmath\sigma}_{3}\mbox{\boldmath\sigma}_{1}=i\mbox{\boldmath\sigma}_{2}, \mbox{\boldmath\sigma}_{3}\mbox{\boldmath\sigma}_{2}=-i\mbox{\boldmath\sigma}_{1} and \mbox{\boldmath\sigma}_{3}\mbox{\boldmath\sigma}_{3}={\bf I}, yields a set of equations which can be recast in the form of equation (1), with $4\times 1$ vectors ${\bf q}$ and ${\bf d}$ and $4\times 4$ operator matrix ${\cal A}$ partitioned as in equation (69), with

[TABLE]

With these definitions of the operator submatrices, matrix ${\cal A}$ obeys symmetry relations (70) – (72), with ${\bf N}$, ${\bf K}$ and ${\bf J}$ defined as

[TABLE]

Although there are no direct applications for geophysics, the Schrödinger and Dirac equations are included in all derivations below, since this comes almost for free. When we speak of the ‘medium’, for the Schrödinger and Dirac equations it should be understood as the ‘potential’.

A.2 Unified propagator matrix

We define the unified $N\times N$ propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ as the solution of wave equation (8), with boundary condition (9) and with the operator matrix ${\cal A}$ being the unified operator matrix discussed in Appendix A.1. Using equation (10), the unified wave field vector ${\bf q}({\bf x})$ can be propagated from $x_{3,F}$ to any depth level $x_{3}$, assuming there are no sources between these depth levels. We partition ${\bf W}({\bf x},{\bf x}_{F})$ as

[TABLE]

where ${\bf W}_{11}$, ${\bf W}_{12}$, ${\bf W}_{21}$ and ${\bf W}_{22}$ are $\frac{N}{2}\times\frac{N}{2}$ submatrices of ${\bf W}$. For each of these submatrices, the second subscript refers to the wave field component (${\bf q}_{1}$ or ${\bf q}_{2}$) it acts on at ${\bf x}_{F}$, whereas the first subscript refers to the wave field component it contributes to at ${\bf x}$. ${\bf W}({\bf x},{\bf x}_{F})$ obeys the recursive relation (12), and from equation (13) it follows that ${\bf W}({\bf x}_{F},{\bf x})$ is the inverse of ${\bf W}({\bf x},{\bf x}_{F})$.

A.3 Relation between unified propagator matrix and Marchenko focusing functions

We assume again that the medium at and above $\partial\mathbb{D}_{F}$ is homogeneous and may be dissipative. The medium below $\partial\mathbb{D}_{F}$ may be inhomogeneous and dissipative, and it is source-free. In preparation for defining the focusing functions, in the upper half-space (i.e., at and above $\partial\mathbb{D}_{F}$) we apply eigenvalue decomposition to matrix $\tilde{\bf A}$ (the spatial Fourier transform of operator matrix ${\cal A}$), as follows

[TABLE]

with

[TABLE]

with ${\bf S}_{3}^{+}$ and ${\bf S}_{3}^{-}$ being diagonal matrices containing vertical slownesses for downgoing and upgoing waves, respectively. We express the Fourier transformed wave field vector $\tilde{\bf q}({\bf s},x_{3})$ in terms of downgoing and upgoing waves $\tilde{\bf p}^{+}({\bf s},x_{3})$ and $\tilde{\bf p}^{-}({\bf s},x_{3})$ via

[TABLE]

with

[TABLE]

Note that these equations imply $\tilde{\bf q}_{1}=\tilde{\bf L}_{1}^{+}\tilde{\bf p}^{+}+\tilde{\bf L}_{1}^{-}\tilde{\bf p}^{-}$. Similar as in section 2.3, equation (24), we continue with downgoing and upgoing waves $\tilde{\bf q}_{1}^{+}$ and $\tilde{\bf q}_{1}^{-}$ which are field-normalised such that $\tilde{\bf q}_{1}=\tilde{\bf q}_{1}^{+}+\tilde{\bf q}_{1}^{-}$. To this end, we define $\tilde{\bf q}_{1}^{\pm}=\tilde{\bf L}_{1}^{\pm}\tilde{\bf p}^{\pm}$ and we replace equation (91) by

[TABLE]

where

[TABLE]

with

[TABLE]

Whereas there is ambiguity in the normalization of the matrices $\tilde{\bf L}_{1}^{\pm}$ and $\tilde{\bf L}_{2}^{\pm}$, the matrix $\tilde{\bf D}_{1}^{\pm}$ is uniquely defined (for each wave phenomenon). Note that for the acoustic situation we have $\tilde{\bf L}_{1}^{\pm}=1$, hence $\tilde{\bf D}_{1}^{\pm}=\tilde{\bf L}_{2}^{\pm}$, $\tilde{\bf D}=\tilde{\bf L}$ and $\tilde{\bf b}=\tilde{\bf p}$. Some other examples of matrix $\tilde{\bf D}_{1}^{\pm}$ (for electromagnetic and elastodynamic waves) are given by Wapenaar (2022). In Appendix C we derive for any wave phenomenon

[TABLE]

with ${\bf J}_{11}$ and ${\bf J}_{22}$ being the $\frac{N}{2}\times\frac{N}{2}$ submatrices of $N\times N$ matrix ${\bf J}$. From equation (73) we have for all wave phenomena except for the Dirac equation ${\bf J}_{11}=-{\bf J}_{22}={\bf I}$, and from equation (87) we have for the Dirac equation {\bf J}_{11}={\bf J}_{22}=\mbox{\boldmath\sigma}_{2}.

We use equation (93) at $\partial\mathbb{D}_{F}$ and the properties of matrix $\tilde{\bar{{\bf D}}}_{1}^{\pm}({\bf s},x_{3})$ to derive unified focusing functions and express them in the components of the unified propagator matrix and vice-versa. First we aim to substitute equation (93) for $x_{3}=x_{3,F}$ into a transformed version of equation (10). This equation contains the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$. For a function of two space variables, $u({{\bf x}},{\bf x}_{F})$ (with ${\bf x}_{F}$ at $\partial\mathbb{D}_{F}$), we define the spatial Fourier transformation along the horizontal components of the second space variable as

[TABLE]

and its inverse as

[TABLE]

Note that the sign in the exponential of equation (98) is opposite to that in equation (14). Using these definitions and Parseval’s theorem, we rewrite equation (10) as

[TABLE]

with $\tilde{\bf W}({{\bf x}},{{\bf s}},{x_{3,F}})$ obeying the boundary condition

[TABLE]

Substitution of equation (93) into equation (100) gives

[TABLE]

for $x_{3}\geq x_{3,F}$, with

[TABLE]

We partition matrix $\tilde{\bf Y}({{\bf x}},{{\bf s}},{x_{3,F}})$ as follows

[TABLE]

Using equation (94) and the spatial Fourier transformation of equation (88), we obtain

[TABLE]

We analyse these expressions one by one. First consider $\tilde{\bf Y}_{1}^{-}({{\bf x}},{{\bf s}},{x_{3,F}})$. Via equation (102) it can be seen that subscript $1$ refers to wavefield component ${\bf q}_{1}$ at ${\bf x}$ and superscript $-$ refers to the upgoing wavefield component $\tilde{\bf q}_{1}^{-}$ at $x_{3,F}$. Moreover, for $x_{3}=x_{3,F}$ we obtain, using equation (101), $\tilde{\bf Y}_{1}^{-}({{\bf x}},{{\bf s}},{x_{3,F}})|_{x_{3}={x_{3,F}}}={\bf I}\exp\{i\omega{{\bf s}}\cdot{{\bf x}_{{\rm H}}}\}$, or, applying the inverse spatial Fourier transformation defined in equation (99), ${\bf Y}_{1}^{-}({\bf x},{\bf x}_{F})|_{x_{3}={x_{3,F}}}={\bf I}\delta({{\bf x}_{\rm H}}-{{\bf x}_{{\rm H},F}})$, which is a focusing condition. Hence, we define

[TABLE]

with $\tilde{\bf F}_{1}({{\bf x}},{{\bf s}},{x_{3,F}})$ denoting the spatial Fourier transform of the focusing function ${\bf F}_{1}({\bf x},{\bf x}_{F})$ for wavefield component ${\bf q}_{1}$, which focuses as an upgoing field at ${\bf x}={\bf x}_{F}$ and continues as an upgoing field in the homogeneous upper half-space. Note that the focusing function is a $\frac{N}{2}\times\frac{N}{2}$ matrix. Next, we consider $\tilde{\bf Y}_{2}^{-}({{\bf x}},{{\bf s}},{x_{3,F}})$. Subscript $2$ refers to wavefield component ${\bf q}_{2}$ at ${\bf x}$ and superscript $-$ refers again to the upgoing wavefield component $\tilde{\bf q}_{1}^{-}$ at $x_{3,F}$. For $x_{3}=x_{3,F}$ we obtain, using equation (101), $\tilde{\bf Y}_{2}^{-}({{\bf x}},{{\bf s}},{x_{3,F}})|_{x_{3}={x_{3,F}}}=\tilde{\bf D}_{1}^{-}({\bf s},x_{3,F})\exp\{i\omega{{\bf s}}\cdot{{\bf x}_{{\rm H}}}\}$, which is a focusing condition, but somewhat more complicated than for $\tilde{\bf Y}_{1}^{-}({{\bf x}},{{\bf s}},{x_{3,F}})$ because of the mix of wavefield components ${\bf q}_{2}$ and $\tilde{\bf q}_{1}^{-}$. Hence, we define

[TABLE]

with $\tilde{\bf F}_{2}({{\bf x}},{{\bf s}},{x_{3,F}})$ denoting the spatial Fourier transform of the focusing function ${\bf F}_{2}({\bf x},{\bf x}_{F})$ for wavefield component ${\bf q}_{2}$, which focuses as an upgoing field at ${\bf x}={\bf x}_{F}$ and continues as an upgoing field in the homogeneous upper half-space (note that the definition of $\tilde{\bf F}_{2}$ is different from that in Wapenaar (2022), to facilitate the derivations below). The focusing functions $\tilde{\bf F}_{1}({{\bf x}},{{\bf s}},{x_{3,F}})$ and $\tilde{\bf F}_{2}({{\bf x}},{{\bf s}},{x_{3,F}})$ together form the right column of matrix $\tilde{\bf Y}({{\bf x}},{{\bf s}},{x_{3,F}})$.

For the analysis of the submatrices in the left column of $\tilde{\bf Y}({{\bf x}},{{\bf s}},{x_{3,F}})$, we use symmetry relation (97) and we need a similar relation for the submatrices of $\tilde{\bf W}({{\bf x}},{{\bf s}},{x_{3,F}})$. In Appendix C we derive

[TABLE]

From the spatial Fourier transform of this equation we obtain for the submatrices of $\tilde{\bf W}({{\bf x}},{{\bf s}},{x_{3,F}})$

[TABLE]

(no summation for repeated subscripts). Substituting equations (97) and (110) into equations (105) and (106) yields

[TABLE]

Hence, using equations (107) and (108), we find for the submatrices in the left column of $\tilde{\bf Y}({{\bf x}},{{\bf s}},{x_{3,F}})$

[TABLE]

Hence, matrix $\tilde{\bf Y}({{\bf x}},{{\bf s}},{x_{3,F}})$ becomes

[TABLE]

or, using the inverse Fourier transformation defined in equation (99),

[TABLE]

This is a generalisation of equation (39). Note that $\tilde{\bf F}_{1}$, $\tilde{\bf F}_{2}$, ${\tilde{\bar{\bf F}}}_{1}^{*}$ and ${\tilde{\bar{\bf F}}}_{2}^{*}$ are expressed in terms of the submatrices of the propagator matrix $\tilde{\bf W}({{\bf x}},{{\bf s}},{x_{3,F}})$ via equations (107), (108), (113) and (114). Conversely, we can express the submatrices of the propagator matrix $\tilde{\bf W}({{\bf x}},{{\bf s}},{x_{3,F}})$ in terms of the focusing functions $\tilde{\bf F}_{1}$, $\tilde{\bf F}_{2}$, ${\tilde{\bar{\bf F}}}_{1}^{*}$ and ${\tilde{\bar{\bf F}}}_{2}^{*}$. To this end, we start with inverting equation (103), according to

[TABLE]

with

[TABLE]

[TABLE]

Using equations (115) and (118), we obtain

[TABLE]

(no summation for repeated subscripts). These expressions are a generalisation of equations (41) – (44). Those equations follow as a special case from equations (120) and (121) by substituting ${\bf J}_{11}=-{\bf J}_{22}=1$, $\tilde{\bf D}_{1}^{\pm}({{\bf s}},x_{3,F})=\pm{s_{3,0}}/{\rho_{0}}$, $\{\tilde{\bf\Delta}_{1}({{\bf s}},x_{3,F})\}^{-1}={\rho_{0}}/{2s_{3,0}}$, and applying an inverse spatial Fourier transformation, which involves replacing $s_{3,0}$ by operator $\frac{1}{\omega}{\cal H}_{1}({\bf x}_{F})$.

A.4 Representations with unified Marchenko focusing functions

Applying Parseval’s theorem to equation (102) and substituting the expressions for ${\bf q}({\bf x})$, ${\bf b}({\bf x}_{F})$ and ${\bf Y}({\bf x},{\bf x}_{F})$ (equations (69), (95) and (116)), gives the following representation for the quantities ${\bf q}_{1}({\bf x})$ and ${\bf q}_{2}({\bf x})$ inside the inhomogeneous medium

[TABLE]

(no summation for repeated subscripts) for $x_{3}\geq x_{3,F}$. This is a generalisation of equations (52) and (53).

We use equation (122) to derive representations for Green’s functions between the boundary $\partial\mathbb{D}_{F}$ and any position ${\bf x}$ inside the medium. We define a unit ${\bf d}_{2}$-type source (see equation (69)) at ${\bf x}_{S}$ just above $\partial\mathbb{D}_{F}$. The $\frac{N}{2}\times\frac{N}{2}$ Green’s matrix ${\bf G}_{12}({\bf x},{\bf x}_{S})$ stands for the ${\bf q}_{1}$-type field observed at ${\bf x}$, in response to this source. The spatial Fourier transform of the downgoing component at $\partial\mathbb{D}_{F}$ (i.e., just below the source) is proportional to the upper-right submatrix of the decomposition operator of equation (118), according to

[TABLE]

(Wapenaar, 2022). To compensate for the effects of the inverse matrix $\{\tilde{\bf\Delta}_{1}({{\bf s}},x_{3,F})\}^{-1}$, we define a modified Green’s matrix as

[TABLE]

such that its downgoing component at $\partial\mathbb{D}_{F}$ is given by

[TABLE]

or, after applying an inverse spatial Fourier transformation

[TABLE]

The upgoing response at $\partial\mathbb{D}_{F}$ to this downgoing source field is by definition the reflection response, hence

[TABLE]

The field at ${\bf x}$ inside the medium consists of ${\bf\Gamma}_{12}({\bf x},{\bf x}_{S})$ and ${\bf\Gamma}_{22}({\bf x},{\bf x}_{S})$, where $\tilde{\bf\Gamma}_{22}({\bf x},{\bf s},x_{3,S})=\tilde{\bf G}_{22}({\bf x},{\bf s},x_{3,S})\tilde{\bf\Delta}_{1}({\bf s},x_{3,F})$, with $\tilde{\bf G}_{22}({\bf x},{\bf s},x_{3,S})$ being the Green’s function for the ${\bf q}_{2}$-type field observed at ${\bf x}$. Substituting ${\bf q}_{\alpha}({\bf x})={\bf\Gamma}_{\alpha 2}({\bf x},{\bf x}_{S})$ and ${\bf q}_{1}^{\pm}({\bf x}_{F})={\bf\Gamma}_{12}^{\pm}({\bf x}_{F},{\bf x}_{S})$ into equation (122), using equations (126) and (127), we obtain

[TABLE]

(no summation for repeated subscripts) for $x_{3}\geq x_{3,F}$. This is a generalisation of equations (54) and (55) and a starting point for developing a unified Marchenko method for full wave fields, accounting for evanescent waves inside the medium. Once the focusing functions are found, they can be used to retrieve the Green’s matrices ${\bf\Gamma}_{\alpha 2}({\bf x},{\bf x}_{S})$ for $\alpha=1,2$ (from equation (128)) and all components of the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ (from equations (120) and (121)).

Appendix B Unified transfer matrix and decomposed focusing functions

In this appendix, we extend the theory of section 3 to unified wave fields.

B.1 Unified transfer matrix

We introduce the unified transfer matrix as follows. Given the downgoing and upgoing fields ${\bf q}^{+}({\bf x}_{F})$ and ${\bf q}^{-}({\bf x}_{F})$ contained in vector ${\bf b}({\bf x}_{F})$ at the boundary $\partial\mathbb{D}_{F}$, we transfer these fields to depth level $x_{3}$ via

[TABLE]

for $x_{3}\geq x_{3,F}$. The unified transfer matrix {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) is partitioned as follows

[TABLE]

with {{{\mbox{\boldmath{\cal T}}}}}^{\pm,\pm} being $\frac{N}{2}\times\frac{N}{2}$ submatrices. Analogous to equation (58), matrix {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) is related to the unified propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$ of equation (88) via

[TABLE]

with {{{{\mbox{\boldmath{\cal D}}}}}}({\bf x}_{F}) and {{{\mbox{\boldmath{\cal D}}}}}^{-1}({\bf x}) being the inverse spatial Fourier transforms of $\tilde{\bf D}({\bf s},x_{3,F})$ and $\{\tilde{\bf D}({\bf s},x_{3})\}^{-1}$, defined in equations (94) and (118), respectively. Unlike in the acoustic situation, where {{{\mbox{\boldmath{\cal L}}}}}({\bf x}_{F}) and {{{\mbox{\boldmath{\cal L}}}}}^{-1}({\bf x}) in equation (58) account for lateral variations of the medium parameters, the unified matrices $\tilde{\bf D}({\bf s},x_{3,F})$ and $\{\tilde{\bf D}({\bf s},x_{3})\}^{-1}$ are defined for laterally invariant medium parameters at depths $x_{3,F}$ and $x_{3}$. For $\tilde{\bf D}({\bf s},x_{3,F})$ this is not a restriction, since $x_{3,F}$ is the depth of the boundary $\partial\mathbb{D}_{F}$ between the inhomogeneous medium and the homogeneous upper half-space. However, for $\{\tilde{\bf D}({\bf s},x_{3})\}^{-1}$ it implies that this operator can only be applied at depths where no lateral variations occur.

B.2 Relation between unified transfer matrix and decomposed Marchenko focusing functions

Assuming there are no lateral variations at a specific depth level $x_{3}$, we use the spatial Fourier transformation defined in equation (14) along the horizontal components of the first space variable to express the transfer matrix (analogous to equation (59)) as

[TABLE]

with $\{\tilde{\bf D}({\bf s},x_{3})\}^{-1}$ defined in equation (118) and $\tilde{\bf Y}({\bf s},x_{3},{\bf x}_{F})$ being the Fourier transform of ${\bf Y}({\bf x},{\bf x}_{F})$ defined in equation (116). Analogous to equation (60), we obtain for the right column of \tilde{{{\mbox{\boldmath{\cal T}}}}}({\bf s},x_{3},{\bf x}_{F})

[TABLE]

with $\tilde{\bf F}^{+}({\bf s},x_{3},{\bf x}_{F})$ and $\tilde{\bf F}^{-}({\bf s},x_{3},{\bf x}_{F})$ being the downgoing and upgoing parts at $x_{3}$ of $\tilde{\bf F}_{1}({\bf s},x_{3},{\bf x}_{F})$. For the left column of \tilde{{{\mbox{\boldmath{\cal T}}}}}({\bf s},x_{3},{\bf x}_{F}) we analyse the following expression

[TABLE]

Using equations (97) and (158) in equation (134) gives

[TABLE]

By comparing this with equation (133) we find that the expression in equation (135) is equal to

[TABLE]

Combining the right column (equation (133)) and the left column (equation (136)), we obtain the following expression for the unified transfer matrix

[TABLE]

or, in the space domain,

[TABLE]

This is the generalisation of equation (63).

B.3 Representations with decomposed unified Marchenko focusing functions

Substituting the expressions for ${\bf b}({\bf x})$ and {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) into equation (129) gives the following representations for the downgoing and upgoing fields, ${\bf q}_{1}^{+}({\bf x})$ and ${\bf q}_{1}^{-}({\bf x})$ respectively, inside the inhomogeneous medium

[TABLE]

for $x_{3}\geq x_{3,F}$. These expressions are exact and hold for dissipative media. Making similar substitutions as in section A.4 we obtain

[TABLE]

for $x_{3}\geq x_{3,F}$. Here ${\bf\Gamma}_{12}^{+}({\bf x},{\bf x}_{S})$ and ${\bf\Gamma}_{12}^{-}({\bf x},{\bf x}_{S})$ stand for the downgoing and upgoing part of the Green’s function ${\bf\Gamma}_{12}({\bf x},{\bf x}_{S})$. These equations are generalisations of equations (67) and (68) and form a starting point for developing a unified Marchenko method for decomposed wave fields. Once the focusing functions are found, they can be used to retrieve the decomposed Green’s functions ${\bf\Gamma}_{12}^{+}({\bf x},{\bf x}_{S})$ and ${\bf\Gamma}_{12}^{-}({\bf x},{\bf x}_{S})$ (from equations (141) and (142)) and all components of the transfer matrix {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}) (from equation (138)). Versions of the Marchenko method based on expressions similar to equations (141) and (142) have already been implemented for the retrieval of decomposed elastodynamic Green’s functions in lossless media, ignoring evanescent waves (Wapenaar & Slob, 2014; da Costa Filho et al., 2014; Reinicke & Wapenaar, 2019; Reinicke et al., 2020).

Appendix C Symmetry properties of the propagator and decomposition matrices

Let $N\times N$ propagator matrices ${\bf W}({\bf x},{\bf x}_{A})$ and ${\bf W}({\bf x},{\bf x}_{B})$ be two independent solutions of equation (8) (with unified operator matrix {{\mbox{\boldmath{\cal A}}}}({\bf x}) defined in equation (69)), with boundary condition (9), modified for coordinate vectors ${\bf x}_{A}$ and ${\bf x}_{B}$. We show that the quantity $\int_{{\mathbb{R}}^{2}}{\bf W}^{t}({\bf x}_{\rm H},x_{3},{\bf x}_{A}){\bf N}{\bf W}({\bf x}_{\rm H},x_{3},{\bf x}_{B}){\rm d}^{2}{\bf x}_{\rm H}$ is a ‘propagation invariant’, meaning that it is independent of $x_{3}$ (Haines, 1988; Kennett et al., 1990; Koketsu et al., 1991; Takenaka et al., 1993). To this end we take the derivative in the $x_{3}$-direction, apply the product rule for differentiation, use equation (8) and symmetry relation (70), according to

[TABLE]

This confirms that the integral is a propagation invariant. In a similar way, using symmetry relation (71), it can be shown that $\int_{{\mathbb{R}}^{2}}\bar{\bf W}^{\dagger}({\bf x}_{\rm H},x_{3},{\bf x}_{A}){\bf K}{\bf W}({\bf x}_{\rm H},x_{3},{\bf x}_{B}){\rm d}^{2}{\bf x}_{\rm H}$ is also a propagation invariant. Using boundary condition (9), modified for $x_{3}=x_{3,A}$ and $x_{3}=x_{3,B}$, we find from the first propagation invariant

[TABLE]

and from the second propagation invariant

[TABLE]

From equations (144) and (145), using ${\bf K}{\bf N}^{-1}={\bf J}$, we find

[TABLE]

This equation is used in section 2.3 and Appendix A.3 in the derivation of the relation between the propagator matrix and the Marchenko focusing functions.

To derive a symmetry property for $\tilde{\bf D}_{1}^{\pm}=\tilde{\bf L}_{2}^{\pm}(\tilde{\bf L}_{1}^{\pm})^{-1}$ (equation (96)), we start by Fourier transforming symmetry relations (70) – (72), assuming the medium is laterally invariant at the depth level where the transformation is applied. This gives

[TABLE]

The eigenvalue decomposition of matrix $\tilde{\bf A}$ is defined in equation (89), with the partitioning of $\tilde{\bf\Lambda}$ and $\tilde{\bf L}$ defined in equation (90). For all wave phenomena mentioned in Appendix A.1, the eigenvalue matrix \tilde{{{\mbox{\boldmath\Lambda}}}} obeys the following symmetry relations

[TABLE]

Given equations (89) and (147) – (151), matrix $\tilde{\bf L}$ can be normalized such that

[TABLE]

From the latter two equations, we obtain

[TABLE]

or, using the partitioning of $\tilde{\bf L}$ as defined in equation (90),

[TABLE]

with ${\bf J}_{11}$ and ${\bf J}_{22}$ being the upper-left and lower-right submatrices of matrix ${\bf J}$, and ${\bf K}_{12}$ being the upper-right (= lower-left) submatrix of matrix ${\bf K}$. From equations (155) and (156) it follows that $\tilde{\bf D}_{1}^{\pm}$ as defined in equation (96) obeys the following symmetry relation

[TABLE]

Finally, using $\tilde{\bf\Delta}_{1}=\tilde{\bf D}_{1}^{+}-\tilde{\bf D}_{1}^{-}$ (equation 119), we obtain

[TABLE]

Captions

Fig 1. Relations between the propagator matrix ${\bf W}({\bf x},{\bf x}_{F})$, the transfer matrix {{{\mbox{\boldmath{\cal T}}}}}({\bf x},{\bf x}_{F}), and the Marchenko focusing functions $F^{p}({\bf x},{\bf x}_{F})$ and $F^{v}({\bf x},{\bf x}_{F})$ (right column of ${\bf Y}({\bf x},{\bf x}_{F})$). The green and yellow double-sided arrows indicate full wave fields (implicitly consisting of downgoing and upgoing components), whereas the red and blue single-sided arrows indicate decomposed downgoing and upgoing wave fields, respectively.

Fig 2. (a) Horizontally layered medium. (b) Propagator matrix component $W^{p,p}(s_{1},x_{3},x_{3,F},\tau)$ (for fixed $s_{1}=1/3500$ m/s). (c) Propagator matrix component $W^{p,v}(s_{1},x_{3},x_{3,F},\tau)$.

Fig 3. (a) Focusing function $F^{p}(s_{1},x_{3},x_{3,F},\tau)$ (for fixed $s_{1}=1/3500$ m/s). (b) Time-reversed focusing function $F^{p}(s_{1},x_{3},x_{3,F},-\tau)$.

Fig 4. (a) Decomposed focusing function $F^{-}(s_{1},x_{3},x_{3,F},\tau)$ (for fixed $s_{1}=1/3500$ m/s). (b) Decomposed focusing function $F^{+}(s_{1},x_{3},x_{3,F},\tau)$ .