# Reciprocity and Representation Theorems for Flux- and Field-Normalised Decomposed Wave Fields

Authors: Kees Wapenaar
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

We consider wave propagation problems in which there is a preferred direction of propagation. To account for propagation in preferred directions, the wave equation is decomposed into a set of coupled equations for waves that propagate in opposite directions along the preferred axis. This decomposition is not unique. We discuss flux-normalised and field-normalised decomposition in a systematic way, analyse the symmetry properties of the decomposition operators, and use these symmetry properties to derive reciprocity theorems for the decomposed wave fields, for both types of normalisation. Based on the field-normalised reciprocity theorems, we derive representation theorems for decomposed wave fields. In particular, we derive double- and single-sided Kirchhoff-Helmholtz integrals for forward and backward propagation of decomposed wave fields. The single-sided Kirchhoff-Helmholtz integrals for backward propagation of field-normalised decomposed wave fields find applications in reflection imaging, accounting for multiple scattering.


---

## Body

## 1. Introduction

In many wave propagation problems, it is possible to define a preferred direction of propagation. For example, in ocean acoustics, waves propagate primarily in the horizontal direction in an acoustic wave guide, bounded by the water surface and the ocean bottom. Similarly, in communication engineering, microwaves or optical waves propagate as beams through electromagnetic or optical wave guides. Wave propagation in preferred directions is not restricted to wave guides. For example, in geophysical reflection imaging applications, seismic or electromagnetic waves propagate mainly in the vertical direction (downward and upward) through a laterally unbounded medium.To account for propagation in preferred directions, the wave equation for the full wave field can be decomposed into a set of coupled equations for waves that propagate in opposite directions along the preferred axis (for example, leftward and rightward in ocean acoustics or downward and upward in reflection imaging). In the literature on electromagnetic wave propagation, these oppositely propagating waves are often called “bidirectional beams” [1, 2] whereas in the acoustic literature they are usually called “one-way wave fields” [3–7]. In this paper, we use the latter terminology.There is a vast amount of literature on the analytical and numerical aspects of one-way wave propagation [8–13]. A discussion of this is beyond the scope of this paper. Instead, we concentrate on the choice of the decomposition operator and the consequences for reciprocity and representation theorems.Decomposition of a wave field into one-way wave fields is not unique. In particular, the amplitudes of the one-way wave fields can be scaled in different ways. In this paper, we distinguish between the so-called “flux-normalised” and “field-normalised” one-way wave fields. The square of the amplitude of a flux-normalised one-way wave field is by definition the power-flux density (or, for quantum-mechanical waves, the probability-flux density) in the direction of preference. Field-normalised one-way wave fields, on the other hand, are scaled such that the sum of the two oppositely propagating components equals the full wave field. These two forms of normalisation have been briefly analysed by de Hoop [14, 15]. From this analysis, it appeared that the operators for flux-normalised decomposition exhibit more symmetry than the operators for field-normalised decomposition. Exploiting the symmetry of the flux-normalised decomposition operators, the author derived reciprocity and representation theorems for flux-normalised one-way wave fields [16, 17].The first aim of this paper is to discuss flux-normalised versus field-normalised decomposition in a systematic way. In particular, it will be shown that reciprocity theorems for field-normalised one-way wave fields can be derived in a similar way as those for flux-normalised one-way wave fields, even though the operators for field-normalised decomposition exhibit less symmetry.The second aim is to discuss representation theorems for field-normalised one-way wave fields in a systematic way. This discussion includes links to “classical” Kirchhoff-Helmholtz integrals for one-way wave fields as well as to recent single-sided representations for backward propagation, used for example in Marchenko imaging [18]. Despite the links to earlier results, the discussed representations are more general. An advantage of the representations for field-normalised one-way wave fields is that a straightforward summation of the one-way wave fields gives the full wave field.We restrict the discussion to scalar wave fields. In Section2, we formulate a unified scalar wave equation for acoustic waves, horizontally polarised shear waves, transverse electric and transverse magnetic EM waves, and finally quantum-mechanical waves. Next, we reformulate the unified wave equation into a matrix-vector form, discuss symmetry properties of the operator matrix, and use this to derive reciprocity theorems in matrix-vector form. In Section 3, we decompose the matrix-vector wave equation into a coupled system of equations for oppositely propagating one-way wave fields. We separately consider flux normalisation and field normalisation and derive reciprocity theorems for one-way wave fields, using both normalisations. In Section 4, we extensively discuss representation theorems for field-normalised one-way wave fields and indicate applications. We end with conclusions in Section 5.


## 2. Unified Wave Equation and Its Symmetry Properties

### 2.1. Unified Scalar Wave Equation

Using a unified notation, wave propagation in a lossless medium (or, for quantum-mechanical waves, in a lossless potential) is governed by the following two equations in the space-frequency domain:(1)−iωαP+∂jQj=B,(2)−iωβQj+∂jP=Cj.Here,i is the imaginary unit and ω the angular frequency (in this paper, we consider positive frequencies only). Operator ∂j stands for the spatial differential operator ∂/∂xj, and Einstein’s summation convention applies to repeated subscripts. Px,ω and Qjx,ω are space- and frequency-dependent wave field quantities, αx and βx are real-valued space-dependent parameters, and Bx,ω and Cjx,ω are space- and frequency-dependent source distributions. Parameters α and β are both assumed to be positive; hence, metamaterials are not considered in this paper. All quantities are specified in Table 1 for different wave phenomena and are discussed in more detail below. As indicated in the first column of Table 1, we consider 3D and 2D wave problems. For the 3D situation, x=x1,x2,x3 is the 3D Cartesian coordinate vector and lowercase Latin subscripts take on the values 1, 2, and 3. For the 2D situation, x=x1,x3 is the 2D Cartesian coordinate vector and lowercase Latin subscripts take on the values 1 and 3 only.Table 1
Quantities in unified equations (1) and (2).
P
Qj
α
β
B
Cj
Acoustic waves (3D)
p
vj
κ
ρ
q
fj
SH waves (2D)
v2
−τ2j
ρ
1μ
f2
2h2j
TE waves (2D)
E2
−ϵ2jkHk
ε
μ
−J2e
ϵ2jkJkm
TM waves (2D)
H2
ϵ2jkEk
μ
ε
−J2m
−ϵ2jkJke
Quantum waves (3D)
Ψ
2ℏmi∂jΨ
4−4Vℏω
m2ℏωThe unified boundary conditions at an interface between two media with different parameters state thatP and njQj are continuous over the interface. Here, nj represents the components of the normal vector n=n1,n2,n3 at the interface for the 3D situation or n=n1,n3 for the 2D situation.We discuss the quantities in Table1 in more detail. The quantities in row 1, associated to 3D acoustic wave propagation in a lossless fluid medium, are acoustic pressure px,ω, particle velocity vjx,ω, compressibility κx, mass density ρx, volume-injection rate density qx,ω, and external force density fjx,ω. For 2D horizontally polarised shear waves in a lossless solid medium, we have in row 2 horizontal particle velocity v2x,ω, shear stress τ2jx,ω, mass density ρx, shear modulus μx, external horizontal force density f2x,ω, and external shear deformation rate density h2jx,ω. Rows 3 and 4 contain the quantities for 2D electromagnetic wave propagation, with TE standing for transverse electric and TM for transverse magnetic. The quantities are electric field strength Ekx,ω, magnetic field strength Hkx,ω, permittivity εx, permeability μx, external electric current density Jkex,ω, and external magnetic current density Jkmx,ω. Furthermore, ϵijk is the alternating tensor (or Levi-Civita tensor), with ϵ123=ϵ312=ϵ231=1, ϵ213=ϵ321=ϵ132=−1, and all other components being zero. In row 5, the quantities related to 3D quantum-mechanical wave propagation are wave function Ψx,ω, potential Vx, particle mass m, and ℏ=h/2π, with h Planck’s constant.By eliminatingQj from equations (1) and (2), we obtain the unified scalar wave equation
(3)β∂j1β∂jP+k2P=β∂j1βCj+iωβB,with wave number k defined via
(4)k2=αβω2.


### 2.2. Unified Wave Equation in Matrix-Vector Form

We define a configuration with a preferred direction and reorganise equations (1) and (2) accordingly.Consider a 3D spatial domainD, enclosed by surface ∂D. This surface consists of two planar surfaces ∂D0 and ∂D1 perpendicular to the x3-axis and a cylindrical surface ∂Dcyl with its axis parallel to the x3-axis, see Figure 1. The surfaces ∂D0 and ∂D1 are situated at x3=x3,0 and x3=x3,1, respectively, with x3,1>x3,0. In general, these surfaces do not coincide with physical boundaries. Surface S in Figure 1 is a cross section of D at arbitrary x3. The parameters αx and βx are piecewise continuous smoothly varying functions in D, with discontinuous jumps only at interfaces that are perpendicular to the x3-axis (hence, P and Q3 are continuous over the interfaces). In the lateral direction, the domain D can be bounded or unbounded. When D is laterally bounded, the configuration in Figure 1 represents a wave guide. For this situation, we assume that homogeneous Dirichlet or Neumann boundary conditions apply, i.e., P=Q3=0 or nν∂νP=nν∂νQ3=0 at ∂Dcyl, where lowercase Greek subscripts take on the values 1 and 2. When D is laterally unbounded (for example, for reflection imaging applications), the cylindrical surface ∂Dcyl has an infinite radius and we assume that P and Q3 have “sufficient decay” at infinity. For the 2D situation, the configuration is a cross section of the 3D situation for x2=0 and lowercase Greek subscripts take on the value 1 only.Figure 1
Configuration with thex3-direction as the preferred direction. In the lateral direction, this configuration can be bounded (for wave guides) or unbounded (for example, for geophysical reflection imaging applications). For the 2D situation, the configuration is a cross section of the 3D situation for x2=0.We reorganise equations (1) and (2) into a matrix-vector wave equation which acknowledges the x3-direction as the direction of preference. By eliminating the lateral components Q1 and Q2 (or, for 2D wave problems, the lateral component Q1), we obtain [8, 15, 19–21]
(5)∂3q=Aq+d,where wave vector q and source vector d are defined as
(6)q=PQ3,d=C3B0,with
(7)B0=B+1iω∂ν1βCνand operator matrix A defined as
(8)A=0A12A210,with
(9)A12=iωβ,(10)A21=iωα−1iω∂ν1β∂ν.The notation in the right-hand side of equations (7) and (10) should be understood in the sense that differential operators act on all factors to the right of it. Hence, operator ∂ν1/β∂ν, applied via equation (5) to P, stands for ∂ν1/β∂νP.Note that the quantities contained in the wave vectorq are continuous over interfaces perpendicular to the x3-axis. Moreover, these quantities constitute the power-flux density (or, for quantum-mechanical waves, the probability-flux density) in the x3-direction via
(11)j=14P∗Q3+Q3∗P,where the asterisk denotes complex conjugation.


### 2.3. Symmetry Properties of the Operator Matrix

