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

Reciprocity and Representation Theorems for Flux- and Field-Normalised Decomposed

Wave Fields

Wapenaar, Kees DOI 10.1155/2020/9540135 Publication date 2020 Document Version Final published version Published in

Advances in Mathematical Physics

Citation (APA)

Wapenaar, K. (2020). Reciprocity and Representation Theorems for Flux- and Field-Normalised Decomposed Wave Fields. Advances in Mathematical Physics, 2020, 1-15. [9540135].

https://doi.org/10.1155/2020/9540135 Important note

To cite this publication, please use the final published version (if applicable). Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Research Article

Reciprocity and Representation Theorems for Flux- and

Field-Normalised Decomposed Wave Fields

Kees Wapenaar

Department of Geoscience and Engineering, Delft University of Technology, 2600 GA Delft, Netherlands Correspondence should be addressed to Kees Wapenaar; c.p.a.wapenaar@tudelft.nl

Received 3 October 2019; Accepted 30 November 2019; Published 13 January 2020 Academic Editor: Remi Léandre

Copyright © 2020 Kees Wapenaar. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 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 wavefields, 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.

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 direc-tion in an acoustic wave guide, bounded by the water surface and the ocean bottom. Similarly, in communication engi-neering, microwaves or optical waves propagate as beams through electromagnetic or optical wave guides. Wave prop-agation in preferred directions is not restricted to wave guides. For example, in geophysical reflection imaging appli-cations, 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 elec-tromagnetic 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 represen-tation theorems.

Decomposition of a wavefield into one-way wave fields is not unique. In particular, the amplitudes of the one-way wavefields 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 aflux-normalised one-way wave field is by def-inition the power-flux density (or, for quantum-mechanical waves, the probability-flux density) in the direction of prefer-ence. 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

Volume 2020, Article ID 9540135, 15 pages https://doi.org/10.1155/2020/9540135

symmetry than the operators for field-normalised decom-position. 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].

Thefirst aim of this paper is to discuss flux-normalised versusfield-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 sim-ilar way as those for flux-normalised one-way wave fields, even though the operators forfield-normalised decomposi-tion 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 propaga-tion, 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 straight-forward summation of the one-way wave fields gives the full wave field.

We restrict the discussion to scalar wavefields. In Sec-tion 2, we formulate a unified scalar wave equaSec-tion for acoustic waves, horizontally polarised shear waves, trans-verse electric and transtrans-verse 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 normal-isation and field normalisation and derive reciprocity the-orems for one-way wavefields, using both normalisations. In Section 4, we extensively discuss representation theorems for field-normalised one-way wave fields and indicate applica-tions. 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:

−iωαP + ∂jQj= B, ð1Þ

−iωβQj+ ∂jP = Cj: ð2Þ

Here, i is the imaginary unit and ω the angular fre-quency (in this paper, we consider positive frequencies only). Operator∂_j stands for the spatial differential opera-tor ∂/∂x_j, and Einstein’s summation convention applies to repeated subscripts. Pðx, ωÞ and Qjðx, ωÞ are space- and

frequency-dependent wave field quantities, αðxÞ and βðxÞ

are real-valued space-dependent parameters, and Bðx, ωÞ andC_jðx, ωÞ are space- and frequency-dependent source dis-tributions. Parametersα and β are both assumed to be posi-tive; hence, metamaterials are not considered in this paper. All quantities are specified in Table 1 for different wave phe-nomena and are discussed in more detail below. As indicated in thefirst column of Table 1, we consider 3D and 2D wave problems. For the 3D situation, x = ðx1, x2, x3Þ is the 3D

Car-tesian 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.

The unified boundary conditions at an interface between two media with different parameters state that P and n_jQ_jare continuous over the interface. Here,n_jrepresents the compo-nents 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 Table 1 in more detail. The quantities in row 1, associated to 3D acoustic wave propaga-tion in a losslessfluid medium, are acoustic pressure pðx, ωÞ, particle velocity vjðx, ωÞ, compressibility κðxÞ, mass density

ρðxÞ, volume-injection rate density qðx, ωÞ, and external

force density f_jðx, ωÞ. For 2D horizontally polarised shear waves in a lossless solid medium, we have in row 2 horizon-tal particle velocityv₂ðx, ωÞ, shear stress τ_2jðx, ωÞ, mass den-sity ρðxÞ, shear modulus μðxÞ, external horizontal force densityf₂ðx, ωÞ, and external shear deformation rate density

h_2jðx, ωÞ. Rows 3 and 4 contain the quantities for 2D electro-magnetic wave propagation, with TE standing for transverse electric and TM for transverse magnetic. The quantities are electric field strength E_kðx, ωÞ, magnetic field strength

H_kðx, ωÞ, permittivity εðxÞ, permeability μðxÞ, external elec-tric current density Je_kðx, ωÞ, and external magnetic current density Jm_kðx, ωÞ. Furthermore, ϵijk is the alternating tensor

(or Levi-Civita tensor), with ϵ₁₂₃= ϵ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 V ðxÞ, particle massm, and ℏ = h/2π, with h Planck’s constant.

By eliminatingQ_jfrom equations (1) and (2), we obtain the unified scalar wave equation

β∂j _β1∂jP



+ k2P = β∂j _β1Cj



+ iωβB, ð3Þ

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 ε μ −Je2 ϵ2jkJmk TM waves (2D) H2 ϵ2jkEk μ ε −Jm2 −ϵ2jkJek Quantum waves (3D) Ψ 2ℏ mi ∂jΨ 4 −4V_ℏω m 2ℏω

with wave numberk defined via

k2

= αβω2: ð4Þ

