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

SUMMARY

Many seismic imaging methods use wavefield extrapolation operators to redatum sources and receivers from the surface into the subsurface. We discuss wavefield 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 wavefield 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 wavefield that focuses at a designated point in space at zero time. Marchenko focusing functions are useful for retrieving the wavefield 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 wavefield vector, a transfer matrix (according to the definition used in this paper) ‘transfers’ a decomposed wavefield 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 wavefield focusing, Green’s function retrieval and imaging.

1 INTRODUCTION

In many seismic imaging methods, wavefield 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 wavefield extrapolation operators, which only account for primaries. To account for internal multiple reflections in redatuming, more advanced wavefield extrapolation operators are required. This paper is not about the redatuming process itself, but about wavefield 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 wavefield 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 wavefield vector that the propagator matrix acts on contains components of the full wavefield (e.g. particle velocity and stress). Here ‘full’ means that the wavefield implicitly consists of downgoing and upgoing, propagating and evanescent waves.

A Marchenko focusing function is a wavefield 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 wavefield 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 wavefield 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 wavefield 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 wavefield 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 wavefield 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 × 2 propagator matrix for acoustic wavefields 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 × 2 transfer matrix for acoustic wavefields 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 × 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

(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, q is a vector containing the wavefield components p (acoustic pressure) and $v_3$ (vertical component of the particle velocity), both as a function of the space coordinate vector x = (x₁, x₂, x₃) (with positive x₃ denoting depth) and the angular frequency ω, hence,

Operator ∂₃ stands for the partial differential operator ∂/∂x₃. The space- and frequency-dependent operator matrix $\mathit{A}$ is defined as

where κ(x, ω) is the compressibility, ρ(x, ω) the mass density and i the imaginary unit. Operator ∂_α stands for the partial differential operator ∂/∂x_α. 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 κ and ρ may be frequency-dependent and complex-valued, with (for positive ω) ℑ(κ) ≥ 0 and ℑ(ρ) ≥ 0, where ℑ denotes the imaginary part. For later convenience we rewrite the operator matrix as follows:

Here ${\mathcal{H}}_2(\mathbf{x},\omega )$ is the Helmholtz operator, defined as

with wavenumber k(x, ω) defined via

(Brekhovskikh 1960; Wapenaar et al. 2001). Finally, vector d in eq. (1) contains source terms, according to

Here ${\hat{f}}_{\alpha }(\mathbf{x},\omega )$ and ${\hat{f}}_3(\mathbf{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}(\mathbf{x},\omega )$ is the volume injection-rate density (where $\hat{q}$ is to be distinguished from the wavefield vector 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₃ = x_{3, F}. We define a coordinate vector x_F at this boundary as x_F = (x_{1, F}, x_{2, F}, x_{3, F}) (with fixed x_{3, F}). We introduce the propagator matrix W(x, x_F) as a solution of wave eq. (1) for the source-free situation, according to

with boundary condition

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

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

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 wavefield component (p or $v_3$⁠) it acts on at x_F, whereas the first superscript refers to the wavefield component it contributes to at x. Eq. (10) is illustrated in the upper-left frame of Fig. 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₃ indicates an arbitrary depth level inside the inhomogeneous medium.

By applying eq. (10) recursively, it follows that W obeys the following recursive expression

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

from which it follows that W(x_F, x) is the inverse of W(x, x_F).

The propagator matrix W(x, x_F) accounts for primaries and multiples between x_{3, F} and x₃ 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 ${\mathcal{H}}_2$ (rather than on its square-root, appearing in one-way wave equations), Kosloff & Baysal (1983) argue that the numerical evaluation of eq. (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 Fig. 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(x, ω) along the horizontal coordinate x_H for constant x₃ as

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

Next, we define the inverse temporal Fourier transformation for constant s and x₃ as

where ℜ denotes the real part and τ is the intercept time (Stoffa 1989). Applying this transformation to $\tilde{\mathbf{W}}(\mathbf{s},x_3,x_{3,F})$ we obtain W(s, x₃, x_{3, F}, τ), with boundary condition W(s, x_{3, F}, x_{3, F}, τ) = Iδ(τ) and with W(s, x₃, x_{3, F}, τ) obeying the recursive expression

where the inline asterisk denotes temporal convolution. Although the numerical modelling is most efficiently done in the slowness-frequency domain (using eq. 15), the results are more conveniently interpreted when displayed in the slowness intercept-time domain. Setting s₂ = 0, the components $W^{p,p}$(s₁, x₃, x_{3, F}, τ) and $W^{p,v}$(s₁, x₃, x_{3, F}, τ), with boundary conditions $W^{p,p}$(s₁, x_{3, F}, x_{3, F}, τ) = δ(τ) and $W^{p,v}$(s₁, x_{3, F}, x_{3, F}, τ) = 0, are shown in Figs 2(b) and (c) for fixed s₁ = 1/3500 s m⁻¹, as a function of intercept time τ and depth x₃. 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₃ = x_{3, F} = 0 m represent the aforementioned boundary conditions. Note that $W^{p,p}$ and $W^{p,v}$ are, for each depth x₃, even and odd functions, respectively, of intercept time τ. The recursive character, described by eq. (17), is manifest in Figs 2(b) and (c). The propagation velocity in the layer between x₃ = 760 m and x₃ = 800 m equals 3600 m s⁻¹, which implies that for the chosen horizontal slowness s₁ = 1/3500 s m⁻¹ we have sin θ = cs₁ > 1 (i.e. θ is complex-valued). Hence, waves become ‘evanescent’ in this layer. The wavefield 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 wavefield 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 ρ₀ and propagation velocity c₀. 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 $\mathit{A}$ (in the space–frequency domain) as follows:

with

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

(Wapenaar et al. 2001), where g(x_H) and h(x_H) are test functions in the horizontal plane with ‘sufficient decay at infinity’. Operator ${\mathcal{H}}_1$⁠, as defined in eq. (19) is not symmetric, but operator $\frac{1}{\rho }{\mathcal{H}}_1$ and its inverse, both appearing in eq. (20), are symmetric. We use the operator matrix $\mathit{L}$ to express the wavefield vector q(x) in terms of downgoing and upgoing waves p⁺(x) and p⁻(x) via

with

Note that these equations imply

hence, the downgoing and upgoing waves p⁺ and p⁻ are pressure-normalized. An advantage of pressure-normalized (or, more generally, field-normalized) decomposition is that the decomposed quantities simply add up to a field quantity [acoustic pressure in the case of eq. (24)]. This property does not apply to flux-normalized decomposed wavefields (Frasier 1970; Kennett et al. 1978; Ursin 1983). On the other hand, an advantage of flux-normalized decomposition is that the underlying equations obey more simple symmetry properties. For a comprehensive discussion on field-normalized versus flux-normalized decomposition in inhomogeneous media, see de Hoop (1996) and Wapenaar (2020). In this paper, we use field-normalized 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 eq. (22) at $\partial {\mathbb{D}}_F$ to derive focusing functions and express them in the components of the propagator matrix and vice versa. Substituting eq. (22), with x replaced by x_F, into the right-hand side of eq. (10) gives

for x₃ ≥ x_{3, F}, with

or

The operators $\pm \frac{1}{\omega {\rho }_0}{\mathcal{H}}_1({\mathbf{x}}_F)$ in eq. (27) act, via eq. (25), on p^±(x_F). However, since these operators are symmetric [in the sense of eq. (21)], we may replace the actions of these operators on p^±(x_F) by actions on the elements $W^{p,v}$(x, x_F) and $W^{v,v}$(x, x_F). To be more specific, if we partition Y(x, x_F) as follows:

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

We analyse these expressions one by one. First, we consider the element $Y^{p,-}$. From eqs (25) and (28) it can be seen that the superscript p refers to the acoustic pressure p(x) contained in q(x) and superscript − refers to the upgoing wavefield component p⁻(x_F) in p(x_F). Using eqs (9), (11) and (29) we obtain

which is a focusing condition. Hence, we define

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

which is also a focusing condition, but somewhat more complicated than eq. (31) because of the mix of the involved wavefield components $v_3$(x) and p⁻(x_F). Hence, we define

with $F^v$(x, x_F) denoting the particle velocity counterpart of the focusing function $F^p$(x, x_F) (note that the definition of $F^v$(x, x_F) is different from that in Wapenaar (2022), to facilitate the derivations below). The focusing functions $F^p$(x, x_F) and $F^v$(x, x_F), which together form the right-hand column of matrix Y(x, x_F), are illustrated in the lower frame of Fig. 1. They resemble the focusing function f₂ 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₁, which focuses inside the medium). However, there are also some notable differences. First, f₂(x, 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^{+}(\mathbf{x},{\mathbf{x}}_F)$ and $f_2^{-}(\mathbf{x},{\mathbf{x}}_F)$ respectively, at 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₃. In contrast, $F^p$(x, x_F) and $F^v$(x, x_F) are defined in the actual (i.e. untruncated) medium and represent the full pressure and vertical particle velocity at x of a field that focuses at 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 eqs (32) and (34), accounts for evanescent waves. Last but not least, $F^p$ and $F^v$ hold for dissipative media and they are normalized differently from f₂.

Before we analyse the elements in the left-hand column of matrix Y(x, x_F), we introduce an adjoint medium, with parameters $\overline{\kappa }(\mathbf{x})={\kappa }^{\ast }(\mathbf{x})$ and $\overline{\rho }(\mathbf{x})={\rho }^{\ast }(\mathbf{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 ω) $\mathrm{\Im }(\overline{\kappa })\le 0$ and $\mathrm{\Im }(\overline{\rho })\le 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 $\overline{\mathit{A}}$ and the Helmholtz operator ${\overline{\mathcal{H}}}_2$ of the adjoint medium are defined similarly as $\mathit{A}$ and ${\mathcal{H}}_2$ in eqs (4) and (5), respectively, but with κ(x) and ρ(x) replaced by $\overline{\kappa }(\mathbf{x})$ and $\overline{\rho }(\mathbf{x})$⁠, respectively. Hence, ${\overline{\mathcal{H}}}_2={\mathcal{H}}_2^{\ast }$⁠. Analogous to eqs (8) and (9), we define the propagator matrix $\overline{\mathbf{W}}(\mathbf{x},{\mathbf{x}}_F)$ of the adjoint medium as the solution of ${\partial }_3\overline{\mathbf{W}}(\mathbf{x},{\mathbf{x}}_F)=\overline{\mathit{A}}(\mathbf{x})\overline{\mathbf{W}}(\mathbf{x},{\mathbf{x}}_F)$⁠, with boundary condition $\overline{\mathbf{W}}(\mathbf{x},{\mathbf{x}}_F){|}_{x_3=x_{3,F}}=\mathbf{I}\delta ({\mathbf{x}}_{\mathrm{H}}-{\mathbf{x}}_{\mathrm{H},F})$⁠. In Appendix  C we derive

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

(Wapenaar et al. 2001). Using eqs (35) and (36) in eqs (29) and (30), we find ${\overline{Y}}^{p,+}(\mathbf{x},{\mathbf{x}}_F)=Y^{p,-\ast }(\mathbf{x},{\mathbf{x}}_F)$ and ${\overline{Y}}^{v,+}(\mathbf{x},{\mathbf{x}}_F)=-Y^{v,-\ast }(\mathbf{x},{\mathbf{x}}_F)$⁠. Hence, using eqs (32) and (34), we find for the elements in the left-hand column of matrix Y(x, x_F)

For matrix Y(x, x_F) we thus obtain

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

with ${\mathit{L}}^{-1}$ defined in eq. (20). Since operator $\frac{1}{\rho }{\mathcal{H}}_1$ is symmetric, its inverse ${\mathcal{H}}_1^{-1}\rho$ is symmetric as well. Hence, in eq. (40) these operators can be taken to act on the elements of matrix Y(x, x_F). This yields

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 eqs (41)–(44) simplify to

We illustrate the focusing function and its relation with the propagator matrix with a numerical example. Applying the transformations of eqs (14) and (16) to eq. (32) (assuming a laterally invariant medium), taking x_F = (0, 0, x_{3, F}) and s₂ = 0, we obtain

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

For the acoustic case all components are symmetric in s₁, that is $F^p$( − s₁, x₃, x_{3, F}, −τ) = $F^p$(s₁, x₃, x_{3, F}, −τ), etc. Hence, eqs (50) and (51) show how the even and odd components $W^{p,p}$(s₁, x₃, x_{3, F}, τ) and $W^{p,v}$(s₁, x₃, x_{3, F}, τ) of Figs 2(b) and (c) are obtained from the focusing function $F^p$(s₁, x₃, x_{3, F}, τ) and its time-reversal $F^p$(s₁, x₃, x_{3, F}, −τ) of Fig. 3.

2.4 Representations with acoustic Marchenko focusing functions

Substituting the expressions for q(x), p(x_F) and Y(x, x_F) (eqs 2, 23 and 39) into eq. (25), gives the following representations for the acoustic pressure p(x) and the vertical particle velocity $v_3$(x) inside the inhomogeneous medium

for x₃ ≥ x_{3, F}. These expressions are exact and hold for dissipative media. Eq. (52) is a generalisation of eq. (17) of Wapenaar & de Ridder (2022) for dissipative media.

We use eqs (52) and (53) to derive representations for Green’s functions between the boundary $\partial {\mathbb{D}}_F$ and any position x inside the medium. To this end, we define a unit point source of vertical force at 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^{+}({\mathbf{x}}_F)=\frac{1}{2}\delta ({\mathbf{x}}_{\mathrm{H},F}-{\mathbf{x}}_{\mathrm{H},S})$⁠, where x_{H, S} denotes the horizontal coordinates of x_S. The upgoing field at $\partial {\mathbb{D}}_F$ is the reflection response to this downgoing source field, hence $p^{-}({\mathbf{x}}_F)=\frac{1}{2}R({\mathbf{x}}_F,{\mathbf{x}}_S)$⁠. The field at x inside the medium is the Green’s response to the source at x_S, hence p(x) = $G^{p,f}$(x, x_S) and $v_3$(x) = $G^{v,f}$(x, x_S). Here the second superscript (f) refers to the vertical force source at x_S, whereas the first superscripts (p and v) refer to the observed quantities (pressure and vertical particle velocity) at x. Substitution of these expressions for p^±(x_F), p(x) and $v_3$(x) into eqs (52) and (53) gives

for x₃ ≥ x_{3, F}. Slob (2016) derived similar representations for decomposed wavefields 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(x_F, 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}$(x, x_S) and $G^{v,f}$(x, x_S) [from eqs (54) and (55)] and all components of the propagator matrix W(x, x_F) [from eqs (45) to (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⁺(x_F) and p⁻(x_F) at the boundary $\partial {\mathbb{D}}_F$⁠, we ‘transfer’ these fields to downgoing and upgoing fields p⁺(x) and p⁻(x) at any depth level x₃ inside the medium using the following expression:

for x₃ ≥ x_{3, F}. Vectors p(x_F) and p(x) contain the downgoing and upgoing fields at depth levels x_{3, F} and x₃ (eq. 23). We call $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ the transfer matrix, which we partition as follows:

For each component of this matrix, the superscripts refer to the propagation direction at x and at x_F, respectively. Eq. (56) is illustrated in the upper-right frame of Fig. 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 $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ for laterally varying media. Substituting eq. (22) into eq. (10) we obtain eq. (56), with

with $\mathit{L}(\mathbf{x})$ and its inverse defined in eq. (20). Eq. (58), which relates the transfer matrix $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ to the propagator matrix W(x, x_F), is illustrated in the upper half of Fig. 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 eqs (26) and (58) we find

According to eq. (39), the right-hand column of Y(x, x_F) contains $F^p$(x, x_F) and $F^v$(x, x_F), that is, the pressure and vertical particle velocity components at x of the focusing function. Hence, analogous to $\mathbf{p}(\mathbf{x})={\mathit{L}}^{-1}(\mathbf{x})\mathbf{q}(\mathbf{x})$⁠, we obtain for the right-hand column of $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$

with F⁺(x, x_F) and F⁻(x, x_F) being the downgoing and upgoing parts at x of the focusing function $F^p$(x, x_F).

According to eq. (39), the left-hand column of Y(x, x_F) contains ${\overline{F}}^{p\ast }(\mathbf{x},{\mathbf{x}}_F)$ and $-{\overline{F}}^{v\ast }(\mathbf{x},{\mathbf{x}}_F)$⁠. Hence, for the left-hand column of $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ we obtain

or, using ${\mathcal{H}}_1={\overline{\mathcal{H}}}_1^{\ast }$ (eq. 36) and $\rho ={\overline{\rho }}^{\ast }$⁠,

Comparing this with eq. (60) we find that this gives a vector with ${\overline{F}}^{-\ast }(\mathbf{x},{\mathbf{x}}_F)$ and ${\overline{F}}^{+\ast }(\mathbf{x},{\mathbf{x}}_F)$⁠. This is the left-hand column of $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$⁠. Hence, we have obtained

see the upper-right frame of Fig. 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 eqs (14) and (16) to eq. (63), taking x_F = (0, 0, x_{3, F}), we obtain

Dukalski et al. (2022a, b) used a recursive approach and obtained an expression similar to eq. (64). In their derivation they used a path-reversal operator $\mathcal{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, $\mathcal{P}\{F^{\pm }(\mathbf{s},x_3,x_{3,F},\tau )\}$ is equivalent with ${\overline{F}}^{\pm }(-\mathbf{s},x_3,x_{3,F},-\tau )$⁠.

For the lossless medium of Fig. 2(a), the decomposed focusing functions F⁻(s₁, x₃, x_{3, F}, τ) and F⁺(s₁, x₃, x_{3, F}, τ) for s₁ = 1/3500 m s⁻¹ and s₂ = 0 are shown in Figs 4(a) and (b), respectively. For each x₃, the function F⁻(s₁, x₃, x_{3, F}, τ) can be seen as the intricate field that needs to be emitted upward from x₃ to arrive as a single upward propagating pulse at the focal depth x_{3, F} at τ = 0. For the same x₃, the function F⁺(s₁, x₃, x_{3, F}, τ) is the downward reflected response to F⁻(s₁, x₃, x_{3, F}, τ). Figs 4(b) and (a) together form the right-hand column of the transformed transfer matrix $\mathit{T}(s_1,x_3,x_{3,F},\tau )$⁠. Their superposition gives the focusing function $F^p$(s₁, x₃, x_{3, F}, τ), shown in Fig. 3(a).

3.3 Representations with decomposed acoustic Marchenko focusing functions

Substituting the expressions for p(x) and $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ (eqs 23 and 63) into eq. (56), gives the following representations for the downgoing and upgoing components of the acoustic pressure, p⁺(x) and p⁻(x), respectively, inside the inhomogeneous medium

for x₃ ≥ x_{3, F}. These expressions are exact and hold for dissipative media. Making similar substitutions as in Section 2.4 we obtain

for x₃ ≥ x_{3, F}. Here G^{±, f}(x, x_S) stands for the downgoing (+) and upgoing (−) part of the Green’s function $G^{p,f}$(x, x_S).

When the medium is lossless and when evanescent waves are neglected at $\partial {\mathbb{D}}_F$ and at depth level x₃ 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(x_F, 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}(x, x_S) and G^{−, f}(x, x_S) [from eqs (67) and (68)] and all components of the transfer matrix $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ [from eq. (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 wavefield 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.

References

APPENDIX A: UNIFIED PROPAGATOR MATRIX AND FOCUSING FUNCTIONS

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

A1 Unified matrix-vector wave equation

Consider matrix-vector wave eq. (1). We partition the N × 1 wave field vector q, the N × 1 source vector d and the N × N operator matrix $\mathit{A}$ as follows

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 $\mathit{A}$ obeys the following symmetry properties

with

where O and 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 † 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, eqs (A1) – (A4) also hold for the matrix-vector form of the Dirac equation for the 4 × 1 spinor $\mathit{\psi }$⁠. The Dirac equation is given by (Sakurai 1967)

(summation over μ from 1 to 4), with the Dirac spinor partitioned as

and with

with V(x) the space-dependent potential, h Planck’s constant, c the speed of light, e the electron charge and ${\mathit{\sigma }}_1$⁠, ${\mathit{\sigma }}_2$ and ${\mathit{\sigma }}_3$ the Pauli matrices, defined as

Assuming a time-dependence exp ( − iEt/ℏ), we replace ∂₄ by −(E − eV)/ℏc. Eq. (A6) can be rewritten as follows

or

Pre-multiplying all terms by $-i{\mathit{\sigma }}_3$ and bringing the ∂₃-operators to the left-hand side, using ${\mathit{\sigma }}_3{\mathit{\sigma }}_1=\mathrm{i}{\mathit{\sigma }}_2$⁠, ${\mathit{\sigma }}_3{\mathit{\sigma }}_2=-\mathrm{i}{\mathit{\sigma }}_1$ and ${\mathit{\sigma }}_3{\mathit{\sigma }}_3=\mathbf{I}$⁠, yields a set of equations which can be recast in the form of eq. (1), with 4 × 1 vectors q and d and 4 × 4 operator matrix $\mathit{A}$ partitioned as in eq. (A1), with

With these definitions of the operator submatrices, matrix $\mathit{A}$ obeys symmetry relations (A2) – (A4), with N, K and J defined as

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

A2 Unified propagator matrix

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

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

A3 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{\mathbf{A}}$ (the spatial Fourier transform of operator matrix $\mathit{A}$⁠), as follows

with

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

with

Note that these equations imply ${\tilde{\mathbf{q}}}_1={\tilde{\mathbf{L}}}_1^{+}{\tilde{\mathbf{p}}}^{+}+{\tilde{\mathbf{L}}}_1^{-}{\tilde{\mathbf{p}}}^{-}$⁠. Similar as in Section 2.3, eq. (24), we continue with downgoing and upgoing waves ${\tilde{\mathbf{q}}}_1^{+}$ and ${\tilde{\mathbf{q}}}_1^{-}$ which are field-normalized such that ${\tilde{\mathbf{q}}}_1={\tilde{\mathbf{q}}}_1^{+}+{\tilde{\mathbf{q}}}_1^{-}$⁠. To this end, we define ${\tilde{\mathbf{q}}}_1^{\pm }={\tilde{\mathbf{L}}}_1^{\pm }{\tilde{\mathbf{p}}}^{\pm }$ and we replace eq. (A23) by

where

with

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

with J₁₁ and J₂₂ being the $\frac{N}{2}\times \frac{N}{2}$ submatrices of N × N matrix J. From eq. (A5) we have for all wave phenomena except for the Dirac equation J₁₁ = −J₂₂ = I, and from eq. (A19) we have for the Dirac equation ${\mathbf{J}}_{11}={\mathbf{J}}_{22}={\mathit{\sigma }}_2$⁠.

We use eq. (A25) at $\partial {\mathbb{D}}_F$ and the properties of matrix ${\widetilde{\overline{\mathbf{D}}}}_1^{\pm }(\mathbf{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 eq. (A25) for x₃ = x_{3, F} into a transformed version of eq. (10). This equation contains the propagator matrix W(x, x_F). For a function of two space variables, u(x, x_F) (with x_F at $\partial {\mathbb{D}}_F$⁠), we define the spatial Fourier transformation along the horizontal components of the second space variable as

and its inverse as

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

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

Substitution of eq. (A25) into eq. (A32) gives

for x₃ ≥ x_{3, F}, with

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

Using equation (A26) and the spatial Fourier transform of eq. (A20), we obtain

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

with ${\tilde{\mathbf{F}}}_1(\mathbf{x},\mathbf{s},x_{3,F})$ denoting the spatial Fourier transform of the focusing function F₁(x, x_F) for wavefield component q₁, which focuses as an upgoing field at x = 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{\mathbf{Y}}}_2^{-}(\mathbf{x},\mathbf{s},x_{3,F})$⁠. Subscript 2 refers to wavefield component q₂ at x and superscript − refers again to the upgoing wavefield component ${\tilde{\mathbf{q}}}_1^{-}$ at x_{3, F}. For x₃ = x_{3, F} we obtain, using eq. (A33), ${\tilde{\mathbf{Y}}}_2^{-}(\mathbf{x},\mathbf{s},x_{3,F}){|}_{x_3=x_{3,F}}={\tilde{\mathbf{D}}}_1^{-}(\mathbf{s},x_{3,F})\exp \{i\omega \mathbf{s}\cdot {\mathbf{x}}_{\mathrm{H}}\}$⁠, which is a focusing condition, but somewhat more complicated than for ${\tilde{\mathbf{Y}}}_1^{-}(\mathbf{x},\mathbf{s},x_{3,F})$ because of the mix of wavefield components q₂ and ${\tilde{\mathbf{q}}}_1^{-}$⁠. Hence, we define

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

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

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

(no summation for repeated subscripts). Substituting eqs (A29) and (A42) into eqs (A37) and (A38) yields

Hence, using eqs (A39) and (A40), we find for the submatrices in the left-hand column of $\tilde{\mathbf{Y}}(\mathbf{x},\mathbf{s},x_{3,F})$

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

or, using the inverse Fourier transformation defined in eq. (A31),

This is a generalisation of eq. (39). Note that ${\tilde{\mathbf{F}}}_1$⁠, ${\tilde{\mathbf{F}}}_2$⁠, ${\widetilde{\overline{\mathbf{F}}}}_1^{\ast }$ and ${\widetilde{\overline{\mathbf{F}}}}_2^{\ast }$ are expressed in terms of the submatrices of the propagator matrix $\tilde{\mathbf{W}}(\mathbf{x},\mathbf{s},x_{3,F})$ via eqs (A37) – (A40), (A45) and (A46). Conversely, we can express the submatrices of the propagator matrix $\tilde{\mathbf{W}}(\mathbf{x},\mathbf{s},x_{3,F})$ in terms of the focusing functions ${\tilde{\mathbf{F}}}_1$⁠, ${\tilde{\mathbf{F}}}_2$⁠, ${\widetilde{\overline{\mathbf{F}}}}_1^{\ast }$ and ${\widetilde{\overline{\mathbf{F}}}}_2^{\ast }$⁠. To this end, we start with inverting eq. (A35), according to

with

Using eqs (A47) and (A50), we obtain

(no summation for repeated subscripts). These expressions are a generalisation of eqs (41) – (44). Those equations follow as a special case from eqs (A52) and (A53) by substituting J₁₁ = −J₂₂ = 1, ${\tilde{\mathbf{D}}}_1^{\pm }(\mathbf{s},x_{3,F})=\pm s_{3,0}/{\rho }_0$⁠, $\{{\tilde{\mathbf{\Delta }}}_1(\mathbf{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 }{\mathcal{H}}_1({\mathbf{x}}_F)$⁠.

A4 Representations with unified Marchenko focusing functions

Applying Parseval’s theorem to eq. (A34) and substituting the expressions for q(x), b(x_F) and Y(x, x_F) (eqs (A1), (A27) and (A48)), gives the following representation for the quantities q₁(x) and q₂(x) inside the inhomogeneous medium

(no summation for repeated subscripts) for x₃ ≥ x_{3, F}. This is a generalisation of eqs (52) and (53).

We use eq. (A54) to derive representations for Green’s functions between the boundary $\partial {\mathbb{D}}_F$ and any position x inside the medium. We define a unit d₂-type source (see eq. (A1)) at x_S just above $\partial {\mathbb{D}}_F$⁠. The $\frac{N}{2}\times \frac{N}{2}$ Green’s matrix G₁₂(x, x_S) stands for the q₁-type field observed at 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 eq. (A50), according to

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

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

or, after applying an inverse spatial Fourier transformation

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

The field at x inside the medium consists of Γ₁₂(x, x_S) and Γ₂₂(x, x_S), where ${\tilde{\mathbf{\Gamma }}}_{22}(\mathbf{x},\mathbf{s},x_{3,S})={\tilde{\mathbf{G}}}_{22}(\mathbf{x},\mathbf{s},x_{3,S}){\tilde{\mathbf{\Delta }}}_1(\mathbf{s},x_{3,F})$⁠, with ${\tilde{\mathbf{G}}}_{22}(\mathbf{x},\mathbf{s},x_{3,S})$ being the Green’s function for the q₂-type field observed at x. Substituting q_α(x) = Γ_α2(x, x_S) and ${\mathbf{q}}_1^{\pm }({\mathbf{x}}_F)={\mathbf{\Gamma }}_{12}^{\pm }({\mathbf{x}}_F,{\mathbf{x}}_S)$ into eq. (A54), using eqs (A58) and (A59), we obtain

(no summation for repeated subscripts) for x₃ ≥ x_{3, F}. This is a generalisation of eqs (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 Γ_α2(x, x_S) for α = 1, 2 (from eq. (A60)) and all components of the propagator matrix W(x, x_F) (from eqs (A52) and (A53)).

APPENDIX B: UNIFIED TRANSFER MATRIX AND DECOMPOSED FOCUSING FUNCTIONS

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

B1 Unified transfer matrix

We introduce the unified transfer matrix as follows. Given the downgoing and upgoing fields q⁺(x_F) and q⁻(x_F) contained in vector b(x_F) at the boundary $\partial {\mathbb{D}}_F$⁠, we transfer these fields to depth level x₃ via

for x₃ ≥ x_{3, F}. The unified transfer matrix $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ is partitioned as follows

with ${\mathit{T}}^{\pm ,\pm }$ being $\frac{N}{2}\times \frac{N}{2}$ submatrices. Analogous to eq. (58), matrix $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ is related to the unified propagator matrix W(x, x_F) of eq. (A20) via

with $\mathit{D}({\mathbf{x}}_F)$ and ${\mathit{D}}^{-1}(\mathbf{x})$ being the inverse spatial Fourier transforms of $\tilde{\mathbf{D}}(\mathbf{s},x_{3,F})$ and $\{\tilde{\mathbf{D}}(\mathbf{s},x_3){\}}^{-1}$⁠, defined in eqs (A26) and (A50), respectively. Unlike in the acoustic situation, where $\mathit{L}({\mathbf{x}}_F)$ and ${\mathit{L}}^{-1}(\mathbf{x})$ in eq. (58) account for lateral variations of the medium parameters, the unified matrices $\tilde{\mathbf{D}}(\mathbf{s},x_{3,F})$ and $\{\tilde{\mathbf{D}}(\mathbf{s},x_3){\}}^{-1}$ are defined for laterally invariant medium parameters at depths x_{3, F} and x₃. For $\tilde{\mathbf{D}}(\mathbf{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{\mathbf{D}}(\mathbf{s},x_3){\}}^{-1}$ it implies that this operator can only be applied at depths where no lateral variations occur.

B2 Relation between unified transfer matrix and decomposed Marchenko focusing functions

Assuming there are no lateral variations at a specific depth level x₃, we use the spatial Fourier transformation defined in eq. (14) along the horizontal components of the first space variable to express the transfer matrix (analogous to eq. (59)) as

with $\{\tilde{\mathbf{D}}(\mathbf{s},x_3){\}}^{-1}$ defined in eq. (A50) and $\tilde{\mathbf{Y}}(\mathbf{s},x_3,{\mathbf{x}}_F)$ being the Fourier transform of Y(x, x_F) defined in eq. (A48). Analogous to eq. (60), we obtain for the right-hand column of $\tilde{\mathit{T}}(\mathbf{s},x_3,{\mathbf{x}}_F)$

with ${\tilde{\mathbf{F}}}^{+}(\mathbf{s},x_3,{\mathbf{x}}_F)$ and ${\tilde{\mathbf{F}}}^{-}(\mathbf{s},x_3,{\mathbf{x}}_F)$ being the downgoing and upgoing parts at x₃ of ${\tilde{\mathbf{F}}}_1(\mathbf{s},x_3,{\mathbf{x}}_F)$⁠. For the left-hand column of $\tilde{\mathit{T}}(\mathbf{s},x_3,{\mathbf{x}}_F)$ we analyse the following expression

Using eqs (A29) and (C16) in eq. (B6) gives

By comparing this with eq. (B5) we find that the expression in eq. (B7) is equal to

Combining the right-hand column (eq. B5) and the left-hand column (eq. B8), we obtain the following expression for the unified transfer matrix

or, in the space domain,

This is the generalisation of eq. (63).

B3 Representations with decomposed unified Marchenko focusing functions

Substituting the expressions for b(x) and $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ into eq. (B1) gives the following representations for the downgoing and upgoing fields, ${\mathbf{q}}_1^{+}(\mathbf{x})$ and ${\mathbf{q}}_1^{-}(\mathbf{x})$ respectively, inside the inhomogeneous medium

for x₃ ≥ x_{3, F}. These expressions are exact and hold for dissipative media. Making similar substitutions as in Section A4 we obtain

for x₃ ≥ x_{3, F}. Here ${\mathbf{\Gamma }}_{12}^{+}(\mathbf{x},{\mathbf{x}}_S)$ and ${\mathbf{\Gamma }}_{12}^{-}(\mathbf{x},{\mathbf{x}}_S)$ stand for the downgoing and upgoing part of the Green’s function Γ₁₂(x, x_S). These equations are generalisations of eqs (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 ${\mathbf{\Gamma }}_{12}^{+}(\mathbf{x},{\mathbf{x}}_S)$ and ${\mathbf{\Gamma }}_{12}^{-}(\mathbf{x},{\mathbf{x}}_S)$ (from eqs (B13) and (B14)) and all components of the transfer matrix $\mathit{T}(\mathbf{x},{\mathbf{x}}_F)$ (from eq. (B10)). Versions of the Marchenko method based on expressions similar to eqs (B13) and (B14) 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 × N propagator matrices W(x, x_A) and W(x, x_B) be two independent solutions of eq. (8) (with unified operator matrix $\mathit{A}(\mathbf{x})$ defined in eq. (A1)), with boundary condition (9), modified for coordinate vectors x_A and x_B. We show that the quantity ${\int }_{{\mathbb{R}}^2}{\mathbf{W}}^t({\mathbf{x}}_{\mathrm{H}},x_3,{\mathbf{x}}_A)\mathbf{N}\mathbf{W}({\mathbf{x}}_{\mathrm{H}},x_3,{\mathbf{x}}_B){\mathrm{d}}^2{\mathbf{x}}_{\mathrm{H}}$ is a ‘propagation invariant’, meaning that it is independent of x₃ (Haines 1988; Kennett et al. 1990; Koketsu et al. 1991; Takenaka et al. 1993). To this end we take the derivative in the x₃-direction, apply the product rule for differentiation, use eq. (8) and symmetry relation (A2), according to

This confirms that the integral is a propagation invariant. In a similar way, using symmetry relation (A3), it can be shown that ${\int }_{{\mathbb{R}}^2}{\overline{\mathbf{W}}}^{†}({\mathbf{x}}_{\mathrm{H}},x_3,{\mathbf{x}}_A)\mathbf{K}\mathbf{W}({\mathbf{x}}_{\mathrm{H}},x_3,{\mathbf{x}}_B){\mathrm{d}}^2{\mathbf{x}}_{\mathrm{H}}$ is also a propagation invariant. Using boundary condition (9), modified for x₃ = x_{3, A} and x₃ = x_{3, B}, we find from the first propagation invariant

and from the second propagation invariant

From eqs (C2) and (C3), using KN⁻¹ = J, we find

This equation is used in Section 2.3 and Appendix A3 in the derivation of the relation between the propagator matrix and the Marchenko focusing functions.

To derive a symmetry property for ${\tilde{\mathbf{D}}}_1^{\pm }={\tilde{\mathbf{L}}}_2^{\pm }({\tilde{\mathbf{L}}}_1^{\pm }{)}^{-1}$ (eq. (A28)), we start by Fourier transforming symmetry relations (A2) – (A4), assuming the medium is laterally invariant at the depth level where the transformation is applied. This gives

The eigenvalue decomposition of matrix $\tilde{\mathbf{A}}$ is defined in eq. (A21), with the partitioning of $\tilde{\mathbf{\Lambda }}$ and $\tilde{\mathbf{L}}$ defined in eq. (A22). For all wave phenomena mentioned in Appendix A1, the eigenvalue matrix $\tilde{\mathbf{\Lambda }}$ obeys the following symmetry relations

Given eqs (A21) and (C5) – (C9), matrix $\tilde{\mathbf{L}}$ can be normalized such that

From the latter two equations, we obtain

or, using the partitioning of $\tilde{\mathbf{L}}$ as defined in eq. (A22),

with J₁₁ and J₂₂ being the upper-left and lower-right submatrices of matrix J, and K₁₂ being the upper-right (= lower-left) submatrix of matrix K. From eqs (C13) and (C14) it follows that ${\tilde{\mathbf{D}}}_1^{\pm }$ as defined in eq. (A28) obeys the following symmetry relation

Finally, using ${\tilde{\mathbf{\Delta }}}_1={\tilde{\mathbf{D}}}_1^{+}-{\tilde{\mathbf{D}}}_1^{-}$ (equation A51), we obtain