We discuss the symmetry properties of the operator matrixA. First, consider a general operator U (which can be a scalar or a matrix), containing space-dependent parameters and differential operators ∂ν. We introduce the transpose operator Ut via the following integral property:
(12)∫SUftgdxL=∫SftUtgdxL.Here,xL is the lateral coordinate vector, with xL=x1,x2 for 3D and xL=x1 for 2D wave problems. S denotes an integration surface perpendicular to the x3-axis at arbitrary x3, with edge ∂S, see Figure 1. The quantities fxL and gxL are space-dependent test functions (scalars or vectors). When these functions are vectors, ft is the transpose of f; when they are scalars, ft is equal to f. When S is bounded, homogeneous Dirichlet or Neumann conditions are imposed at ∂S. When S is unbounded, ∂S has an infinite radius and fxL and gxL are assumed to have sufficient decay along S towards infinity. Operator U is said to be symmetric when Ut=U and skew-symmetric when Ut=−U. For the special case that U=∂ν, equation (12) implies ∂νt=−∂ν (via integration by parts along S). Hence, operator ∂ν is skew-symmetric.We introduce the adjoint operatorU† (i.e., the complex conjugate transpose of U) via the integral property
(13)∫SUf†gdxL=∫Sf†U†gdxL.When the test functions are vectors,f† is the complex conjugate transpose of f; when they are scalars, f† is the complex conjugate of f. Operator U is said to be Hermitian (or self-adjoint) when U†=U and skew-Hermitian when U†=−U. For the operators A12 and A21, defined in equations (9) and (10), we find A12t=A12, A21t=A21, A12†=−A12, and A21†=−A21. Hence, operators A12 and A21 are symmetric and skew-Hermitian. With these relations, we find for the operator matrix A
(14)AtN=−NA,(15)A†K=−KA,with
(16)N=01−10,K=0110.Note that, using the expressions forq and K in equations (6) and (16), we can rewrite equation (11) for the power-flux density (or, for quantum-mechanical waves, the probability-flux density) as
(17)j=14q†Kq.


### 2.4. Reciprocity Theorems

We derive reciprocity theorems between two independent solutions of wave equation (5) for the configuration of Figure 1. We consider two states A and B, characterised by wave vectors qAx,ω and qBx,ω, obeying wave equation (5), with source vectors dAx,ω and dBx,ω. In domain D, the parameters α and β, and hence the matrix operator A, are chosen the same in the two states (outside ∂D they may be different in the two states). Consider the quantity ∂3qAtNqB in domain D. Applying the product rule for differentiation, using equation (5) for both states, integrating the result over D and applying the theorem of Gauss yields
(18)∫DAqAt+dAtNqB+qAtNAqB+dBdx=∫∂DqAtNqBn3dx.Here,n3 is the component parallel to the x3-axis of the outward pointing normal vector on ∂D, with n3=−1 at ∂D0, n3=+1 at ∂D1, and n3=0 at ∂Dcyl, see Figure 1. In the following, the integral on the right-hand side is restricted to the horizontal surfaces ∂D0 and ∂D1, which together are denoted by ∂D0,1. The integral on the left-hand side can be written as ∫D⋯dx=∫x3,0x3,1dx3∫S⋯dxL. Using equation (12) for the integral along S and symmetry property (14), it follows that the two terms in equation (18) containing operator A cancel each other. Hence, we are left with
(19)∫DdAtNqB+qAtNdBdx=∫∂D0,1qAtNqBn3dxL.This is a convolution-type reciprocity theorem [22–24], because products like qAtx,ωNqBx,ω in the frequency domain correspond to convolutions in the time domain. A more familiar form is obtained by substituting the expressions for q, d, and N (equations (6) and (16)), choosing Cj=0 and using equation (2) to eliminate Q3, which gives
(20)∫D−BAPB+PABBdx=∫∂D0,11iωβPA∂3PB−∂3PAPBn3dxL.Next, consider the quantity∂3qA†KqB in domain D. Following the same steps as above, using equations (13) and (15) instead of (12) and (14), we obtain
(21)∫DdA†KqB+qA†KdBdx=∫∂D0,1qA†KqBn3dxL.This is a correlation-type reciprocity theorem [25], because products like qA†x,ωKqBx,ω in the frequency domain correspond to correlations in the time domain. Substituting the expressions for q, d, and K and choosing Cj=0 yield the more familiar form
(22)∫DBA∗PB+PA∗BBdx=∫∂D0,11iωβPA∗∂3PB−∂3PA∗PBn3dxL.We obtain a special case by choosing statesA and B identical. Dropping the subscripts A and B in equations (21) and (22) and multiplying both sides of these equations by 1/4 give
(23)14∫Dd†Kq+q†Kddx=14∫∂D0,1q†Kqn3dxL,(24)14∫DB∗P+P∗Bdx=14∫∂D0,11iωβP∗∂3P−∂3P∗Pn3dxL,respectively. These equations quantify conservation of power (or, for quantum-mechanical waves, probability).


## 2.1. Unified Scalar Wave Equation

Using a unified notation, wave propagation in a lossless medium (or, for quantum-mechanical waves, in a lossless potential) is governed by the following two equations in the space-frequency domain:(1)−iωαP+∂jQj=B,(2)−iωβQj+∂jP=Cj.Here,i is the imaginary unit and ω the angular frequency (in this paper, we consider positive frequencies only). Operator ∂j stands for the spatial differential operator ∂/∂xj, and Einstein’s summation convention applies to repeated subscripts. Px,ω and Qjx,ω are space- and frequency-dependent wave field quantities, αx and βx are real-valued space-dependent parameters, and Bx,ω and Cjx,ω are space- and frequency-dependent source distributions. Parameters α and β are both assumed to be positive; hence, metamaterials are not considered in this paper. All quantities are specified in Table 1 for different wave phenomena and are discussed in more detail below. As indicated in the first column of Table 1, we consider 3D and 2D wave problems. For the 3D situation, x=x1,x2,x3 is the 3D Cartesian coordinate vector and lowercase Latin subscripts take on the values 1, 2, and 3. For the 2D situation, x=x1,x3 is the 2D Cartesian coordinate vector and lowercase Latin subscripts take on the values 1 and 3 only.Table 1
Quantities in unified equations (1) and (2).
P
Qj
α
β
B
Cj
Acoustic waves (3D)
p
vj
κ
ρ
q
fj
SH waves (2D)
v2
−τ2j
ρ
1μ
f2
2h2j
TE waves (2D)
E2
−ϵ2jkHk
ε
μ
−J2e
ϵ2jkJkm
TM waves (2D)
H2
ϵ2jkEk
μ
ε
−J2m
−ϵ2jkJke
Quantum waves (3D)
Ψ
2ℏmi∂jΨ
4−4Vℏω
m2ℏωThe unified boundary conditions at an interface between two media with different parameters state thatP and njQj are continuous over the interface. Here, nj represents the components of the normal vector n=n1,n2,n3 at the interface for the 3D situation or n=n1,n3 for the 2D situation.We discuss the quantities in Table1 in more detail. The quantities in row 1, associated to 3D acoustic wave propagation in a lossless fluid medium, are acoustic pressure px,ω, particle velocity vjx,ω, compressibility κx, mass density ρx, volume-injection rate density qx,ω, and external force density fjx,ω. For 2D horizontally polarised shear waves in a lossless solid medium, we have in row 2 horizontal particle velocity v2x,ω, shear stress τ2jx,ω, mass density ρx, shear modulus μx, external horizontal force density f2x,ω, and external shear deformation rate density h2jx,ω. Rows 3 and 4 contain the quantities for 2D electromagnetic wave propagation, with TE standing for transverse electric and TM for transverse magnetic. The quantities are electric field strength Ekx,ω, magnetic field strength Hkx,ω, permittivity εx, permeability μx, external electric current density Jkex,ω, and external magnetic current density Jkmx,ω. Furthermore, ϵijk is the alternating tensor (or Levi-Civita tensor), with ϵ123=ϵ312=ϵ231=1, ϵ213=ϵ321=ϵ132=−1, and all other components being zero. In row 5, the quantities related to 3D quantum-mechanical wave propagation are wave function Ψx,ω, potential Vx, particle mass m, and ℏ=h/2π, with h Planck’s constant.By eliminatingQj from equations (1) and (2), we obtain the unified scalar wave equation
(3)β∂j1β∂jP+k2P=β∂j1βCj+iωβB,with wave number k defined via
(4)k2=αβω2.


## 2.2. Unified Wave Equation in Matrix-Vector Form

We define a configuration with a preferred direction and reorganise equations (1) and (2) accordingly.Consider a 3D spatial domainD, enclosed by surface ∂D. This surface consists of two planar surfaces ∂D0 and ∂D1 perpendicular to the x3-axis and a cylindrical surface ∂Dcyl with its axis parallel to the x3-axis, see Figure 1. The surfaces ∂D0 and ∂D1 are situated at x3=x3,0 and x3=x3,1, respectively, with x3,1>x3,0. In general, these surfaces do not coincide with physical boundaries. Surface S in Figure 1 is a cross section of D at arbitrary x3. The parameters αx and βx are piecewise continuous smoothly varying functions in D, with discontinuous jumps only at interfaces that are perpendicular to the x3-axis (hence, P and Q3 are continuous over the interfaces). In the lateral direction, the domain D can be bounded or unbounded. When D is laterally bounded, the configuration in Figure 1 represents a wave guide. For this situation, we assume that homogeneous Dirichlet or Neumann boundary conditions apply, i.e., P=Q3=0 or nν∂νP=nν∂νQ3=0 at ∂Dcyl, where lowercase Greek subscripts take on the values 1 and 2. When D is laterally unbounded (for example, for reflection imaging applications), the cylindrical surface ∂Dcyl has an infinite radius and we assume that P and Q3 have “sufficient decay” at infinity. For the 2D situation, the configuration is a cross section of the 3D situation for x2=0 and lowercase Greek subscripts take on the value 1 only.Figure 1
Configuration with thex3-direction as the preferred direction. In the lateral direction, this configuration can be bounded (for wave guides) or unbounded (for example, for geophysical reflection imaging applications). For the 2D situation, the configuration is a cross section of the 3D situation for x2=0.We reorganise equations (1) and (2) into a matrix-vector wave equation which acknowledges the x3-direction as the direction of preference. By eliminating the lateral components Q1 and Q2 (or, for 2D wave problems, the lateral component Q1), we obtain [8, 15, 19–21]
(5)∂3q=Aq+d,where wave vector q and source vector d are defined as
(6)q=PQ3,d=C3B0,with
(7)B0=B+1iω∂ν1βCνand operator matrix A defined as
(8)A=0A12A210,with
(9)A12=iωβ,(10)A21=iωα−1iω∂ν1β∂ν.The notation in the right-hand side of equations (7) and (10) should be understood in the sense that differential operators act on all factors to the right of it. Hence, operator ∂ν1/β∂ν, applied via equation (5) to P, stands for ∂ν1/β∂νP.Note that the quantities contained in the wave vectorq are continuous over interfaces perpendicular to the x3-axis. Moreover, these quantities constitute the power-flux density (or, for quantum-mechanical waves, the probability-flux density) in the x3-direction via
(11)j=14P∗Q3+Q3∗P,where the asterisk denotes complex conjugation.


## 2.3. Symmetry Properties of the Operator Matrix