2.2. Unified Wave Equation in Matrix-Vector Form. We define a configuration with a preferred direction and reorga-nise equations (1) and (2) accordingly.

Consider a 3D spatial domainD, enclosed by surface ∂D. This surface consists of two planar surfaces∂D₀and∂D₁ per-pendicular to thex₃-axis and a cylindrical surface∂D_cylwith its axis parallel to thex₃-axis, see Figure 1. The surfaces∂D₀ and ∂D₁ are situated atx₃= x3,0 andx3= x3,1, respectively,

withx_3,1> x3,0. In general, these surfaces do not coincide with

physical boundaries. SurfaceS in Figure 1 is a cross section of D at arbitrary x3. The parametersαðxÞ and βðxÞ are piecewise

continuous smoothly varying functions inD, with discontin-uous jumps only at interfaces that are perpendicular to thex₃ -axis (hence,P and Q₃are continuous over the interfaces). In the lateral direction, the domain D can be bounded or unbounded. WhenD 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 = Q₃= 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 x₂= 0 and lowercase Greek subscripts take on the value 1 only.

We reorganise equations (1) and (2) into a matrix-vector wave equation which acknowledges thex₃-direction as the direction of preference. By eliminating the lateral compo-nentsQ₁andQ₂(or, for 2D wave problems, the lateral com-ponentQ₁), we obtain [8, 15, 19–21]

∂3q = Aq + d, ð5Þ

where wave vector q and source vector d are defined as q = P Q₃ ! , d = C3 B₀ ! , ð6Þ with B₀= B + 1 iω∂ν 1 βCν ð7Þ

and operator matrixA defined as A = 0 A12 A21 0 ! , ð8Þ with A12= iωβ, ð9Þ A21= iωα − 1 iω∂ν_β1∂ν: ð10Þ

The notation in the right-hand side of equations (7) and (10) should be understood in the sense that di fferen-tial 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 vector q are continuous over interfaces perpendicular to thex₃-axis. Moreover, these quantities constitute the power-flux density (or, for quantum-mechanical waves, the probability-flux density) in thex₃-direction via

j =1

4 P

∗_Q 3+ Q∗3P

f g, ð11Þ

where the asterisk denotes complex conjugation.

2.3. Symmetry Properties of the Operator Matrix. We discuss the symmetry properties of the operator matrix A. First, consider a general operator U (which can be a scalar or a matrix), containing space-dependent parameters and differ-ential operators ∂_ν. We introduce the transpose operator Ut_{via the following integral property:}

ð SðUfÞ t_{g dx} L= ð Sf t_Ut_g_dx L: ð12Þ

Here, xLis the lateral coordinate vector, with xL= ðx1, x2Þ

for 3D and xL= x1 for 2D wave problems. S denotes an

integration surface perpendicular to thex₃-axis at arbitrary

x₃, with edge ∂S, see Figure 1. The quantities f ðx_LÞ and

gðxLÞ 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 tof . When S is bounded,

x1 x2 _x 3 𝜕 1 𝜕 0 𝜕 cyl x3,0 x3,1 n = (0, 0, −1) n = (0, 0, 1) n = (n1, n2, 0) 𝜕

Figure 1: Configuration with the x3-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.

homogeneous Dirichlet or Neumann conditions are imposed at ∂S. When S is unbounded, ∂S has an infinite radius and f ðx_LÞ and gðx_LÞ are assumed to have sufficient decay alongS towards infinity. Operator U is said to be symmetric when Ut= U and skew-symmetric when Ut= −U. For the special case thatU = ∂_ν, 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 ofU) via the integral property

ð SðUfÞ †_{g dx} L= ð Sf †_U†_g_dx L: ð13Þ

When the test functions are vectors, f† is the complex conjugate transpose off ; when they are scalars, f†is the com-plex conjugate off . Operator U is said to be Hermitian (or self-adjoint) whenU†= U and skew-Hermitian when U†= −U. For the operators A12 and A21, defined in equations

(9) and (10), we find At₁₂= A12, A21t = A21, A†12= −A12,

andA†₂₁= −A21. Hence, operatorsA12andA21are

symmet-ric and skew-Hermitian. With these relations, wefind for the operator matrixA At N = −NA, ð14Þ A†_{K = −KA, ð15Þ} with N = 0 1 −1 0 ! , K = 0 1 1 0 ! : ð16Þ

Note that, using the expressions for q 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

j = 1

4q

†_{Kq: ð17Þ}

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 andB, characterised by wave vectors qAðx, ωÞ and qBðx, ωÞ,

obeying wave equation (5), with source vectors dAðx, ωÞ

and dBðx, ωÞ. 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 ∂₃ðqt_ANqBÞ in domain D.

Applying the product rule for differentiation, using equation (5) for both states, integrating the result overD and applying the theorem of Gauss yields

ð D  AqA ð Þt_{+ d}t A   NqB+ qtAN Aqð B+ dBÞ  dx = ð ∂Dq t ANqBn3dx: ð18Þ

Here,n₃is the component parallel to thex₃-axis of the outward pointing normal vector on ∂D, with n₃= −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 ∂D₀ and∂D₁, which together are denoted by∂D_0,1. The integral on the left-hand side can be written as Ð_Dð⋯Þdx =Ðx3,1

x3,0dx3

Ð

Sð⋯ÞdxL. Using equation

(12) for the integral alongS and symmetry property (14), it follows that the two terms in equation (18) containing oper-atorA cancel each other. Hence, we are left with

ð D d t ANqB+ qtANdB   dx = ð ∂D0,1 qtANqBn3dxL: ð19Þ

This is a convolution-type reciprocity theorem [22–24], because products like qtAðx, ωÞNqBðx, ωÞ in the frequency