We discuss the symmetry properties of the operator matrixA. First, consider a general operator U (which can be a scalar or a matrix), containing space-dependent parameters and differential operators ∂ν. We introduce the transpose operator Ut via the following integral property:
(12)∫SUftgdxL=∫SftUtgdxL.Here,xL is the lateral coordinate vector, with xL=x1,x2 for 3D and xL=x1 for 2D wave problems. S denotes an integration surface perpendicular to the x3-axis at arbitrary x3, with edge ∂S, see Figure 1. The quantities fxL and gxL are space-dependent test functions (scalars or vectors). When these functions are vectors, ft is the transpose of f; when they are scalars, ft is equal to f. When S is bounded, homogeneous Dirichlet or Neumann conditions are imposed at ∂S. When S is unbounded, ∂S has an infinite radius and fxL and gxL are assumed to have sufficient decay along S towards infinity. Operator U is said to be symmetric when Ut=U and skew-symmetric when Ut=−U. For the special case that U=∂ν, equation (12) implies ∂νt=−∂ν (via integration by parts along S). Hence, operator ∂ν is skew-symmetric.We introduce the adjoint operatorU† (i.e., the complex conjugate transpose of U) via the integral property
(13)∫SUf†gdxL=∫Sf†U†gdxL.When the test functions are vectors,f† is the complex conjugate transpose of f; when they are scalars, f† is the complex conjugate of f. Operator U is said to be Hermitian (or self-adjoint) when U†=U and skew-Hermitian when U†=−U. For the operators A12 and A21, defined in equations (9) and (10), we find A12t=A12, A21t=A21, A12†=−A12, and A21†=−A21. Hence, operators A12 and A21 are symmetric and skew-Hermitian. With these relations, we find for the operator matrix A
(14)AtN=−NA,(15)A†K=−KA,with
(16)N=01−10,K=0110.Note that, using the expressions forq and K in equations (6) and (16), we can rewrite equation (11) for the power-flux density (or, for quantum-mechanical waves, the probability-flux density) as
(17)j=14q†Kq.


## 2.4. Reciprocity Theorems

We derive reciprocity theorems between two independent solutions of wave equation (5) for the configuration of Figure 1. We consider two states A and B, characterised by wave vectors qAx,ω and qBx,ω, obeying wave equation (5), with source vectors dAx,ω and dBx,ω. In domain D, the parameters α and β, and hence the matrix operator A, are chosen the same in the two states (outside ∂D they may be different in the two states). Consider the quantity ∂3qAtNqB in domain D. Applying the product rule for differentiation, using equation (5) for both states, integrating the result over D and applying the theorem of Gauss yields
(18)∫DAqAt+dAtNqB+qAtNAqB+dBdx=∫∂DqAtNqBn3dx.Here,n3 is the component parallel to the x3-axis of the outward pointing normal vector on ∂D, with n3=−1 at ∂D0, n3=+1 at ∂D1, and n3=0 at ∂Dcyl, see Figure 1. In the following, the integral on the right-hand side is restricted to the horizontal surfaces ∂D0 and ∂D1, which together are denoted by ∂D0,1. The integral on the left-hand side can be written as ∫D⋯dx=∫x3,0x3,1dx3∫S⋯dxL. Using equation (12) for the integral along S and symmetry property (14), it follows that the two terms in equation (18) containing operator A cancel each other. Hence, we are left with
(19)∫DdAtNqB+qAtNdBdx=∫∂D0,1qAtNqBn3dxL.This is a convolution-type reciprocity theorem [22–24], because products like qAtx,ωNqBx,ω in the frequency domain correspond to convolutions in the time domain. A more familiar form is obtained by substituting the expressions for q, d, and N (equations (6) and (16)), choosing Cj=0 and using equation (2) to eliminate Q3, which gives
(20)∫D−BAPB+PABBdx=∫∂D0,11iωβPA∂3PB−∂3PAPBn3dxL.Next, consider the quantity∂3qA†KqB in domain D. Following the same steps as above, using equations (13) and (15) instead of (12) and (14), we obtain
(21)∫DdA†KqB+qA†KdBdx=∫∂D0,1qA†KqBn3dxL.This is a correlation-type reciprocity theorem [25], because products like qA†x,ωKqBx,ω in the frequency domain correspond to correlations in the time domain. Substituting the expressions for q, d, and K and choosing Cj=0 yield the more familiar form
(22)∫DBA∗PB+PA∗BBdx=∫∂D0,11iωβPA∗∂3PB−∂3PA∗PBn3dxL.We obtain a special case by choosing statesA and B identical. Dropping the subscripts A and B in equations (21) and (22) and multiplying both sides of these equations by 1/4 give
(23)14∫Dd†Kq+q†Kddx=14∫∂D0,1q†Kqn3dxL,(24)14∫DB∗P+P∗Bdx=14∫∂D0,11iωβP∗∂3P−∂3P∗Pn3dxL,respectively. These equations quantify conservation of power (or, for quantum-mechanical waves, probability).


## 3. Decomposed Wave Equation and Its Symmetry Properties

### 3.1. General Decomposition of the Matrix-Vector Wave Equation

To facilitate the decomposition of the matrix-vector wave equation (equation (5)), we recast the operator matrix A into a somewhat different form. To this end, we introduce an operator H2, according to
(25)H2=−iωβA21β=k2+β∂ν1β∂νβ,with operator A21 defined in equation (10) and wavenumber k in equation (4). Operator H2 can be rewritten as a Helmholtz operator [14, 21]
(26)H2=ks2+∂ν∂ν,with the scaled wavenumber ks defined as [26]
(27)ks2=k2−3∂νβ∂νβ4β2+∂ν∂νβ2β.Note thatH2t=H2 and H2†=H2; hence, operator H2 is symmetric and self-adjoint and its spectrum is real-valued (with positive and negative eigenvalues). Using equation (25), we rewrite operator matrix A, defined in equation (8), as
(28)A=0iωβ−1iωβH21β0.Next, we decompose this operator matrix as follows(29)A=LHL−1,with
(30)H=iH100−iH1,(31)L=L1L1L2−L2,(32)L−1=12L1−1L2−1L1−1−L2−1.OperatorsH1, L1, and L2 are pseudodifferential operators [7, 8, 14, 16, 21, 27–30]. The decomposition expressed by equation (29) is not unique; hence, different choices for operators H1, L1, and L2 are possible. We discuss two of these choices in detail in the next two sections. Here, we derive some general relations that are independent of these choices.By substituting equations (28), (30), (31), and (32) into equation (29), we obtain the following relations
(33)ωβ=L1H1L2−1,(34)1ωβH21β=L2H1L1−1.We introduce a decomposed field vectorp and a decomposed source vector s via
(35)q=Lp,p=L−1q,(36)d=Ls,s=L−1d,where
(37)p=P+P−,s=S+S−.Substitution of equations (29), (35), and (36) into the matrix-vector wave equation (5) yields
(38)∂3p=H−L−1∂3Lp+s.Substituting equations (30), (31), (32), and (37) into equation (38) gives
(39)∂3P+P−=iH100−iH1P+P−−12L1−1L2−1L1−1−L2−1∂3L1∂3L1∂3L2−∂3L2P+P−+S+S−.This is a system of coupled one-way wave equations. From the first term on the right-hand side, it follows that the one-way wave fieldsP+ and P− propagate in the positive and negative x3-direction, respectively. The second term on the right-hand side accounts for coupling between P+ and P−. The last term on the right-hand side contains sources S+ and S− which emit waves in the positive and negative x3-direction, respectively.We conclude this section by substituting equations (35) and (36) into equations (19), (21), and (23). Using equations (12) and (13) for the integration along the lateral coordinates, this yields
(40)∫DsAtLtNLpB+pAtLtNLsBdx=∫∂D0,1pAtLtNLpBn3dxL,(41)∫DsA†L†KLpB+pA†L†KLsBdx=∫∂D0,1pA†L†KLpBn3dxL,(42)14∫Ds†L†KLp+p†L†KLsdx=14∫∂D0,1p†L†KLpn3dxL.These equations form the basis for reciprocity theorems for the decomposed field and source vectorsp and s in the next two sections.


### 3.2. Flux-Normalised Decomposition and Reciprocity Theorems

The first choice of operatorsH1, L1, and L2 obeying equations (33) and (34) is [14–16]
(43)H1=H21/2,(44)L1=ω/21/2β1/2H1−1/2,(45)L2=2ω−1/2β−1/2H11/2.OperatorH1, which is the square root of the Helmholtz operator H2, is commonly known as the square root operator [3, 4, 8]. Like the Helmholtz operator H2, the square root operator H1 is a symmetric operator [16], hence H1t=H1. For the adjoint square root operator, we have H1†=H1t∗=H1∗. The spectrum of H1 is real-valued for propagating waves and imaginary-valued for evanescent waves. Hence, unlike the Helmholtz operator, the square root operator is not self-adjoint. If we neglect evanescent waves, we may approximate the adjoint square root operator as H1†≈H1. Similar relations hold for the square root of the square root operator and its inverse; hence, H1±1/2t=H1±1/2, and neglecting evanescent waves, H1±1/2†≈H1±1/2. From here onward, we replace ≈ by = when the only approximation is the negligence of evanescent waves. Using these symmetry relations for H1 and equations (16), (31), (44), and (45), we obtain
(46)LtNL=−N,and neglecting evanescent waves,
(47)L†KL=J,with
(48)J=100−1.Hence, equations (40), (41), and (42) simplify to
(49)−∫DsAtNpB+pAtNsBdx=−∫∂D0,1pAtNpBn3dxL,∫DsA†JpB+pA†JsBdx=∫∂D0,1pA†JpBn3dxL,14∫Ds†Jp+p†Jsdx=14∫∂D0,1p†Jpn3dxL.By substituting the expressions forp, s, N, and J (equations (37), (16), and (48)), we obtain
(50)∫D−SA+PB−+SA−PB+−PA+SB−+PA−SB+dx=∫∂D0,1−PA+PB−+PA−PB+n3dxL,(51)∫DSA+∗PB+−SA−∗PB−+PA+∗SB+−PA−∗SB−dx=∫∂D0,1PA+∗PB+−PA−∗PB−n3dxL,(52)14∫DS+∗P+−S−∗P−+P+∗S+−P−∗S−dx=14∫∂D0,1P+2−P−2n3dxL.Note that, since the right-hand side of equation (52) is equal to the right-hand side of equation (24), it quantifies the power flux (or the probability flux for quantum-mechanical waves) through the surface ∂D0,1. Therefore, we call P+ and P− flux-normalised one-way wave fields. Consequently, equations (50) and (51) are reciprocity theorems of the convolution type and correlation type, respectively, for flux-normalised one-way wave fields. These theorems have been derived previously [16] and have found applications in advanced wave field imaging methods for active and passive data [31–42].


### 3.3. Field-Normalised Decomposition and Reciprocity Theorems

The second choice of operatorsH1, L1, and L2 obeying equations (33) and (34) is [21]
(53)H1=β1/2H21/2β−1/2,(54)L1=1,(55)L2=ωβ−1H1.Only the Helmholtz operatorH2 is the same as in the previous section (it is defined in equation (26)). The operators H1, L1, and L2 are different from those in the previous section, but for convenience, we use the same symbols. Using q=Lp (equation (35)) and equations (6), (31), (37), and (54), we find
(56)P=P++P−.Hence, P+ and P− have the same physical dimension as the full field variable P (which is defined in Table 1 for different wave phenomena). Therefore, we call P+ and P− field-normalised one-way wave fields (for convenience, we use the same symbols as in the previous section).The square root operatorH21/2 is symmetric, but H1 defined in equation (53) is not. From this equation, it easily follows that H1 premultiplied by β−1 is symmetric, hence
(57)1βH1t=1βH1,and neglecting evanescent waves,
(58)1βH1†=1βH1.Using these symmetry relations for1/βH1 and equations (16), (31), (54), and (55), we obtain
(59)LtNL=0−2L22L20=−N2ωβH1=−2ωβH1tN,and neglecting evanescent waves,
(60)L†KL=2L200−2L2=J2ωβH1=2ωβH1†J.Using this in equations (40) and (41) yields
(61)−∫DsAt2ωβH1tNpB+pAtN2ωβH1sBdx=−∫∂D0,1pAt2ωβH1tNpBn3dxL,∫DsA†2ωβH1†JpB+pA†J2ωβH1sBdx=∫∂D0,1pA†2ωβH1†JpBn3dxL.By substituting the expressions forp, s, N, and J (equations (37), (16), and (48)), using equations (12) and (13), we obtain
(62)−∫D2ωβH1SA+PB−−H1SA−PB++PA+H1SB−−PA−H1SB+dx=−∫∂D0.12ωβH1PA+PB−−H1PA−PB+n3dxL,(63)∫D2ωβH1SA+∗PB+−H1SA−∗PB−+PA+∗H1SB+−PA−∗H1SB−dx=∫∂D0,12ωβH1PA+∗PB+−H1PA−∗PB−n3dxL.We aim to remove the operatorH1 from these equations. From equations (39) and (54), we obtain
(64)∂3P+=+iH1P+−12L2−1∂3L2P+−P−+S+,(65)∂3P−=−iH1P−+12L2−1∂3L2P+−P−+S−,with L2 defined in equation (55). Assuming that in state A the derivatives in the x3-direction of the parameters α and β at ∂D0,1 vanish and there are no sources at ∂D0,1, we find from equations (64) and (65)
(66)∂3PA±=±iH1PA±at∂D0,1.Below we use this to removeH1 from the right-hand sides of equations (62) and (63). Next, we aim to remove H1 from the left-hand sides of these equations. From s=L−1d (equation (36)) and equations (6), (32), (37), (54), and (55), we find
(67)S±=±121ωβH1−1B0+12C3,or
(68)±2ωβH1S±=B0±1ωβH1C3.We define new decomposed sourcesB0+ and B0−, according to
(69)B0±=B0±1ωβH1C3=±2ωβH1S±.Using equations (66) and (69) in the right- and left-hand sides of equations (62) and (63), we obtain
(70)∫D−B0,A+PB−−B0,A−PB++PA+B0,B−+PA−B0,B+dx=∫∂D0,1−2iωβ∂3PA+PB−+∂3PA−PB+n3dxL,(71)∫DB0,A+∗PB++B0,A−∗PB−+PA+∗B0,B++PA−∗B0,B−dx=∫∂D0,1−2iωβ∂3PA+∗PB++∂3PA−∗PB−n3dxL.Equations (70) and (71) are reciprocity theorems of the convolution type and correlation type, respectively, for field-normalised one-way wave fields. These theorems are modifications of previously obtained results [43, 44]. The main modification is that we applied decomposition at both sides of the equations instead of at the right-hand sides only. Moreover, in the present derivation, the condition for the validity of equation (66) is only imposed for state A. In the next section, we use equations (70) and (71) to derive representation theorems for field-normalised one-way wave fields and we indicate applications.


## 3.1. General Decomposition of the Matrix-Vector Wave Equation

To facilitate the decomposition of the matrix-vector wave equation (equation (5)), we recast the operator matrix A into a somewhat different form. To this end, we introduce an operator H2, according to
(25)H2=−iωβA21β=k2+β∂ν1β∂νβ,with operator A21 defined in equation (10) and wavenumber k in equation (4). Operator H2 can be rewritten as a Helmholtz operator [14, 21]
(26)H2=ks2+∂ν∂ν,with the scaled wavenumber ks defined as [26]
(27)ks2=k2−3∂νβ∂νβ4β2+∂ν∂νβ2β.Note thatH2t=H2 and H2†=H2; hence, operator H2 is symmetric and self-adjoint and its spectrum is real-valued (with positive and negative eigenvalues). Using equation (25), we rewrite operator matrix A, defined in equation (8), as
(28)A=0iωβ−1iωβH21β0.Next, we decompose this operator matrix as follows(29)A=LHL−1,with
(30)H=iH100−iH1,(31)L=L1L1L2−L2,(32)L−1=12L1−1L2−1L1−1−L2−1.OperatorsH1, L1, and L2 are pseudodifferential operators [7, 8, 14, 16, 21, 27–30]. The decomposition expressed by equation (29) is not unique; hence, different choices for operators H1, L1, and L2 are possible. We discuss two of these choices in detail in the next two sections. Here, we derive some general relations that are independent of these choices.By substituting equations (28), (30), (31), and (32) into equation (29), we obtain the following relations
(33)ωβ=L1H1L2−1,(34)1ωβH21β=L2H1L1−1.We introduce a decomposed field vectorp and a decomposed source vector s via
(35)q=Lp,p=L−1q,(36)d=Ls,s=L−1d,where
(37)p=P+P−,s=S+S−.Substitution of equations (29), (35), and (36) into the matrix-vector wave equation (5) yields
(38)∂3p=H−L−1∂3Lp+s.Substituting equations (30), (31), (32), and (37) into equation (38) gives
(39)∂3P+P−=iH100−iH1P+P−−12L1−1L2−1L1−1−L2−1∂3L1∂3L1∂3L2−∂3L2P+P−+S+S−.This is a system of coupled one-way wave equations. From the first term on the right-hand side, it follows that the one-way wave fieldsP+ and P− propagate in the positive and negative x3-direction, respectively. The second term on the right-hand side accounts for coupling between P+ and P−. The last term on the right-hand side contains sources S+ and S− which emit waves in the positive and negative x3-direction, respectively.We conclude this section by substituting equations (35) and (36) into equations (19), (21), and (23). Using equations (12) and (13) for the integration along the lateral coordinates, this yields
(40)∫DsAtLtNLpB+pAtLtNLsBdx=∫∂D0,1pAtLtNLpBn3dxL,(41)∫DsA†L†KLpB+pA†L†KLsBdx=∫∂D0,1pA†L†KLpBn3dxL,(42)14∫Ds†L†KLp+p†L†KLsdx=14∫∂D0,1p†L†KLpn3dxL.These equations form the basis for reciprocity theorems for the decomposed field and source vectorsp and s in the next two sections.


## 3.2. Flux-Normalised Decomposition and Reciprocity Theorems

The first choice of operatorsH1, L1, and L2 obeying equations (33) and (34) is [14–16]
(43)H1=H21/2,(44)L1=ω/21/2β1/2H1−1/2,(45)L2=2ω−1/2β−1/2H11/2.OperatorH1, which is the square root of the Helmholtz operator H2, is commonly known as the square root operator [3, 4, 8]. Like the Helmholtz operator H2, the square root operator H1 is a symmetric operator [16], hence H1t=H1. For the adjoint square root operator, we have H1†=H1t∗=H1∗. The spectrum of H1 is real-valued for propagating waves and imaginary-valued for evanescent waves. Hence, unlike the Helmholtz operator, the square root operator is not self-adjoint. If we neglect evanescent waves, we may approximate the adjoint square root operator as H1†≈H1. Similar relations hold for the square root of the square root operator and its inverse; hence, H1±1/2t=H1±1/2, and neglecting evanescent waves, H1±1/2†≈H1±1/2. From here onward, we replace ≈ by = when the only approximation is the negligence of evanescent waves. Using these symmetry relations for H1 and equations (16), (31), (44), and (45), we obtain
(46)LtNL=−N,and neglecting evanescent waves,
(47)L†KL=J,with
(48)J=100−1.Hence, equations (40), (41), and (42) simplify to
(49)−∫DsAtNpB+pAtNsBdx=−∫∂D0,1pAtNpBn3dxL,∫DsA†JpB+pA†JsBdx=∫∂D0,1pA†JpBn3dxL,14∫Ds†Jp+p†Jsdx=14∫∂D0,1p†Jpn3dxL.By substituting the expressions forp, s, N, and J (equations (37), (16), and (48)), we obtain
(50)∫D−SA+PB−+SA−PB+−PA+SB−+PA−SB+dx=∫∂D0,1−PA+PB−+PA−PB+n3dxL,(51)∫DSA+∗PB+−SA−∗PB−+PA+∗SB+−PA−∗SB−dx=∫∂D0,1PA+∗PB+−PA−∗PB−n3dxL,(52)14∫DS+∗P+−S−∗P−+P+∗S+−P−∗S−dx=14∫∂D0,1P+2−P−2n3dxL.Note that, since the right-hand side of equation (52) is equal to the right-hand side of equation (24), it quantifies the power flux (or the probability flux for quantum-mechanical waves) through the surface ∂D0,1. Therefore, we call P+ and P− flux-normalised one-way wave fields. Consequently, equations (50) and (51) are reciprocity theorems of the convolution type and correlation type, respectively, for flux-normalised one-way wave fields. These theorems have been derived previously [16] and have found applications in advanced wave field imaging methods for active and passive data [31–42].


## 3.3. Field-Normalised Decomposition and Reciprocity Theorems

The second choice of operatorsH1, L1, and L2 obeying equations (33) and (34) is [21]
(53)H1=β1/2H21/2β−1/2,(54)L1=1,(55)L2=ωβ−1H1.Only the Helmholtz operatorH2 is the same as in the previous section (it is defined in equation (26)). The operators H1, L1, and L2 are different from those in the previous section, but for convenience, we use the same symbols. Using q=Lp (equation (35)) and equations (6), (31), (37), and (54), we find
(56)P=P++P−.Hence, P+ and P− have the same physical dimension as the full field variable P (which is defined in Table 1 for different wave phenomena). Therefore, we call P+ and P− field-normalised one-way wave fields (for convenience, we use the same symbols as in the previous section).The square root operatorH21/2 is symmetric, but H1 defined in equation (53) is not. From this equation, it easily follows that H1 premultiplied by β−1 is symmetric, hence
(57)1βH1t=1βH1,and neglecting evanescent waves,
(58)1βH1†=1βH1.Using these symmetry relations for1/βH1 and equations (16), (31), (54), and (55), we obtain
(59)LtNL=0−2L22L20=−N2ωβH1=−2ωβH1tN,and neglecting evanescent waves,
(60)L†KL=2L200−2L2=J2ωβH1=2ωβH1†J.Using this in equations (40) and (41) yields
(61)−∫DsAt2ωβH1tNpB+pAtN2ωβH1sBdx=−∫∂D0,1pAt2ωβH1tNpBn3dxL,∫DsA†2ωβH1†JpB+pA†J2ωβH1sBdx=∫∂D0,1pA†2ωβH1†JpBn3dxL.By substituting the expressions forp, s, N, and J (equations (37), (16), and (48)), using equations (12) and (13), we obtain
(62)−∫D2ωβH1SA+PB−−H1SA−PB++PA+H1SB−−PA−H1SB+dx=−∫∂D0.12ωβH1PA+PB−−H1PA−PB+n3dxL,(63)∫D2ωβH1SA+∗PB+−H1SA−∗PB−+PA+∗H1SB+−PA−∗H1SB−dx=∫∂D0,12ωβH1PA+∗PB+−H1PA−∗PB−n3dxL.We aim to remove the operatorH1 from these equations. From equations (39) and (54), we obtain
(64)∂3P+=+iH1P+−12L2−1∂3L2P+−P−+S+,(65)∂3P−=−iH1P−+12L2−1∂3L2P+−P−+S−,with L2 defined in equation (55). Assuming that in state A the derivatives in the x3-direction of the parameters α and β at ∂D0,1 vanish and there are no sources at ∂D0,1, we find from equations (64) and (65)
(66)∂3PA±=±iH1PA±at∂D0,1.Below we use this to removeH1 from the right-hand sides of equations (62) and (63). Next, we aim to remove H1 from the left-hand sides of these equations. From s=L−1d (equation (36)) and equations (6), (32), (37), (54), and (55), we find
(67)S±=±121ωβH1−1B0+12C3,or
(68)±2ωβH1S±=B0±1ωβH1C3.We define new decomposed sourcesB0+ and B0−, according to
(69)B0±=B0±1ωβH1C3=±2ωβH1S±.Using equations (66) and (69) in the right- and left-hand sides of equations (62) and (63), we obtain
(70)∫D−B0,A+PB−−B0,A−PB++PA+B0,B−+PA−B0,B+dx=∫∂D0,1−2iωβ∂3PA+PB−+∂3PA−PB+n3dxL,(71)∫DB0,A+∗PB++B0,A−∗PB−+PA+∗B0,B++PA−∗B0,B−dx=∫∂D0,1−2iωβ∂3PA+∗PB++∂3PA−∗PB−n3dxL.Equations (70) and (71) are reciprocity theorems of the convolution type and correlation type, respectively, for field-normalised one-way wave fields. These theorems are modifications of previously obtained results [43, 44]. The main modification is that we applied decomposition at both sides of the equations instead of at the right-hand sides only. Moreover, in the present derivation, the condition for the validity of equation (66) is only imposed for state A. In the next section, we use equations (70) and (71) to derive representation theorems for field-normalised one-way wave fields and we indicate applications.