domain correspond to convolutions in the time domain. A more familiar form is obtained by substituting the expres-sions for q, d, and N (equations (6) and (16)), choosing Cj

= 0 and using equation (2) to eliminate Q3, which gives ð Dð−BAPB+ PABBÞdx = ð ∂D0,1 1 iωβðPA∂3PB− ∂ð 3PAÞPBÞn3dxL: ð20Þ Next, consider the quantity∂₃ðq†AKqBÞ in domain D.

Fol-lowing the same steps as above, using equations (13) and (15) instead of (12) and (14), we obtain

ð D d † AKqB+ q†AKdB   dx = ð ∂D0,1 q†AKqBn3dxL: ð21Þ

This is a correlation-type reciprocity theorem [25], because products like q†Aðx, ωÞKqBðx, ωÞ in the frequency

domain correspond to correlations in the time domain. Substituting the expressions for q, d, and K and choosing

C_j= 0 yield the more familiar form

ð D B ∗ APB+ P∗ABB ð Þdx = ð ∂D0,1 1 iωβðP∗A∂3PB− ∂ð 3PAÞ∗PBÞn3dxL: ð22Þ We obtain a special case by choosing statesA and B iden-tical. Dropping the subscriptsA and B in equations (21) and (22) and multiplying both sides of these equations by 1/4 give

1 4 ð D d †_{Kq + q}†_Kd   dx =1 4 ð ∂D0,1 q†Kqn3dxL, ð23Þ

1 4 ð D B ∗_{P + P}∗_B ð Þdx =1 4 ð ∂D_0,1 1 iωβðP∗∂3P − ∂ð 3PÞ∗PÞn3dxL, ð24Þ 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 matrixA into a somewhat different form. To this end, we introduce an operatorH₂, according to

H2= −iω ffiffiffi β p A21 ffiffiffi β p = k2+pffiffiffiβ∂ν_β1∂ν ffiffiffi β p , ð25Þ

with operatorA₂₁defined in equation (10) and wavenumber

k in equation (4). Operator H2can be rewritten as a

Helm-holtz operator [14, 21]

H2= k2s + ∂ν∂ν, ð26Þ

with the scaled wavenumberk_sdefined as [26]

k2 s= k2−3 ∂ð νβÞð∂νβÞ 4β2 + ∂ν∂νβ ð Þ 2β : ð27Þ

Note thatHt₂= H2andH†2= H2; hence, operatorH2is

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 A = 0 iωβ − 1 iωp Hffiffiffiβ 2 1ffiffiffi β p 0 0 B @ 1 C A: ð28Þ

Next, we decompose this operator matrix as follows

A = LHL−1, ð29Þ with H = iH1 0 0 −iH1 ! , ð30Þ L = L1 L1 L2 −L2 ! , ð31Þ L−1₌ 1 2 L−1 1 L−12 L−1 1 −L−12 ! : ð32Þ

OperatorsH₁,L₁, andL₂are pseudodifferential oper-ators [7, 8, 14, 16, 21, 27–30]. The decomposition expressed by equation (29) is not unique; hence, different choices for operatorsH₁,L₁, andL₂ 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

ωβ = L1H1L−12 , ð33Þ

1

ωp Hffiffiffiβ 2 1ffiffiffi β

p = L2H1L−11 : ð34Þ

We introduce a decomposedfield vector p and a decom-posed source vector s via

q = Lp,  p = L−1q, ð35Þ d = Ls,  s = L−1d, ð36Þ where p = P + P− ! , s = S + S− ! : ð37Þ

Substitution of equations (29), (35), and (36) into the matrix-vector wave equation (5) yields

∂3p = H − L−1∂3L

 

p + s: ð38Þ Substituting equations (30), (31), (32), and (37) into equation (38) gives ∂3 P+ P− ! = iH1 0 0 −iH1 ! _P+ P− ! −1 2 L−1 1 L−12 L−1 1 −L−12 !  ∂3L1 ∂3L1 ∂3L2 −∂3L2 ! _P+ P− ! + S + S− ! : ð39Þ 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 wavefields P+_andP−_{propagate in the positive}

and negative x₃-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} x₃-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

ð D s t ALtNLpB+ ptALtNLsB   dx = ð ∂D0,1 ptALtNLpBn3dxL, ð40Þ ð D s † AL†KLpB+ p†AL†KLsB   dx = ð ∂D0,1 p†AL†KLpBn3dxL, ð41Þ 1 4 ð D s †_L†_{KLp + p}†_L†_KLs   dx =1 4 ð ∂D0,1 p†L†KLpn3dxL: ð42Þ

These equations form the basis for reciprocity theorems for the decomposedfield and source vectors p and s in the next two sections.

3.2. Flux-Normalised Decomposition and Reciprocity Theorems. The first choice of operators H₁, L₁, and L₂ obeying equations (33) and (34) is [14–16]

H1= H1/22 , ð43Þ

L1= ω/2ð Þ1/2β1/2H−1/21 , ð44Þ

L2= 2ωð Þ−1/2β−1/2H1/21 : ð45Þ

OperatorH₁, which is the square root of the Helmholtz operatorH₂, is commonly known as the square root opera-tor [3, 4, 8]. Like the Helmholtz operaopera-torH₂, the square root operatorH₁is a symmetric operator [16], henceHt₁= H1.

For the adjoint square root operator, we have H†₁= ðHt₁Þ∗ = H∗1. The spectrum ofH1 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 H†₁≈ H₁. Similar relations hold for the square root of the square root operator and its inverse; hence, ðH±1/2

1 Þ

t

= H±1/2 1 , and

neglecting evanescent waves, ðH±1/2₁ Þ†≈ H±1/2₁ . From here onward, we replace≈ by = when the only approximation is the negligence of evanescent waves. Using these symmetry relations for H₁ and equations (16), (31), (44), and (45), we obtain

Lt

NL = −N, ð46Þ

and neglecting evanescent waves,

L†_{KL = J, ð47Þ} with J = 1 0 0 −1 ! : ð48Þ

Hence, equations (40), (41), and (42) simplify to −ð D s t ANpB+ ptANsB   dx = − ð ∂D0,1 ptANpBn3dxL, ð D s † AJpB+ p†AJsB   dx = ð ∂D0,1 p†AJpBn3dxL, 1 4 ð D s †_{Jp + p}†_Js   dx =1 4 ð ∂D0,1 p†Jpn3dxL: ð49Þ

By substituting the expressions for p, s, N, and J (equa-tions (37), (16), and (48)), we obtain

ð D −S + AP−B+ S−AP+B− P+AS−B+ P−AS+B ð Þdx = ð ∂D0,1 −P+ AP−B+ P−AP+B ð Þn₃dxL, ð50Þ ð D S +∗ A P+B−S−∗A P−B+P+∗A S+B−P−∗A S−B ð Þdx = ð ∂D0,1 P+∗ A P+B−P−∗A P−B ð Þn₃dxL, ð51Þ 1 4 ð D S +∗_P+_−S−∗_P−_+P+∗_S+_−P−∗_S− ð Þdx =1 4 ð ∂D0,1 P+ 2 _{− P} − 2   n₃dxL: ð52Þ 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∂D_0,1. Therefore, we callP+andP− flux-normalised one-way wave fields. Conse-quently, 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 wavefield imaging methods for active and pas-sive data [31–42].

3.3. Field-Normalised Decomposition and Reciprocity Theorems. The second choice of operatorsH₁,L1, andL2

obeying equations (33) and (34) is [21]

H1= β1/2H1/22 β−1/2, ð53Þ

L1= 1, ð54Þ

L2= ωβð Þ−1H1: ð55Þ

Only the Helmholtz operator H₂ is the same as in the previous section (it is defined in equation (26)). The

operatorsH₁, L₁, and L₂ are different from those in the previous section, but for convenience, we use the same sym-bols. Using q = Lp (equation (35)) and equations (6), (31), (37), and (54), wefind

P = P+_{+ P}−_{: ð56Þ}

Hence,P+_{and P}− _{have the same physical dimension as the}

fullfield 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 operator H1/2

2 is symmetric, but H1

defined in equation (53) is not. From this equation, it easily follows thatH₁premultiplied byβ−1is symmetric, hence

1

βH1

t

= 1

βH1, ð57Þ

and neglecting evanescent waves, 1

βH1

†

= 1

βH1: ð58Þ

Using these symmetry relations for ð1/βÞH1 and

equa-tions (16), (31), (54), and (55), we obtain Lt NL = 0 −2L2 2L2 0 ! = −N _ωβ2 H1
 = − _ωβ2 H1
t N, ð59Þ

and neglecting evanescent waves, L†KL = 2L2 0 0 −2L2 ! = J 2 ωβH1
 = _ωβ2 H1
† J: ð60Þ

Using this in equations (40) and (41) yields −ð D s t A _ωβ2 H1
t NpB+ ptAN _ωβ2 H1
 sB  dx = − ð ∂D0,1 ptA _ωβ2 H1
t NpBn3dxL, ð D s † A _ωβ2 H1
† JpB+ p†AJ _ωβ2 H1
 sB " # dx = ð ∂D0,1 p†A _ωβ2 H1
† JpBn3dxL: ð61Þ

By substituting the expressions for p, s, N, and J (equations (37), (16), and (48)), using equations (12) and (13), we obtain −ð D 2 ωβ  H1S+A ð ÞP− B− Hð 1S−AÞP+B + P+AðH1S−BÞ− P−AðH1S+BÞdx = − ð ∂D0:1 2 ωβ  H1P+A ð ÞP−B− Hð 1P−AÞP+Bn3dxL, ð62Þ ð D 2 ωβ  H1S+A ð Þ∗P+ B− Hð 1S−AÞ∗P−B + P+∗A ðH1S+BÞ−P−∗A ðH1S−BÞdx = ð ∂D0,1 2 ωβðH1P+AÞ∗P+B− Hð 1P−AÞ∗P−Bn3dxL: ð63Þ

We aim to remove the operator H₁ from these equa-tions. From equations (39) and (54), we obtain

∂3P+= +iH1P+− 1 2 L −1 2 ∂3L2   P+_{− P}− ð Þ + S+_{, ð64Þ} ∂3P−= −iH1P−+ 1 2 L −1 2 ∂3L2  _P+ − P− ð Þ + S−, ð65Þ with L₂ defined in equation (55). Assuming that in state

A the derivatives in the x₃-direction of the parameters α and β at ∂D_0,1 vanish and there are no sources at ∂D_0,1, we find from equations (64) and (65)

∂3P±A= ±iH1P±Aat ∂D0,1: ð66Þ

Below we use this to remove H₁ 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), wefind S±_{= ±}1 2 1 ωβH1
−1 B₀+1 2C3, ð67Þ or ±_ωβ2 H1S±= B0± 1 ωβH1C3: ð68Þ

We define new decomposed sources B+

0 andB−0, accord-ing to B± 0= B0± 1 ωβH1C3= ± 2 ωβH1S±: ð69Þ

Using equations (66) and (69) in the right- and left-hand sides of equations (62) and (63), we obtain