## 4. Field-Normalised Representation Theorems

### 4.1. Green’s Functions

Representation theorems are obtained by substituting Green’s functions in reciprocity theorems. Our aim is to introduce one-way Green’s functions, to be used in the reciprocity theorems for field-normalised one-way wave fields (equations (70) and (71)). First, we introduce the full Green’s function Gx,xA,ω as a solution of the unified wave equation (3) for a unit monopole point source at xA, with Bx,ω=δx−xA and Cjx,ω=0. Hence,
(72)β∂j1β∂jG+k2G=iωβδx−xA.As boundary condition, we impose the radiation condition (i.e., outward propagating waves at infinity). Next, we introduce one-way Green’s function as solutions of the coupled one-way equations (64) and (65) for a unit monopole point source at xA. Hence, we choose again Bx,ω=δx−xA and Cjx,ω=0. Using equations (69) and (7), we define decomposed sources as B0±=B±=B=±2L2S±, with L2 defined in equation (55), or
(73)S±x,ω=±12L2−1B±x,ω=±12L2−1Bx,ω=±12L2−1δx−xA.We consider two sets of one-way Green’s functions. For the first set, we choose a point sourceS+x,ω=1/2L2−1B+x,ω, with B+x,ω=δx−xA, which emits waves from xA in the positive x3-direction, and we set S−x,ω equal to zero. Hence, for this first set, one-way equations (64) and (65) become
(74)∂3G+,+=+iH1G+,+−12L2−1∂3L2G+,+−G−,++12L2−1δx−xA,(75)∂3G−,+=−iH1G−,++12L2−1∂3L2G+,+−G−,+.Here,G±,+ stands for G±,+x,xA,ω. The second superscript (+) indicates that the source at xA emits waves in the positive x3-direction. The first superscript (±) denotes the propagation direction at x. For the second set of one-way Green’s functions, we choose a point source S−x,ω=−1/2L2−1B−x,ω, with B−x,ω=δx−xA, which emits waves from xA in the negative x3-direction, and we set S+x,ω equal to zero. Hence, for this second set, one-way equations (64) and (65) become
(76)∂3G+,−=+iH1G+,−−12L2−1∂3L2G+,−−G−,−,(77)∂3G−,−=−iH1G−,−+12L2−1∂3L2G+,−−G−,−−12L2−1δx−xA.Here,G±,− stands for G±,−x,xA,ω, with the second superscript (−) indicating that the source at xA emits waves in the negative x3-direction. Like for the full Green’s function Gx,xA,ω, we impose radiation conditions for both sets of one-way Green’s functions.To find a relation between the full Green’s function and the one-way Green’s functions, we evaluateβ∂31/β∂3G+,++G−,++G+,−+G−,− using equations (74), (75), (76), (77), (25), (53), and (55). This gives equation (72), with G replaced by G+,++G−,++G+,−+G−,−. Since the full Green’s function and the one-way Green’s functions obey the same radiation conditions, we thus find
(78)G=G+,++G−,++G+,−+G−,−.This very simple relation is a consequence of the field-normalised decomposition, introduced in Section3.3.


### 4.2. Source-Receiver Reciprocity

We derive source-receiver reciprocity relations for the field-normalised one-way Green’s functions introduced in the previous section. To this end, we make use of the reciprocity theorem of the convolution type for field-normalised one-way wave fields (equation (70)). This theorem was derived for the configuration of Figure 1, assuming that in domain D, the parameters α and β are the same in the two states (see Section 2.4). Outside D, these parameters may be different in the two states. For the Green’s state, we choose the parameters for x3≤x3,0 and for x3≥x3,1 independent of the x3-coordinate, according to αxL and βxL. Hence, if we let Green’s state (with a point source at xA in D) take the role of state A, then the condition for the validity of equation (66) is fulfilled. Moreover, Green’s functions are purely outward propagating at ∂D0,1 (because outside D no scattering occurs along the x3-coordinate). Hence, G+,±x,xA,ω=0 at ∂D0 and G−,±x,xA,ω=0 at ∂D1. We let a second Green’s state (with a point source at xB in D and the same parameters α and β as in state A, inside as well as outside D) take the role of state B. Hence, G+,±x,xB,ω=0 at ∂D0 and G−,±x,xB,ω=0 at ∂D1. With only outward propagating waves at ∂D0,1, the surface integral on the right-hand side of equation (70) vanishes. Hence, taking into account that B0±=B± (since Cj=0), equation (70) simplifies to
(79)∫D−BA+PB−−BA−PB++PA+BB−+PA−BB+dx=0.First, we consider sources emitting waves in the positivex3-direction in both Green’s states, hence BA+=δx−xA, BA−=0, PA±=G±,+x,xA,ω, BB+=δx−xB, BB−=0, and PB±=G±,+x,xB,ω. Substituting this into equation (79) yields
(80)G−,+xB,xA,ω=G−,+xA,xB,ω,see Figure 2(a). Next, we replace the source in state B by one emitting waves in the negative x3-direction, hence BB+=0, BB−=δx−xB, and PB±=G±,−x,xB,ω. This gives
(81)G+,+xB,xA,ω=G−,−xA,xB,ω,see Figure 2(b). By replacing also the source in state A by one emitting waves in the negative x3-direction, according to BA+=0, BA−=δx−xA, and PA±=G±,−x,xA,ω, we obtain
(82)G+,−xB,xA,ω=G+,−xA,xB,ω,see Figure 2(c). Finally, changing the source in state B back to the one emitting waves in the positive x3-direction yields
(83)G−,−xB,xA,ω=G+,+xA,xB,ω,see Figure 2(d).Figure 2
Visualisation of the source-receiver reciprocity relations for the field-normalised one-way Green’s functions, formulated by equations (80), (81), (82), and (83). The “rays” in this and subsequent figures are strong simplifications of the complete one-way wave fields, which include primary and multiple scattering.
(a)
(b)
(c)
(d)Source-receiver reciprocity relations similar to equations (80), (81), (82), and (83) were previously derived for flux-normalised one-way Green’s functions [17], except that two of those relations involve a change of sign when interchanging the source and the receiver. The absence of sign changes in equations (80), (81), (82), and (83) is due to the definition of B0± in equation (69). Moreover, unlike the flux-normalised reciprocity relations, the field-normalised source-receiver reciprocity relations of equations (80), (81), (82), and (83) have a very straightforward relation with the well-known source-receiver reciprocity relation for the full Green’s function. By separately summing the left- and right-hand sides of equations (80), (81), (82), and (83) and using equation (78), we simply obtain
(84)GxB,xA,ω=GxA,xB,ω.


### 4.3. Kirchhoff-Helmholtz Integrals for Forward Propagation

We derive Kirchhoff-Helmholtz integrals of the convolution type for field-normalised one-way wave fields. For stateB, we consider the decomposed actual field, with sources only outside D; hence, B0,B±=0 in D and PB±=P±x,ω. The parameters α and β are the actual parameters inside as well as outside D. For state A, we choose the Green’s state with a unit point source at xA in D. The parameters α and β in D are the same as those in state B, but for x3≤x3,0 and for x3≥x3,1, they are chosen independent of the x3-coordinate. Hence, the condition for the validity of equation (66) is fulfilled. First, we consider a source in state A which emits waves in the positive x3-direction, hence BA+=δx−xA, BA−=0, and PA±=G±,+x,xA,ω. Substituting all this into equation (70) (with B0,A±=BA±) gives
(85)P−xA,ω=∫∂D0,12iωβx∂3G+,+x,xA,ωP−x,ω+∂3G−,+x,xA,ωP+x,ωn3dxL.Next, we replace the source in stateA by one which emits waves in the negative x3-direction, hence BA+=0, BA−=δx−xA, and PA±=G±,−x,xA,ω. Equation (70) thus gives
(86)P+xA,ω=∫∂D0,12iωβx∂3G+,−x,xA,ωP−x,ω+∂3G−,−x,xA,ωP+x,ωn3dxL.Recall that∂D0,1 consists of ∂D0 (with n3=−1) and ∂D1 (with n3=+1), see Figure 1. Since G+,±x,xA,ω=0 at ∂D0 and G−,±x,xA,ω=0 at ∂D1 (because outside D no scattering occurs along the x3-coordinate in state A), the first term under the integral in equations (85) and (86) gives a contribution only at ∂D1 and the second term only at ∂D0.Hence,(87)P±xA,ω=∫∂D0−2iωβx∂3G−,∓x,xA,ωP+x,ωdxL+∫∂D12iωβx∂3G+,∓x,xA,ωP−x,ωdxL.Note that there is no contribution fromP−x,ω at ∂D0 nor from P+x,ω at ∂D1, see Figure 3.Figure 3
Visualisation of the different terms in the field-normalised one-way Kirchhoff-Helmholtz integral for forward propagation, formulated by equation (87). The solid Green’s functions contribute to P+xA,ω, the dashed Green’s functions to P−xA,ω.We conclude this section by considering a special case. Suppose the source of the actual field (stateB) is located at xB in the half-space x3<x3,0. Then, by taking x3,1→∞, the field P− at ∂D1 vanishes. This leaves the single-sided representation
(88)P±xA,xB,ω=∫∂D0−2iωβx∂3G−,∓x,xA,ωP+x,xB,ωdxL.Note that we included the source coordinate vectorxB in the argument list of P±xA,xB,ω. This representation is an extension of a previously derived result [43], in which the fields were decomposed at ∂D0 but not at xA. It describes forward propagation of the one-way field P+x,xB,ω from the surface ∂D0 to xA (with xA and xB defined at opposite sides of ∂D0). In the following two sections, we discuss representations for backward propagation of one-way wave fields.


### 4.4. Kirchhoff-Helmholtz Integrals for Backward Propagation (Double-Sided)