ð D −B + 0,AP−B− B−0,AP+B+ P+AB−0,B+ P−AB+0,B   dx = ð ∂D0,1 −2 iωβ  ∂3P+A ð ÞP− B+ ∂ð 3P−AÞP+Bn3dxL, ð70Þ ð D B +∗ 0,AP+B+ B−∗0,AP−B+P+∗A B+0,B+P−∗A B−0,B   dx = ð ∂D0,1 −2 iωβð∂3P+AÞ∗P+B+ ∂ð 3P−AÞ∗P−Bn3dxL: ð71Þ 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 stateA. In the next section, we use equations (70) and (71) to derive repre-sentation theorems forfield-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 wavefields (equations (70) and (71)). First, we introduce the full Green’s function Gðx, x_A, ωÞ as a solution of the uni-fied wave equation (3) for a unit monopole point source at xA,

withBðx, ωÞ = δðx − x_AÞ and C_jðx, ωÞ = 0. Hence,

β∂j 1_β∂jG



+ k2G = iωβδ x − xð AÞ: ð72Þ

As boundary condition, we impose the radiation condition (i.e., outward propagating waves at infinity). Next, we intro-duce 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 Bðx, ωÞ = δðx − xAÞ

and Cjðx, ωÞ = 0. Using equations (69) and (7), we define

decomposed sources as B±₀= B±= B = ±2L2S±, with L2

defined in equation (55), or S±_{ðx, ωÞ = ±}1 2L −1 2 B±ðx, ωÞ = ± 1 2L −1 2 B x, ωð Þ = ±1 2L −1 2 δ x − xð AÞ: ð73Þ

We consider two sets of one-way Green’s functions. For the first set, we choose a point source S+ðx, ωÞ = ð1/2ÞL−1₂

B+_{ðx, ωÞ, with B}+_{ðx, ωÞ = δðx − x}

AÞ, 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 ∂3G+,+= +iH1G+,+− 1 2 L −1 2 ∂3L2   _G+,+ −G−,+ ð Þ +1 2L −1 2 δ x − xð AÞ, ð74Þ ∂3G−,+= −iH1G−,++ 1 2 L −1 2 ∂3L2  _G+,+ −G−,+ ð Þ: ð75Þ

Here, G±,+ _{stands for G}±,+_{ðx, x}

A, ωÞ. The second

super-script (+) indicates that the source at xA emits waves in the

positive x₃-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/2Þ}

L−1

2 B−ð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 ∂3G+,−= +iH1G+,−− 1 2 L −1 2 ∂3L2  _G+,− − G−,− ð Þ, ð76Þ ∂3G−,−= −iH1G−,−+ 1 2 L −1 2 ∂3L2   _G+,− − G−,− ð Þ −1 2L −1 2 δ x − xð AÞ: ð77Þ Here, G±,− stands for G±,−ðx, xA, ωÞ, with the second

superscript (−) indicating that the source at xAemits waves

in the negativex₃-direction. Like for the full Green’s function

Gðx, xA, ωÞ, we impose radiation conditions for both sets of

one-way Green’s functions.

Tofind a relation between the full Green’s function and the one-way Green’s functions, we evaluate β∂3ð1/βÞ∂3

ðG+,+_{+ 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 thusfind

G = G+,+_{+ G}−,+_{+ G}+,−_{+ G}−,−_{: ð78Þ}

This very simple relation is a consequence of the field-normalised decomposition, introduced in Section 3.3. 4.2. Source-Receiver Reciprocity. We derive source-receiver reciprocity relations for thefield-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 domainD, the parameters α and β are the same in the two states (see Section 2.4). OutsideD, these parameters may be different in the two states. For the Green’s state, we choose the parameters forx₃≤ x_3,0and forx₃≥ x_3,1 independent of the x₃-coordinate, according to αðx_LÞ 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∂D_0,1 (because outside D no scattering occurs along the x₃-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, x}

B, ωÞ = 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) van-ishes. Hence, taking into account thatB±₀= B± (sinceCj= 0),

equation (70) simplifies to ð D −B + AP−B− B−AP+B+ P+AB−B+ P−AB+B ð Þdx = 0: ð79Þ

First, we consider sources emitting waves in the positive

x₃-direction in both Green’s states, hence B+

A= δðx − xAÞ,

B−

A= 0, P±A= G±,+ðx, xA, ωÞ, B+B= δðx − xBÞ, B−B= 0, and P±B

= G±,+ðx, xB, ωÞ. Substituting this into equation (79) yields

G−,+ _x

B, xA, ω

ð Þ = G−,+ _x

A, xB, ω

ð Þ, ð80Þ

see Figure 2(a). Next, we replace the source in stateB by one emitting waves in the negative x₃-direction, hence B+

B= 0,

B−

B= δðx − xBÞ, and P±B= G±,−ðx, xB, ωÞ. This gives

G+,+ _x

B, xA, ω

ð Þ = G−,− _x

A, xB, ω

ð Þ, ð81Þ

see Figure 2(b). By replacing also the source in stateA by one emitting waves in the negative x3-direction, according to B+ A= 0, B−A= δðx − xAÞ, and P±A= G±,−ðx, xA, ωÞ, we obtain G+,− xB, xA, ω ð Þ = G+,− xA, xB, ω ð Þ, ð82Þ

see Figure 2(c). Finally, changing the source in state B back to the one emitting waves in the positive x₃ -direc-tion yields G−,− _x B, xA, ω ð Þ = G+,+ _x A, xB, ω ð Þ, ð83Þ see Figure 2(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 interchan-ging the source and the receiver. The absence of sign changes in equations (80), (81), (82), and (83) is due to the definition ofB±

0in equation (69). Moreover, unlike theflux-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 func-tion. By separately summing the left- and right-hand sides of equations (80), (81), (82), and (83) and using equation (78), we simply obtain

G xð B, xA, ωÞ = G xð A, xB, ωÞ: ð84Þ

4.3. Kirchhoff-Helmholtz Integrals for Forward Propagation. We derive Kirchhoff-Helmholtz integrals of the convolution type forfield-normalised one-way wave fields. For state B, we

G−,+_(x B, xA, 𝜔) G−,+(xA, xB, 𝜔) xB xA xB xA = x1 x2 _x 3 (a) xB xA xA xB G+,+_(x B, xA, 𝜔) G−,−(xA, xB, 𝜔) = x1 x2 _x 3 (b) xB xA xB xA = G+,−_(x A, xB, 𝜔) G+,−_(x B, xA, 𝜔) x1 x2 _x 3 (c) xB xA xA xB = G−,−_(x B, xA, 𝜔) G+,+(xA, xB, 𝜔) x1 x2 _x 3 (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.

consider the decomposed actualfield, with sources only out-sideD; hence, B±

0,B= 0 in D and P±B= P±ðx, ωÞ. The

parame-ters α and β are the actual parameters inside as well as outsideD. For state A, we choose the Green’s state with a unit point source at xAinD. The parameters α and β in D are the

same as those in stateB, but for x₃≤ x_3,0 and for x₃≥ x_3,1, they are chosen independent of the x₃-coordinate. Hence, the condition for the validity of equation (66) is fulfilled. First, we consider a source in stateA which emits waves in the positivex₃-direction, henceB+A= δðx − xAÞ, B−A= 0, and

P±

A= G±,+ðx, xA, ωÞ. Substituting all this into equation (70)

(withB±_0,A= B±A) gives

P− _x A, ω ð Þ = ð ∂D0,1 2 iωβ xð Þ  ∂3G+,+ðx, xA, ωÞ ð ÞP−ðx, ωÞ + ∂ð 3G−,+ðx, xA, ωÞÞP+ðx, ωÞ _n 3dxL: ð85Þ Next, we replace the source in state A by one which emits waves in the negative x₃-direction, hence B+

A= 0,

B−

A= δðx − xAÞ, and P±A= G±,−ðx, xA, ωÞ. Equation (70) thus

gives P+ _x A, ω ð Þ = ð ∂D0,1 2 iωβ xð Þ  ∂3G+,−ðx, xA, ωÞ ð ÞP−ðx, ωÞ + ∂ð 3G−,−ðx, xA, ωÞÞP+ðx, ωÞ _n 3dxL: ð86Þ Recall that ∂D_0,1 consists of ∂D₀ (with n₃= −1) and

∂D1 (with n3= +1), see Figure 1. Since G+,±ðx, xA, ωÞ = 0

at ∂D₀ and G−,±ðx, xA, ωÞ = 0 at ∂D1 (because outside D

no scattering occurs along the x₃-coordinate in state A), the first term under the integral in equations (85) and (86) gives a contribution only at∂D₁ and the second term only at ∂D₀. Hence, P± _x A, ω ð Þ =ð ∂D0 −2 iωβ xð Þ ∂3G−,∓ðx, xA, ωÞ  _P+ x, ω ð ÞdxL + ð ∂D1 2 iωβ xð Þ∂3G+,∓ðx, xA, ωÞP−ðx, ωÞdxL: ð87Þ

Note that there is no contribution fromP−ðx, ωÞ at ∂D₀ nor fromP+_{ðx, ωÞ at ∂D}

1, see Figure 3.

We conclude this section by considering a special case. Suppose the source of the actual field (state B) is located at xB in the half-space x3< x3,0. Then, by taking x3,1→

∞, the field P− _{at ∂D}

1 vanishes. This leaves the

single-sided representation P± _x A, xB, ω ð Þ = ð ∂D0 −2 iωβ xð Þ ∂3G−,∓ðx, xA, ωÞ   P+ _{x, x} B, ω ð Þdx_L: ð88Þ Note that we included the source coordinate vector xB

in the argument list ofP±_ðx

A, xB, ωÞ. This representation is

an extension of a previously derived result [43], in which the fields were decomposed at ∂D₀ 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 ∂D₀). 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 state B, 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 inD. The parameters α and β in

D are the same as those in state B, but for x3≤ x3,0 and for x₃≥ x3,1, they are chosen independent of thex3-coordinate.

Hence, the condition for the validity of equation (66) is ful-filled. First, we consider sources emitting waves in the posi-tive x₃-direction in both states, hence B+

A= δðx − xAÞ,

B−

A= 0, P±A= G±,+ðx, xA, ωÞ, B+B= δðx − xBÞsðωÞ, B−B= 0, and

P±

B= P±,+ðx, xB, ωÞ. Substituting this into equation (71) (with

B±

0,A= B±AandB±0,B= B±B) gives

P+,+ _x A, xB, ω ð Þ+χ xð Þ GB f +,+ðxB, xA, ωÞg∗s ωð Þ = ð ∂D0,1 −2 iωβ xð Þ  ∂3G+,+ðx, xA,ωÞ f g∗P+,+ _{x, x} B, ω ð Þ + ∂f 3G−,+ðx, xA,ωÞg∗P−,+ðx, xB, ωÞn3dxL, ð89Þ

whereχ is the characteristic function of the domain D. It is defined as χ xð Þ =B 1, for xBin D, 1 2, for xBon ∂D0,1, 0, for xBoutside D: 8 > > > < > > > : ð90Þ Since G+,+_{ðx, x} A, ωÞ = 0 at ∂D0 and G−,+ðx, xA, ωÞ = 0

at ∂D₁ (because outside D no scattering occurs along the

x₃-coordinate in stateA), the first term under the integral

xA x x P+_{(x, 𝜔)} P−_{(x, 𝜔)} x1 x2 _x 3 x3,0 x3,1 G−,−_{(x, x} A, 𝜔) G−,+_{(x, x} A, 𝜔) G+,+_{(x, x} A, 𝜔) G+,−_{(x, x} A, 𝜔) 𝜕 1 𝜕 0

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 toP+_ðx

in equation (89) gives a contribution only at ∂D₁ and the second term only at ∂D₀. Hence,

P+,+ _x A, xB, ω ð Þ + χ xð Þ GB f +,+ðxB, xA, ωÞg∗s ωð Þ = ð ∂D0 2 iωβ xð Þf∂3G−,+ðx, xA,ωÞg∗P−,+ðx, xB, ωÞdxL −ð ∂D1 2 iωβ xð Þf∂3G+,+ðx, xA,ωÞg∗P+,+ðx, xB, ωÞdxL: ð91Þ Next, we replace the source in state B by one emit-ting waves in the negative x₃-direction, hence B+B= 0,

B−

B= δðx − xBÞsðωÞ and P±B= P±,−ðx, xB, ωÞ. This gives

P+,− _x A, xB, ω ð Þ + χ xð Þ GB f −,+ðxB, xA, ωÞg∗s ωð Þ = ð ∂D0 2 iωβ xð Þf∂3G−,+ðx, xA,ωÞg∗P−,−ðx, xB, ωÞdxL −ð ∂D1 2 iωβ xð Þf∂3G+,+ðx, xA,ωÞg∗P+,−ðx, xB, ωÞdxL: ð92Þ By replacing also the source in state A by one emitting waves in the negative x₃-direction, according toB+ A= 0, B−A= δðx − xAÞ, and P±A= G±,−ðx, xA, ωÞ, we obtain P−,− _x A, xB, ω ð Þ + χ xð Þ GB f −,−ðxB, xA, ωÞg∗s ωð Þ = ð ∂D0 2 iωβ xð Þf∂3G−,−ðx, xA,ωÞg∗P−,−ðx, xB, ωÞdxL −ð ∂D1 2 iωβ xð Þf∂3G+,−ðx, xA,ωÞg∗P+,−ðx, xB, ωÞdxL: ð93Þ Finally, changing the source in stateB back to the one emitting waves in the positivex₃-direction yields

P−,+ _x A, xB, ω ð Þ + χ xð Þ GB f +,−ðxB, xA, ωÞg∗s ωð Þ = ð ∂D0 2 iωβ xð Þf∂3G−,−ðx, xA,ωÞg∗P−,+ðx, xB, ωÞdxL −ð ∂D1 2 iωβ xð Þf∂3G+,−ðx, xA,ωÞg∗P+,+ðx, xB, ωÞdxL: ð94Þ Equation (93) is an extension of a previously derived result [44], in which the fields were decomposed at ∂D_0,1 but not at xAand xB. Equations (91), (92), and (94) are

fur-ther 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, x}