We derive Kirchhoff-Helmholtz integrals of the correlation type for field-normalised one-way wave fields. For stateB, we consider the decomposed actual field, with a point source at xB and source spectrum sω. The parameters α and β are the actual parameters inside as well as outside D. For state A, we choose the Green’s state with a unit point source at xA in D. The parameters α and β in D are the same as those in state B, but for x3≤x3,0 and for x3≥x3,1, they are chosen independent of the x3-coordinate. Hence, the condition for the validity of equation (66) is fulfilled. First, we consider sources emitting waves in the positive x3-direction in both states, hence BA+=δx−xA, BA−=0, PA±=G±,+x,xA,ω, BB+=δx−xBsω, BB−=0, and PB±=P±,+x,xB,ω. Substituting this into equation (71) (with B0,A±=BA± and B0,B±=BB±) gives
(89)P+,+xA,xB,ω+χxBG+,+xB,xA,ω∗sω=∫∂D0,1−2iωβx∂3G+,+x,xA,ω∗P+,+x,xB,ω+∂3G−,+x,xA,ω∗P−,+x,xB,ωn3dxL,where χ is the characteristic function of the domain D. It is defined as
(90)χxB=1,forxBinD,12,forxBon∂D0,1,0,forxBoutsideD.SinceG+,+x,xA,ω=0 at ∂D0 and G−,+x,xA,ω=0 at ∂D1 (because outside D no scattering occurs along the x3-coordinate in state A), the first term under the integral in equation (89) gives a contribution only at ∂D1 and the second term only at ∂D0. Hence,
(91)P+,+xA,xB,ω+χxBG+,+xB,xA,ω∗sω=∫∂D02iωβx∂3G−,+x,xA,ω∗P−,+x,xB,ωdxL−∫∂D12iωβx∂3G+,+x,xA,ω∗P+,+x,xB,ωdxL.Next, we replace the source in stateB by one emitting waves in the negative x3-direction, hence BB+=0,BB−=δx−xBsω and PB±=P±,−x,xB,ω. This gives
(92)P+,−xA,xB,ω+χxBG−,+xB,xA,ω∗sω=∫∂D02iωβx∂3G−,+x,xA,ω∗P−,−x,xB,ωdxL−∫∂D12iωβx∂3G+,+x,xA,ω∗P+,−x,xB,ωdxL.By replacing also the source in stateA by one emitting waves in the negative x3-direction, according to BA+=0, BA−=δx−xA, and PA±=G±,−x,xA,ω, we obtain
(93)P−,−xA,xB,ω+χxBG−,−xB,xA,ω∗sω=∫∂D02iωβx∂3G−,−x,xA,ω∗P−,−x,xB,ωdxL−∫∂D12iωβx∂3G+,−x,xA,ω∗P+,−x,xB,ωdxL.Finally, changing the source in stateB back to the one emitting waves in the positive x3-direction yields
(94)P−,+xA,xB,ω+χxBG+,−xB,xA,ω∗sω=∫∂D02iωβx∂3G−,−x,xA,ω∗P−,+x,xB,ωdxL−∫∂D12iωβx∂3G+,−x,xA,ω∗P+,+x,xB,ωdxL.Equation (93) is an extension of a previously derived result [44], in which the fields were decomposed at ∂D0,1 but not at xA and xB. Equations (91), (92), and (94) are further variations. Equation (94) is visualised in Figure 4. Together, these equations describe backward propagation of the one-way wave fields P−,±x,xB,ω from ∂D0 and P+,±x,xB,ω from ∂D1 to xA. Except for some special cases, the integrals along ∂D1 do not vanish by taking x3,1→∞. Hence, unlike the forward propagation representation (87), the double-sided backward propagation representations (91), (92), (93), and (94) in general do not simplify to single-sided representations. In the next section, we discuss an alternative method to derive single-sided representations for backward propagation.Figure 4
Visualisation of the different terms in the field-normalised one-way Kirchhoff-Helmholtz integral for backward propagation, formulated by equation (94).We conclude this section by considering a special case. Suppose that in stateB the parameters α and β are the same as in state A not only in D but also outside D. Then, P±,±x,xB,ω=G±,±x,xB,ωsω for all x. Substituting this into representations (91), (92), (93), and (94), summing the left- and right-hand sides of these representations separately and dividing both sides by sω, using equations (78) and (84) and assuming that xB is located in D, we obtain
(95)GhxA,xB,ω=∫∂D02iωβx∂3G−x,xA,ω∗G−x,xB,ωdxL−∫∂D12iωβx∂3G+x,xA,ω∗G+x,xB,ωdxL,where the so-called homogeneous Green’s function GhxA,xB,ω is defined as
(96)GhxA,xB,ω=GxA,xB,ω+G∗xA,xB,ω=2RGxA,xB,ω,(withR denoting the real part) and where G±x,xA,ω=G±,+x,xA,ω+G±,−x,xA,ω (and a similar expression for G±x,xB,ω). Equation (95) is akin to the well-known representation for the homogeneous Green’s function [45, 46], but with decomposed Green’s functions under the integrals. The simple relation between representations (91), (92), (93), and (94) on the one hand and the homogeneous Green’s function representation (95) on the other hand is a consequence of the field-normalised decomposition, introduced in Section 3.3.


### 4.5. Kirchhoff-Helmholtz Integrals for Backward Propagation (Single-Sided)

The complex-conjugated Green’s functions∂3G±,±x,xA,ω∗ under the integrals in equations (91), (92), (93), and (94) can be seen as focusing functions, which focus the wave fields P±,±x,xB,ω onto a focal point xA. However, this focusing process requires that these wave fields are available at two boundaries ∂D0 and ∂D1, enclosing the focal point xA. Here, we discuss single-sided field-normalised focusing functions f1±x,xA,ω and we use these in modifications of reciprocity theorems (70) and (71) to derive single-sided Kirchhoff-Helmholtz integrals for backward propagation.We start by defining a new domainDA, enclosed by two surfaces ∂D0 and ∂DA perpendicular to the x3-axis at x3=x3,0 and x3=x3,A, respectively, with x3,A>x3,0, see Figure 5. Hence, ∂DA is chosen such that it contains the focal point xA. The two surfaces ∂D0 and ∂DA are together denoted by ∂D0,A. The focusing functions f1±x,xA,ω, which will play the role of state A in the reciprocity theorems, obey the one-way wave equations (64) and (65) (but without the source terms S±), with parameters α and β in DA equal to those in the actual state B, and independent of the x3-coordinate for x3≤x3,0 and for x3≥x3,A. Hence, the condition for the validity of equation (66) is fulfilled. Analogous to equation (56), the field-normalised focusing functions f1±x,xA,ω are related to the full focusing function f1x,xA,ω, according to
(97)f1x,xA,ω=f1+x,xA,ω+f1−x,xA,ω.Figure 5
Configuration for the derivation of the single-sided Kirchhoff-Helmholtz integrals for backward propagation.The focusing functionf1+x,xA,ω is incident to the domain DA from the half-space x3<x3,0 (see Figure 5). It propagates and scatters in the inhomogeneous domain DA, focuses at xA on surface ∂DA, and continues as f1+x,xA,ω in the half-space x3>x3,A. The back-scattered field leaves DA via surface ∂D0 and continues as f1−x,xA,ω in the half-space x3<x3,0. The focusing conditions at the focal plane ∂DA are [18]
(98)∂3f1+x,xA,ωx3=x3,A=12iωβxAδxL−xL,A,(99)∂3f1−x,xA,ωx3=x3,A=0.Here,xL,A denotes the lateral coordinates of xA. The operators ∂3 and the factor 1/2iωβxA are not necessary to define the focusing conditions but are chosen for later convenience. To avoid instability, evanescent waves are excluded from the focusing functions. This implies that the delta function in equation (98) should be interpreted as a spatially band-limited delta function. Note that the sifting property of the delta function, hxL,A=∫SδxL−xL,AhxLdxL, remains valid for a spatially band-limited delta function, assuming hxL is also spatially band-limited.We now derive single-sided Kirchhoff-Helmholtz integrals for backward propagation. We consider the reciprocity theorems for field-normalised one-way wave fields (equations (70) and (71)), with D and ∂D0,1 replaced by DA and ∂D0,A, respectively. For state A, we consider the focusing functions discussed above; hence, BA+x,ω=BA−x,ω=0 and PA±x,ω=f1±x,xA,ω. For state B, we consider the decomposed actual field, with a point source at xB in the half-space x3>x3,0 and source spectrum sω. The parameters α and β in state B are the actual parameters inside as well as outside ∂D0,A. First, we consider a source in state B which emits waves in the positive x3-direction, hence BB+x,ω=δx−xBsω, BB−x,ω=0, and PB±x,ω=P±,+x,xB,ω. Substituting all this into equations (70) and (71) (with B0±=B±), using equations (98) and (99) in the integrals along ∂DA, gives
(100)P−,+xA,xB,ω+χAxBf1−xB,xA,ωsω=∫∂D02iωβx∂3f1+x,xA,ωP−,+x,xB,ω+∂3f1−x,xA,ωP+,+x,xB,ωdxL,(101)P+,+xA,xB,ω−χAxBf1+xB,xA,ω∗sω=∫∂D0−2iωβx∂3f1+x,xA,ω∗P+,+x,xB,ω+∂3f1−x,xA,ω∗P−,+x,xB,ωdxL,where χA is the characteristic function of the domain DA. It is defined by equation (90), with D and ∂D0,1 replaced by DA and ∂D0,A, respectively. Next, we replace the source in state B by one which emits waves in the negative x3-direction, hence BB+x,ω=0, BB−x,ω=δx−xBsω, and PB±x,ω=P±,−x,xB,ω. This gives
(102)P−,−xA,xB,ω+χAxBf1+xB,xA,ωsω=∫∂D02iωβx∂3f1+x,xA,ωP−,−x,xB,ω+∂3f1−x,xA,ωP+,−x,xB,ωdxL,(103)P+,−xA,xB,ω−χAxBf1−xB,xA,ω∗sω=∫∂D0−2iωβx∂3f1+x,xA,ω∗P+,−x,xB,ω+∂3f1−x,xA,ω∗P−,−x,xB,ωdxL.Equations (100), (101), (102), and (103) are single-sided representations for backward propagation of the one-way wave fields P±,±x,xB,ω from ∂D0 to xA. Similar results have been previously obtained [47, 48], but without decomposition at xB. An advantage of these equations over equations (91), (92), (93), and (94) is that the backward propagated fields P±,±xA,xB,ω are expressed entirely in terms of integrals along the surface ∂D0.Single-sided representations containing the field-normalised focusing functionsf1±x,xA,ω find applications, for example, in reflection imaging methods, which account for multiple scattering. In these methods, the focusing functions are retrieved from the reflection response at the surface ∂D0, using the Marchenko method [18, 49–51].We conclude this section by considering a special case. Suppose that in stateB the parameters α and β are the same as in the Green’s state. Then, P±,±x,xB,ω=G±,±x,xB,ωsω for all x. Moreover, P+,±x,xB,ω=0 for x at ∂D0. Substituting this into representations (100), (101), (102), and (103), summing the left- and right-hand sides of these representations separately, dividing both sides by sω and using equation (97), we obtain
(104)GxA,xB,ω+χAxB2iIf1xB,xA,ω=∫∂D02iωβx∂3f1+x,xA,ω−f1−x,xAω∗G−x,xB,ωdxL,(withI denoting the imaginary part), where G−x,xB,ω=G−,+x,xB,ω+G−,−x,xB,ω. Taking the real part of both sides gives
(105)GhxA,xB,ω=R∫∂D04iωβx∂3f1+x,xA,ω−f1−x,xA,ω∗G−x,xB,ωdx,where GhxA,xB,ω is the homogeneous Green’s function, defined in equation (96). Unlike in equation (95), here the homogeneous Green’s function is represented by a single integral along the surface ∂D0, containing field normalised one-way focusing and Green’s functions.


## 4.1. Green’s Functions