B, ωÞ from ∂D1to xA. Except for some special cases,

the integrals along ∂D₁ do not vanish by takingx_3,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.

We conclude this section by considering a special case. Suppose that in state B the parameters α and β are the same as in stateA not only in D but also outside D. Then,

P±,±_{ðx, x}

B, ωÞ = 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 xBis located inD, we obtain

G_hðxA, xB, ωÞ = ð ∂D0 2 iωβ xð Þf∂3G−ðx, xA, ωÞg∗G−ðx, xB, ωÞdxL − ð ∂D1 2 iωβ xð Þf∂3G+ðx, xA, ωÞg∗G+ðx, xB, ωÞdxL, ð95Þ where the so-called homogeneous Green’s function G_hðx_A, xB, ωÞ is defined as

G_hðxA, xB, ωÞ = G xð A, xB, ωÞ + G∗ðxA, xB, ωÞ

= 2R G xf ð A, xB, ωÞg, ð96Þ

(withR denoting the real part) and where G±_{ðx, x}

A, ωÞ =

G±,+_{ðx, x}

A, ωÞ+G±,−ðx, xA, ωÞ (and a similar expression for

G±_{ðx, x}

B, ωÞ). Equation (95) is akin to the well-known

repre-sentation 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 conse-quence 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 f∂3G±,±ðx, xA, ωÞg∗ 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 ∂D0and ∂D1,

enclos-ing the focal point xA. Here, we discuss single-sided

field-normalised focusing functionsf±₁ðx, xA, ωÞ and we use these

in modifications of reciprocity theorems (70) and (71) to

xA x x xB {G−,−_{(x, x} A, 𝜔)} ⁎ P−,+_{(x, x} B, 𝜔) P+,+_{(x, x} B, 𝜔) {G+,−_{(x, x} A, 𝜔)} ⁎ 𝜕 1 𝜕 0 _x 3,0 x3,1 x1 x2 _x 3

Figure 4: Visualisation of the different terms in the field-normalised one-way Kirchhoff-Helmholtz integral for backward propagation, formulated by equation (94).

derive single-sided Kirchhoff-Helmholtz integrals for back-ward propagation.

We start by defining a new domain D_A, enclosed by two surfaces∂D₀ and ∂DAperpendicular to the x3-axis atx3= x_3,0 andx₃= x3,A, respectively, withx3,A> x3,0, see Figure 5.

Hence,∂D_A is chosen such that it contains the focal point xA. The two surfaces∂D0and∂DAare together denoted by

∂D0,A. The focusing functions f±1ð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

x₃-coordinate for x₃≤ x_3,0 and for x₃≥ x_3,A. Hence, the condition for the validity of equation (66) is fulfilled. Analogous to equation (56), the field-normalised focusing functions f±₁ðx, xA, ωÞ are related to the full focusing