Representation theorems are obtained by substituting Green’s functions in reciprocity theorems. Our aim is to introduce one-way Green’s functions, to be used in the reciprocity theorems for field-normalised one-way wave fields (equations (70) and (71)). First, we introduce the full Green’s function Gx,xA,ω as a solution of the unified wave equation (3) for a unit monopole point source at xA, with Bx,ω=δx−xA and Cjx,ω=0. Hence,
(72)β∂j1β∂jG+k2G=iωβδx−xA.As boundary condition, we impose the radiation condition (i.e., outward propagating waves at infinity). Next, we introduce one-way Green’s function as solutions of the coupled one-way equations (64) and (65) for a unit monopole point source at xA. Hence, we choose again Bx,ω=δx−xA and Cjx,ω=0. Using equations (69) and (7), we define decomposed sources as B0±=B±=B=±2L2S±, with L2 defined in equation (55), or
(73)S±x,ω=±12L2−1B±x,ω=±12L2−1Bx,ω=±12L2−1δx−xA.We consider two sets of one-way Green’s functions. For the first set, we choose a point sourceS+x,ω=1/2L2−1B+x,ω, with B+x,ω=δx−xA, which emits waves from xA in the positive x3-direction, and we set S−x,ω equal to zero. Hence, for this first set, one-way equations (64) and (65) become
(74)∂3G+,+=+iH1G+,+−12L2−1∂3L2G+,+−G−,++12L2−1δx−xA,(75)∂3G−,+=−iH1G−,++12L2−1∂3L2G+,+−G−,+.Here,G±,+ stands for G±,+x,xA,ω. The second superscript (+) indicates that the source at xA emits waves in the positive x3-direction. The first superscript (±) denotes the propagation direction at x. For the second set of one-way Green’s functions, we choose a point source S−x,ω=−1/2L2−1B−x,ω, with B−x,ω=δx−xA, which emits waves from xA in the negative x3-direction, and we set S+x,ω equal to zero. Hence, for this second set, one-way equations (64) and (65) become
(76)∂3G+,−=+iH1G+,−−12L2−1∂3L2G+,−−G−,−,(77)∂3G−,−=−iH1G−,−+12L2−1∂3L2G+,−−G−,−−12L2−1δx−xA.Here,G±,− stands for G±,−x,xA,ω, with the second superscript (−) indicating that the source at xA emits waves in the negative x3-direction. Like for the full Green’s function Gx,xA,ω, we impose radiation conditions for both sets of one-way Green’s functions.To find a relation between the full Green’s function and the one-way Green’s functions, we evaluateβ∂31/β∂3G+,++G−,++G+,−+G−,− using equations (74), (75), (76), (77), (25), (53), and (55). This gives equation (72), with G replaced by G+,++G−,++G+,−+G−,−. Since the full Green’s function and the one-way Green’s functions obey the same radiation conditions, we thus find
(78)G=G+,++G−,++G+,−+G−,−.This very simple relation is a consequence of the field-normalised decomposition, introduced in Section3.3.


## 4.2. Source-Receiver Reciprocity

We derive source-receiver reciprocity relations for the field-normalised one-way Green’s functions introduced in the previous section. To this end, we make use of the reciprocity theorem of the convolution type for field-normalised one-way wave fields (equation (70)). This theorem was derived for the configuration of Figure 1, assuming that in domain D, the parameters α and β are the same in the two states (see Section 2.4). Outside D, these parameters may be different in the two states. For the Green’s state, we choose the parameters for x3≤x3,0 and for x3≥x3,1 independent of the x3-coordinate, according to αxL and βxL. Hence, if we let Green’s state (with a point source at xA in D) take the role of state A, then the condition for the validity of equation (66) is fulfilled. Moreover, Green’s functions are purely outward propagating at ∂D0,1 (because outside D no scattering occurs along the x3-coordinate). Hence, G+,±x,xA,ω=0 at ∂D0 and G−,±x,xA,ω=0 at ∂D1. We let a second Green’s state (with a point source at xB in D and the same parameters α and β as in state A, inside as well as outside D) take the role of state B. Hence, G+,±x,xB,ω=0 at ∂D0 and G−,±x,xB,ω=0 at ∂D1. With only outward propagating waves at ∂D0,1, the surface integral on the right-hand side of equation (70) vanishes. Hence, taking into account that B0±=B± (since Cj=0), equation (70) simplifies to
(79)∫D−BA+PB−−BA−PB++PA+BB−+PA−BB+dx=0.First, we consider sources emitting waves in the positivex3-direction in both Green’s states, hence BA+=δx−xA, BA−=0, PA±=G±,+x,xA,ω, BB+=δx−xB, BB−=0, and PB±=G±,+x,xB,ω. Substituting this into equation (79) yields
(80)G−,+xB,xA,ω=G−,+xA,xB,ω,see Figure 2(a). Next, we replace the source in state B by one emitting waves in the negative x3-direction, hence BB+=0, BB−=δx−xB, and PB±=G±,−x,xB,ω. This gives
(81)G+,+xB,xA,ω=G−,−xA,xB,ω,see Figure 2(b). By replacing also the source in state A by one emitting waves in the negative x3-direction, according to BA+=0, BA−=δx−xA, and PA±=G±,−x,xA,ω, we obtain
(82)G+,−xB,xA,ω=G+,−xA,xB,ω,see Figure 2(c). Finally, changing the source in state B back to the one emitting waves in the positive x3-direction yields
(83)G−,−xB,xA,ω=G+,+xA,xB,ω,see Figure 2(d).Figure 2
Visualisation of the source-receiver reciprocity relations for the field-normalised one-way Green’s functions, formulated by equations (80), (81), (82), and (83). The “rays” in this and subsequent figures are strong simplifications of the complete one-way wave fields, which include primary and multiple scattering.
(a)
(b)
(c)
(d)Source-receiver reciprocity relations similar to equations (80), (81), (82), and (83) were previously derived for flux-normalised one-way Green’s functions [17], except that two of those relations involve a change of sign when interchanging the source and the receiver. The absence of sign changes in equations (80), (81), (82), and (83) is due to the definition of B0± in equation (69). Moreover, unlike the flux-normalised reciprocity relations, the field-normalised source-receiver reciprocity relations of equations (80), (81), (82), and (83) have a very straightforward relation with the well-known source-receiver reciprocity relation for the full Green’s function. By separately summing the left- and right-hand sides of equations (80), (81), (82), and (83) and using equation (78), we simply obtain
(84)GxB,xA,ω=GxA,xB,ω.


## 4.3. Kirchhoff-Helmholtz Integrals for Forward Propagation

We derive Kirchhoff-Helmholtz integrals of the convolution type for field-normalised one-way wave fields. For stateB, we consider the decomposed actual field, with sources only outside D; hence, B0,B±=0 in D and PB±=P±x,ω. The parameters α and β are the actual parameters inside as well as outside D. For state A, we choose the Green’s state with a unit point source at xA in D. The parameters α and β in D are the same as those in state B, but for x3≤x3,0 and for x3≥x3,1, they are chosen independent of the x3-coordinate. Hence, the condition for the validity of equation (66) is fulfilled. First, we consider a source in state A which emits waves in the positive x3-direction, hence BA+=δx−xA, BA−=0, and PA±=G±,+x,xA,ω. Substituting all this into equation (70) (with B0,A±=BA±) gives
(85)P−xA,ω=∫∂D0,12iωβx∂3G+,+x,xA,ωP−x,ω+∂3G−,+x,xA,ωP+x,ωn3dxL.Next, we replace the source in stateA by one which emits waves in the negative x3-direction, hence BA+=0, BA−=δx−xA, and PA±=G±,−x,xA,ω. Equation (70) thus gives
(86)P+xA,ω=∫∂D0,12iωβx∂3G+,−x,xA,ωP−x,ω+∂3G−,−x,xA,ωP+x,ωn3dxL.Recall that∂D0,1 consists of ∂D0 (with n3=−1) and ∂D1 (with n3=+1), see Figure 1. Since G+,±x,xA,ω=0 at ∂D0 and G−,±x,xA,ω=0 at ∂D1 (because outside D no scattering occurs along the x3-coordinate in state A), the first term under the integral in equations (85) and (86) gives a contribution only at ∂D1 and the second term only at ∂D0.Hence,(87)P±xA,ω=∫∂D0−2iωβx∂3G−,∓x,xA,ωP+x,ωdxL+∫∂D12iωβx∂3G+,∓x,xA,ωP−x,ωdxL.Note that there is no contribution fromP−x,ω at ∂D0 nor from P+x,ω at ∂D1, see Figure 3.Figure 3
Visualisation of the different terms in the field-normalised one-way Kirchhoff-Helmholtz integral for forward propagation, formulated by equation (87). The solid Green’s functions contribute to P+xA,ω, the dashed Green’s functions to P−xA,ω.We conclude this section by considering a special case. Suppose the source of the actual field (stateB) is located at xB in the half-space x3<x3,0. Then, by taking x3,1→∞, the field P− at ∂D1 vanishes. This leaves the single-sided representation
(88)P±xA,xB,ω=∫∂D0−2iωβx∂3G−,∓x,xA,ωP+x,xB,ωdxL.Note that we included the source coordinate vectorxB in the argument list of P±xA,xB,ω. This representation is an extension of a previously derived result [43], in which the fields were decomposed at ∂D0 but not at xA. It describes forward propagation of the one-way field P+x,xB,ω from the surface ∂D0 to xA (with xA and xB defined at opposite sides of ∂D0). In the following two sections, we discuss representations for backward propagation of one-way wave fields.


## 4.4. Kirchhoff-Helmholtz Integrals for Backward Propagation (Double-Sided)

We derive Kirchhoff-Helmholtz integrals of the correlation type for field-normalised one-way wave fields. For stateB, we consider the decomposed actual field, with a point source at xB and source spectrum sω. The parameters α and β are the actual parameters inside as well as outside D. For state A, we choose the Green’s state with a unit point source at xA in D. The parameters α and β in D are the same as those in state B, but for x3≤x3,0 and for x3≥x3,1, they are chosen independent of the x3-coordinate. Hence, the condition for the validity of equation (66) is fulfilled. First, we consider sources emitting waves in the positive x3-direction in both states, hence BA+=δx−xA, BA−=0, PA±=G±,+x,xA,ω, BB+=δx−xBsω, BB−=0, and PB±=P±,+x,xB,ω. Substituting this into equation (71) (with B0,A±=BA± and B0,B±=BB±) gives
(89)P+,+xA,xB,ω+χxBG+,+xB,xA,ω∗sω=∫∂D0,1−2iωβx∂3G+,+x,xA,ω∗P+,+x,xB,ω+∂3G−,+x,xA,ω∗P−,+x,xB,ωn3dxL,where χ is the characteristic function of the domain D. It is defined as
(90)χxB=1,forxBinD,12,forxBon∂D0,1,0,forxBoutsideD.SinceG+,+x,xA,ω=0 at ∂D0 and G−,+x,xA,ω=0 at ∂D1 (because outside D no scattering occurs along the x3-coordinate in state A), the first term under the integral in equation (89) gives a contribution only at ∂D1 and the second term only at ∂D0. Hence,
(91)P+,+xA,xB,ω+χxBG+,+xB,xA,ω∗sω=∫∂D02iωβx∂3G−,+x,xA,ω∗P−,+x,xB,ωdxL−∫∂D12iωβx∂3G+,+x,xA,ω∗P+,+x,xB,ωdxL.Next, we replace the source in stateB by one emitting waves in the negative x3-direction, hence BB+=0,BB−=δx−xBsω and PB±=P±,−x,xB,ω. This gives
(92)P+,−xA,xB,ω+χxBG−,+xB,xA,ω∗sω=∫∂D02iωβx∂3G−,+x,xA,ω∗P−,−x,xB,ωdxL−∫∂D12iωβx∂3G+,+x,xA,ω∗P+,−x,xB,ωdxL.By replacing also the source in stateA by one emitting waves in the negative x3-direction, according to BA+=0, BA−=δx−xA, and PA±=G±,−x,xA,ω, we obtain
(93)P−,−xA,xB,ω+χxBG−,−xB,xA,ω∗sω=∫∂D02iωβx∂3G−,−x,xA,ω∗P−,−x,xB,ωdxL−∫∂D12iωβx∂3G+,−x,xA,ω∗P+,−x,xB,ωdxL.Finally, changing the source in stateB back to the one emitting waves in the positive x3-direction yields
(94)P−,+xA,xB,ω+χxBG+,−xB,xA,ω∗sω=∫∂D02iωβx∂3G−,−x,xA,ω∗P−,+x,xB,ωdxL−∫∂D12iωβx∂3G+,−x,xA,ω∗P+,+x,xB,ωdxL.Equation (93) is an extension of a previously derived result [44], in which the fields were decomposed at ∂D0,1 but not at xA and xB. Equations (91), (92), and (94) are further variations. Equation (94) is visualised in Figure 4. Together, these equations describe backward propagation of the one-way wave fields P−,±x,xB,ω from ∂D0 and P+,±x,xB,ω from ∂D1 to xA. Except for some special cases, the integrals along ∂D1 do not vanish by taking x3,1→∞. Hence, unlike the forward propagation representation (87), the double-sided backward propagation representations (91), (92), (93), and (94) in general do not simplify to single-sided representations. In the next section, we discuss an alternative method to derive single-sided representations for backward propagation.Figure 4
Visualisation of the different terms in the field-normalised one-way Kirchhoff-Helmholtz integral for backward propagation, formulated by equation (94).We conclude this section by considering a special case. Suppose that in stateB the parameters α and β are the same as in state A not only in D but also outside D. Then, P±,±x,xB,ω=G±,±x,xB,ωsω for all x. Substituting this into representations (91), (92), (93), and (94), summing the left- and right-hand sides of these representations separately and dividing both sides by sω, using equations (78) and (84) and assuming that xB is located in D, we obtain
(95)GhxA,xB,ω=∫∂D02iωβx∂3G−x,xA,ω∗G−x,xB,ωdxL−∫∂D12iωβx∂3G+x,xA,ω∗G+x,xB,ωdxL,where the so-called homogeneous Green’s function GhxA,xB,ω is defined as
(96)GhxA,xB,ω=GxA,xB,ω+G∗xA,xB,ω=2RGxA,xB,ω,(withR denoting the real part) and where G±x,xA,ω=G±,+x,xA,ω+G±,−x,xA,ω (and a similar expression for G±x,xB,ω). Equation (95) is akin to the well-known representation for the homogeneous Green’s function [45, 46], but with decomposed Green’s functions under the integrals. The simple relation between representations (91), (92), (93), and (94) on the one hand and the homogeneous Green’s function representation (95) on the other hand is a consequence of the field-normalised decomposition, introduced in Section 3.3.


## 4.5. Kirchhoff-Helmholtz Integrals for Backward Propagation (Single-Sided)

The complex-conjugated Green’s functions∂3G±,±x,xA,ω∗ under the integrals in equations (91), (92), (93), and (94) can be seen as focusing functions, which focus the wave fields P±,±x,xB,ω onto a focal point xA. However, this focusing process requires that these wave fields are available at two boundaries ∂D0 and ∂D1, enclosing the focal point xA. Here, we discuss single-sided field-normalised focusing functions f1±x,xA,ω and we use these in modifications of reciprocity theorems (70) and (71) to derive single-sided Kirchhoff-Helmholtz integrals for backward propagation.We start by defining a new domainDA, enclosed by two surfaces ∂D0 and ∂DA perpendicular to the x3-axis at x3=x3,0 and x3=x3,A, respectively, with x3,A>x3,0, see Figure 5. Hence, ∂DA is chosen such that it contains the focal point xA. The two surfaces ∂D0 and ∂DA are together denoted by ∂D0,A. The focusing functions f1±x,xA,ω, which will play the role of state A in the reciprocity theorems, obey the one-way wave equations (64) and (65) (but without the source terms S±), with parameters α and β in DA equal to those in the actual state B, and independent of the x3-coordinate for x3≤x3,0 and for x3≥x3,A. Hence, the condition for the validity of equation (66) is fulfilled. Analogous to equation (56), the field-normalised focusing functions f1±x,xA,ω are related to the full focusing function f1x,xA,ω, according to
(97)f1x,xA,ω=f1+x,xA,ω+f1−x,xA,ω.Figure 5
Configuration for the derivation of the single-sided Kirchhoff-Helmholtz integrals for backward propagation.The focusing functionf1+x,xA,ω is incident to the domain DA from the half-space x3<x3,0 (see Figure 5). It propagates and scatters in the inhomogeneous domain DA, focuses at xA on surface ∂DA, and continues as f1+x,xA,ω in the half-space x3>x3,A. The back-scattered field leaves DA via surface ∂D0 and continues as f1−x,xA,ω in the half-space x3<x3,0. The focusing conditions at the focal plane ∂DA are [18]
(98)∂3f1+x,xA,ωx3=x3,A=12iωβxAδxL−xL,A,(99)∂3f1−x,xA,ωx3=x3,A=0.Here,xL,A denotes the lateral coordinates of xA. The operators ∂3 and the factor 1/2iωβxA are not necessary to define the focusing conditions but are chosen for later convenience. To avoid instability, evanescent waves are excluded from the focusing functions. This implies that the delta function in equation (98) should be interpreted as a spatially band-limited delta function. Note that the sifting property of the delta function, hxL,A=∫SδxL−xL,AhxLdxL, remains valid for a spatially band-limited delta function, assuming hxL is also spatially band-limited.We now derive single-sided Kirchhoff-Helmholtz integrals for backward propagation. We consider the reciprocity theorems for field-normalised one-way wave fields (equations (70) and (71)), with D and ∂D0,1 replaced by DA and ∂D0,A, respectively. For state A, we consider the focusing functions discussed above; hence, BA+x,ω=BA−x,ω=0 and PA±x,ω=f1±x,xA,ω. For state B, we consider the decomposed actual field, with a point source at xB in the half-space x3>x3,0 and source spectrum sω. The parameters α and β in state B are the actual parameters inside as well as outside ∂D0,A. First, we consider a source in state B which emits waves in the positive x3-direction, hence BB+x,ω=δx−xBsω, BB−x,ω=0, and PB±x,ω=P±,+x,xB,ω. Substituting all this into equations (70) and (71) (with B0±=B±), using equations (98) and (99) in the integrals along ∂DA, gives
(100)P−,+xA,xB,ω+χAxBf1−xB,xA,ωsω=∫∂D02iωβx∂3f1+x,xA,ωP−,+x,xB,ω+∂3f1−x,xA,ωP+,+x,xB,ωdxL,(101)P+,+xA,xB,ω−χAxBf1+xB,xA,ω∗sω=∫∂D0−2iωβx∂3f1+x,xA,ω∗P+,+x,xB,ω+∂3f1−x,xA,ω∗P−,+x,xB,ωdxL,where χA is the characteristic function of the domain DA. It is defined by equation (90), with D and ∂D0,1 replaced by DA and ∂D0,A, respectively. Next, we replace the source in state B by one which emits waves in the negative x3-direction, hence BB+x,ω=0, BB−x,ω=δx−xBsω, and PB±x,ω=P±,−x,xB,ω. This gives
(102)P−,−xA,xB,ω+χAxBf1+xB,xA,ωsω=∫∂D02iωβx∂3f1+x,xA,ωP−,−x,xB,ω+∂3f1−x,xA,ωP+,−x,xB,ωdxL,(103)P+,−xA,xB,ω−χAxBf1−xB,xA,ω∗sω=∫∂D0−2iωβx∂3f1+x,xA,ω∗P+,−x,xB,ω+∂3f1−x,xA,ω∗P−,−x,xB,ωdxL.Equations (100), (101), (102), and (103) are single-sided representations for backward propagation of the one-way wave fields P±,±x,xB,ω from ∂D0 to xA. Similar results have been previously obtained [47, 48], but without decomposition at xB. An advantage of these equations over equations (91), (92), (93), and (94) is that the backward propagated fields P±,±xA,xB,ω are expressed entirely in terms of integrals along the surface ∂D0.Single-sided representations containing the field-normalised focusing functionsf1±x,xA,ω find applications, for example, in reflection imaging methods, which account for multiple scattering. In these methods, the focusing functions are retrieved from the reflection response at the surface ∂D0, using the Marchenko method [18, 49–51].We conclude this section by considering a special case. Suppose that in stateB the parameters α and β are the same as in the Green’s state. Then, P±,±x,xB,ω=G±,±x,xB,ωsω for all x. Moreover, P+,±x,xB,ω=0 for x at ∂D0. Substituting this into representations (100), (101), (102), and (103), summing the left- and right-hand sides of these representations separately, dividing both sides by sω and using equation (97), we obtain
(104)GxA,xB,ω+χAxB2iIf1xB,xA,ω=∫∂D02iωβx∂3f1+x,xA,ω−f1−x,xAω∗G−x,xB,ωdxL,(withI denoting the imaginary part), where G−x,xB,ω=G−,+x,xB,ω+G−,−x,xB,ω. Taking the real part of both sides gives
(105)GhxA,xB,ω=R∫∂D04iωβx∂3f1+x,xA,ω−f1−x,xA,ω∗G−x,xB,ωdx,where GhxA,xB,ω is the homogeneous Green’s function, defined in equation (96). Unlike in equation (95), here the homogeneous Green’s function is represented by a single integral along the surface ∂D0, containing field normalised one-way focusing and Green’s functions.


## 5. Conclusions

We have considered flux-normalised and field-normalised decomposition of scalar wave fields into coupled one-way wave fields. The operators for field-normalised decomposition exhibit less symmetry than those for flux-normalised decomposition. Nevertheless, we have shown that reciprocity theorems can be derived for field-normalised one-way wave fields in a similar way as those for flux-normalised one-way wave fields. An additional condition for the reciprocity theorems for field-normalised one-way wave fields is that in one of the states the derivatives in thex3-direction of the parameters α and β vanish at the boundary of the considered domain. This condition is easily fulfilled when one of the states is a Green’s function or a focusing function, for which the parameters α and β can be freely chosen at and outside the boundary of the domain.We have used the reciprocity theorems for field-normalised one-way wave fields as a starting point for deriving representation theorems for field-normalised one-way wave fields in a systematic way. We obtained representations for forward and for backward propagations of one-way wave fields. These representations account for multiple scattering. Whereas the Kirchhoff-Helmholtz integrals for forward propagation can be easily transformed into single-sided representations, this transformation is less straightforward for the Kirchhoff-Helmholtz integrals for backward propagation. By replacing the Green’s functions by focusing functions, we obtained single-sided representations for backward propagation of field-normalised one-way wave fields. These representations are particularly useful to retrieve wave fields in the interior of a domain in situations where measurements can be carried out only at a single surface. An important application is reflection imaging, accounting for multiple scattering.



---

*Source: 9540135-2020-01-13.xml*