function f₁ðx, x_A, ωÞ, according to

f₁ðx, xA, ωÞ = f+1ðx, xA, ωÞ + f−1ðx, xA, ωÞ: ð97Þ

The focusing function f+₁ðx, xA, ωÞ is incident to the

domain D_A from the half-space x₃< x3,0 (see Figure 5).

It propagates and scatters in the inhomogeneous domain DA, focuses at xA on surface ∂DA, and continues as

f+

1ðx, xA, ωÞ in the half-space x3> x3,A. The back-scattered

field leaves DAvia surface∂D0and continues asf−1ðx, xA, ωÞ

in the half-space x₃< x3,0. The focusing conditions at the

focal plane ∂DA are [18]

∂3f+1ðx, xA, ωÞ  x3=x3,A= 1 2iωβ xð ÞA δ xð L− xL,AÞ, ð98Þ ∂3f−1ðx, xA, ωÞ ½ _x 3=x3,A= 0: ð99Þ

Here, xL,Adenotes the lateral coordinates of xA. The

oper-ators ∂₃ and the factor ð1/2ÞiωβðxAÞ are not necessary to

define the focusing conditions but are chosen for later conve-nience. To avoid instability, evanescent waves are excluded from the focusing functions. This implies that the delta func-tion in equafunc-tion (98) should be interpreted as a spatially band-limited delta function. Note that the sifting property of the delta function,hðx_L,AÞ =Ð_Sδðx_L− x_L,AÞhðx_LÞdx_L, remains

valid for a spatially band-limited delta function, assuming

hðxLÞ is also spatially band-limited.

We now derive single-sided Kirchhoff-Helmholtz inte-grals for backward propagation. We consider the reciprocity theorems for field-normalised one-way wave fields (equa-tions (70) and (71)), with D and ∂D_0,1 replaced byD_Aand

∂D0,A, respectively. For state A, we consider the focusing

functions discussed above; hence, B+_Aðx, ωÞ = B−_Aðx, ωÞ = 0 and P±

Aðx, ωÞ = f±1ðx, xA, ωÞ. For state B, we consider the

decomposed actual field, with a point source at xB in the

half-space x₃> x3,0and source spectrum sðωÞ. The

parame-tersα and β in state B are the actual parameters inside as well as outside∂D_0,A. First, we consider a source in stateB which emits waves in the positive x₃-direction, hence B+_Bðx, ωÞ =

δðx − xBÞsðωÞ, B−Bðx, ωÞ = 0, and P±Bðx, ωÞ = P±,+ðx, xB, ωÞ.

Substituting all this into equations (70) and (71) (with

B±

0= B±), using equations (98) and (99) in the integrals

along ∂DA, gives P−,+ _x A, xB, ω ð Þ + χAð ÞfxB −1ðxB, xA, ωÞsð Þω = ð ∂D0 2 iωβ xð Þ  ∂3f+1ðx, xA,ωÞ  _P−,+ x, xB, ω ð Þ + ∂ð 3f−1ðx, xA,ωÞÞP+,+ðx, xB, ωÞ  dxL, ð100Þ P+,+ _x A, xB, ω ð Þ− χAð Þ fxB +1ðxB, xA, ωÞ  ∗_{s ωð Þ} = ð ∂D0 −2 iωβ xð Þ  ∂3f+1ðx, xA,ωÞ  ∗_P+,+ _{x, x} B, ω ð Þ + ∂f 3f−1ðx, xA,ωÞg∗P−,+ðx, xB, ωÞ  dxL, ð101Þ

where χ_A is the characteristic function of the domain DA.

It is defined by equation (90), with D and ∂D_0,1 replaced by DA and ∂D0,A, respectively. Next, we replace the source

in state B by one which emits waves in the negative

x₃-direction, hence B+

Bðx, ωÞ = 0, B−Bðx, ωÞ = δðx − xBÞsðωÞ,

andP±_Bðx, ωÞ = P±,−ðx, xB, ωÞ. This gives

P−,− _x A, xB, ω ð Þ + χAð ÞfxB +1ðxB, xA, ωÞsð Þω = ð ∂D0 2 iωβ xð Þ  ∂3f+1ðx, xA,ωÞ  _P−,− x, xB, ω ð Þ + ∂ð 3f−1ðx, xA,ωÞÞP+,−ðx, xB, ωÞ  dxL, ð102Þ P+,− _x A, xB, ω ð Þ− χAð Þ fxB f −1ðxB, xA, ωÞg∗s ωð Þ = ð ∂D0 −2 iωβ xð Þ  ∂3f+1ðx, xA,ωÞ  ∗_P+,− x, xB, ω ð Þ + ∂f 3f−1ðx, xA,ωÞg∗P−,−ðx, xB, ωÞ  dxL: ð103Þ

Equations (100), (101), (102), and (103) are single-sided representations for backward propagation of the one-way wave fields P±,±_{ðx, x}

B, ωÞ from ∂D0 to xA. Similar results

have been previously obtained [47, 48], but without decom-position at xB. An advantage of these equations over

equa-tions (91), (92), (93), and (94) is that the backward propagated fields P±,±ðxA, xB, ωÞ are expressed entirely in

terms of integrals along the surface∂D₀.

f+_{(x, x} A, 𝜔) 1 x1 x2 x3 𝜕 A 𝜕 0 A f+_{(x, x} A, 𝜔) 1 f−1(x, xA, 𝜔) x3,0 x3,A xA

Figure 5: Configuration for the derivation of the single-sided Kirchhoff-Helmholtz integrals for backward propagation.