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

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 ﬂ ux-normalised and ﬁ eld-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 ﬁ elds, for both types of normalisation. Based on the ﬁ eld-normalised reciprocity theorems, we derive representation theorems for decomposed wave ﬁ elds. In particular, we derive double- and single-sided Kirchho ﬀ -Helmholtz integrals for forward and backward propagation of decomposed wave ﬁ elds. The single-sided Kirchho ﬀ -Helmholtz integrals for backward propagation of ﬁ eld-normalised decomposed wave ﬁ elds ﬁ nd applications in re ﬂ ection imaging, accounting for multiple scattering.


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][4][5][6][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][9][10][11][12][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 Section 2, 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.

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 quantummechanical waves, in a lossless potential) is governed by the following two equations in the space-frequency domain: 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 ∂/∂x j , and Einstein's summation convention applies to repeated subscripts. Pðx, ωÞ and Q j ðx, ωÞ are space-and frequency-dependent wave field quantities, αðxÞ and βðxÞ are real-valued space-dependent parameters, and Bðx, ωÞ and C j ðx, ωÞ 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 = ðx 1 , x 2 , x 3 Þ is the 3D Cartesian coordinate vector and lowercase Latin subscripts take on the values 1, 2, and 3. For the 2D situation, x = ðx 1 , x 3 Þ 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 j Q j are continuous over the interface. Here, n j represents the components of the normal vector n = ðn 1 , n 2 , n 3 Þ at the interface for the 3D situation or n = ðn 1 , n 3 Þ for the 2D situation.
Advances in Mathematical Physics with wave number k defined via

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 domain D, enclosed by surface ∂D. This surface consists of two planar surfaces ∂D 0 and ∂D 1 perpendicular to the x 3 -axis and a cylindrical surface ∂D cyl with its axis parallel to the x 3 -axis, see Figure 1. The surfaces ∂D 0 and ∂D 1 are situated at x 3 = x 3,0 and x 3 = x 3,1 , respectively, with x 3,1 > x 3,0 . In general, these surfaces do not coincide with physical boundaries. Surface S in Figure 1 is a cross section of D at arbitrary x 3 . The parameters αðxÞ and βðxÞ are piecewise continuous smoothly varying functions in D, with discontinuous jumps only at interfaces that are perpendicular to the x 3 -axis (hence, P and Q 3 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 = Q 3 = 0 or n ν ∂ ν P = n ν ∂ ν Q 3 = 0 at ∂D cyl , 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 ∂D cyl has an infinite radius and we assume that P and Q 3 have "sufficient decay" at infinity. For the 2D situation, the configuration is a cross section of the 3D situation for x 2 = 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 the x 3 -direction as the direction of preference. By eliminating the lateral components Q 1 and Q 2 (or, for 2D wave problems, the lateral component Q 1 ), we obtain [8,15,[19][20][21] where wave vector q and source vector d are defined as with and operator matrix A defined as with 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 vector q are continuous over interfaces perpendicular to the x 3 -axis. Moreover, these quantities constitute the power-flux density (or, for quantum-mechanical waves, the probability-flux density) in the x 3 -direction via where the asterisk denotes complex conjugation.

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 differential operators ∂ ν . We introduce the transpose operator U t via the following integral property: Here, x L is the lateral coordinate vector, with x L = ðx 1 , x 2 Þ for 3D and x L = x 1 for 2D wave problems. S denotes an integration surface perpendicular to the x 3 -axis at arbitrary x 3 , with edge ∂S, see Figure 1. The quantities f ðx L Þ and gðx L Þ are space-dependent test functions (scalars or vectors). When these functions are vectors, f t is the transpose of f ; when they are scalars, f t is equal to f . When S is bounded,  Figure 1: Configuration with the x 3 -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 x 2 = 0.

Advances in Mathematical Physics
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 along S towards infinity. Operator U is said to be symmetric when U t = U and skew-symmetric when U t = −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 operator U † (i.e., the complex conjugate transpose of U) via the integral property 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 A 12 and A 21 , defined in equations (9) and (10) Hence, operators A 12 and A 21 are symmetric and skew-Hermitian. With these relations, we find for the operator matrix A with 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 probabilityflux density) as 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 q A ðx, ωÞ and q B ðx, ωÞ, obeying wave equation (5), with source vectors d A ðx, ωÞ and d B ð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 ∂ 3 ðq t A Nq B Þ 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 Here, n 3 is the component parallel to the x 3 -axis of the outward pointing normal vector on ∂D, with n 3 = −1 at ∂D 0 , n 3 = +1 at ∂D 1 , and n 3 = 0 at ∂D cyl , see Figure 1. In the following, the integral on the right-hand side is restricted to the horizontal surfaces ∂D 0 and ∂D 1 , which together are denoted by ∂D 0,1 . The integral on the left-hand side can be written as 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 This is a convolution-type reciprocity theorem [22][23][24], because products like q t A ðx, ωÞNq B ðx, ωÞ 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 C j = 0 and using equation (2) to eliminate Q 3 , which gives Next, consider the quantity ∂ 3 ðq † A Kq B Þ in domain D. Following the same steps as above, using equations (13) and (15) instead of (12) and (14), we obtain This is a correlation-type reciprocity theorem [25], because products like q † A ðx, ωÞKq B ð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 We obtain a special case by choosing states A 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 Advances in Mathematical Physics respectively. These equations quantify conservation of power (or, for quantum-mechanical waves, probability).

Decomposed Wave Equation and Its Symmetry Properties
To facilitate the decomposition of the matrixvector wave equation (equation (5)), we recast the operator matrix A into a somewhat different form. To this end, we introduce an operator H 2 , according to with operator A 21 defined in equation (10) and wavenumber k in equation (4). Operator H 2 can be rewritten as a Helmholtz operator [14,21] with the scaled wavenumber k s defined as [26] Note that H t 2 = H 2 and H † 2 = H 2 ; hence, operator H 2 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 Next, we decompose this operator matrix as follows with Operators H 1 , L 1 , and L 2 are pseudodifferential operators [7,8,14,16,21,[27][28][29][30]. The decomposition expressed by equation (29) is not unique; hence, different choices for operators H 1 , L 1 , and L 2 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 We introduce a decomposed field vector p and a decomposed source vector s via where Substitution of equations (29), (35), and (36) into the matrix-vector wave equation (5) yields Substituting equations (30), (31), (32), and (37) into equation (38) gives 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 fields P + and P − propagate in the positive and negative x 3 -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 3 -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 5 Advances in Mathematical Physics These equations form the basis for reciprocity theorems for the decomposed field and source vectors p and s in the next two sections.

Field-Normalised Decomposition and Reciprocity
Theorems. The second choice of operators H 1 , L 1 , and L 2 obeying equations (33) and (34) is [21] Only the Helmholtz operator H 2 is the same as in the previous section (it is defined in equation (26)). The 6 Advances in Mathematical Physics operators H 1 , L 1 , and L 2 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 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 − fieldnormalised one-way wave fields (for convenience, we use the same symbols as in the previous section). The square root operator H 1/2 2 is symmetric, but H 1 defined in equation (53) is not. From this equation, it easily follows that H 1 premultiplied by β −1 is symmetric, hence and neglecting evanescent waves, Using these symmetry relations for ð1/βÞH 1 and equations (16), (31), (54), and (55), we obtain and neglecting evanescent waves, Using this in equations (40) and (41) yields By substituting the expressions for p, s, N, and J (equations (37), (16), and (48)), using equations (12) and (13), we obtain We aim to remove the operator H 1 from these equations. From equations (39) and (54), we obtain with L 2 defined in equation (55). Assuming that in state A the derivatives in the x 3 -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) Below we use this to remove H 1 from the right-hand sides of equations (62) and (63). Next, we aim to remove H 1 from the left-hand sides of these equations. From s = L −1 d (equation (36)) and equations (6), (32), (37), (54), and (55), we find or We define new decomposed sources B + 0 and B − 0 , according to

Advances in Mathematical Physics
Using equations (66) and (69) in the right-and lefthand sides of equations (62) and (63), we obtain ð 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.

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 oneway wave fields (equations (70) and (71)). First, we introduce the full Green's function Gðx, x A , ωÞ as a solution of the unified wave equation (3) for a unit monopole point source at x A , with Bðx, ωÞ = δðx − x A Þ and C j ðx, ωÞ = 0. Hence, 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 x A . Hence, we choose again Bðx, ωÞ = δðx − x A Þ and C j ðx, ωÞ = 0. Using equations (69) and (7), we define decomposed sources as B ± 0 = B ± = B = ±2L 2 S ± , with L 2 defined in equation (55), or 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 2 B + ðx, ωÞ, with B + ðx, ωÞ = δðx − x A Þ, which emits waves from x A in the positive x 3 -direction, and we set S − ðx, ωÞ equal to zero. Hence, for this first set, one-way equations (64) and (65) become Here, G ±,+ stands for G ±,+ ðx, x A , ωÞ. The second superscript (+) indicates that the source at x A emits waves in the positive x 3 -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 − x A Þ, which emits waves from x A in the negative x 3 -direction, and we set S + ðx, ωÞ equal to zero. Hence, for this second set, one-way equations (64) and (65) become Here, G ±,− stands for G ±,− ðx, x A , ωÞ, with the second superscript (−) indicating that the source at x A emits waves in the negative x 3 -direction. Like for the full Green's function Gðx, x A , ωÞ, we impose radiation conditions for both sets of one-way Green's functions.

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 x 3 ≤ x 3,0 and for x 3 ≥ x 3,1 independent of the x 3 -coordinate, according to αðx L Þ and First, we consider sources emitting waves in the positive see Figure 2(a). Next, we replace the source in state B by one emitting waves in the negative x 3 -direction, hence B + B = 0, B − B = δðx − x B Þ, and P ± B = G ±,− ðx, x B , ωÞ. This gives see Figure 2(b). By replacing also the source in state A by one emitting waves in the negative x 3 -direction, according to see Figure 2(c). Finally, changing the source in state B back to the one emitting waves in the positive x 3 -direction yields see Figure 2(d). Source-receiver reciprocity relations similar to equations (80), (81), (82), and (83) were previously derived for fluxnormalised 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)

Kirchhoff-Helmholtz Integrals for Forward Propagation.
We derive Kirchhoff-Helmholtz integrals of the convolution type for field-normalised one-way wave fields. For state B, we x 1 x 2 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. 9 Advances in Mathematical Physics consider the decomposed actual field, with sources only outside D; hence, B ± 0,B = 0 in D and P ± B = 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 x A in D. The parameters α and β in D are the same as those in state B, but for x 3 ≤ x 3,0 and for x 3 ≥ x 3,1 , they are chosen independent of the x 3 -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 Next, we replace the source in state A by one which emits waves in the negative x 3 -direction, hence B + A = 0, B − A = δðx − x A Þ, and P ± A = G ±,− ðx, x A , ωÞ. Equation (70) thus gives Recall that ∂D 0,1 consists of ∂D 0 (with n 3 = −1) and ∂D 1 (with n 3 = +1), see Figure 1. Since G +,± ðx, x A , ωÞ = 0 at ∂D 0 and G −,± ðx, x A , ωÞ = 0 at ∂D 1 (because outside D no scattering occurs along the x 3 -coordinate in state A), the first term under the integral in equations (85) and (86) gives a contribution only at ∂D 1 and the second term only at ∂D 0 . Hence, Note that there is no contribution from P − ðx, ωÞ at ∂D 0 nor from P + ð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 x B in the half-space x 3 < x 3,0 . Then, by taking x 3,1 → ∞, the field P − at ∂D 1 vanishes. This leaves the singlesided representation Note that we included the source coordinate vector x B in the argument list of P ± ðx A , x B , ωÞ. This representation is an extension of a previously derived result [43], in which the fields were decomposed at ∂D 0 but not at x A . It describes forward propagation of the one-way field P + ðx, x B , ωÞ from the surface ∂D 0 to x A (with x A and x B defined at opposite sides of ∂D 0 ). In the following two sections, we discuss representations for backward propagation of oneway wave fields.

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 x B 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 x A in D. The parameters α and β in D are the same as those in state B, but for x 3 ≤ x 3,0 and for x 3 ≥ x 3,1 , they are chosen independent of the x 3 -coordinate. Hence, the condition for the validity of equation (66) is fulfilled. First, we consider sources emitting waves in the positive x 3 -direction in both states, hence where χ is the characteristic function of the domain D. It is defined as , for x B on ∂D 0,1 , 0, for x B outside D: Since G +,+ ðx, x A , ωÞ = 0 at ∂D 0 and G −,+ ðx, x A , ωÞ = 0 at ∂D 1 (because outside D no scattering occurs along the x 3 -coordinate in state A), the first term under the integral 10 Advances in Mathematical Physics in equation (89) gives a contribution only at ∂D 1 and the second term only at ∂D 0 . Hence, Next, we replace the source in state B by one emitting waves in the negative ÞsðωÞ and P ± B = P ±,− ðx, x B , ωÞ. This gives By replacing also the source in state A by one emitting waves in the negative x 3 -direction, according Finally, changing the source in state B back to the one emitting waves in the positive x 3 -direction yields Equation (93) is an extension of a previously derived result [44], in which the fields were decomposed at ∂D 0,1 but not at x A and x B . 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, x B , ωÞ from ∂D 0 and P +,± ðx, x B , ωÞ from ∂D 1 to x A . Except for some special cases, the integrals along ∂D 1 do not vanish by taking x 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 state A not only in D but also outside D. Then, P ±,± ðx, x B , ωÞ = G ±,± ðx, x B , ωÞ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 x B is located in D, we obtain where the so-called homogeneous Green's function G h ðx A , x B , ωÞ is defined as (with R denoting the real part) and where G ± ðx, x A , ωÞ = G ±,+ ðx, x A , ωÞ+G ±,− ðx, x A , ωÞ (and a similar expression for G ± ðx, x B , ωÞ). 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. (Single-Sided). The complex-conjugated Green's functions f∂ 3 G ±,± ðx, x A , ωÞg * under the integrals in equations (91), (92), (93), and (94) can be seen as focusing functions, which focus the wave fields P ±,± ðx, x B , ωÞ onto a focal point x A . However, this focusing process requires that these wave fields are available at two boundaries ∂D 0 and ∂D 1 , enclosing the focal point x A . Here, we discuss single-sided fieldnormalised focusing functions f ± 1 ðx, x A , ωÞ and we use these in modifications of reciprocity theorems (70) and (71) to We start by defining a new domain D A , enclosed by two surfaces ∂D 0 and ∂D A perpendicular to the x 3 -axis at x 3 = x 3,0 and x 3 = x 3,A , respectively, with x 3,A > x 3,0 , see Figure 5. Hence, ∂D A is chosen such that it contains the focal point x A . The two surfaces ∂D 0 and ∂D A are together denoted by ∂D 0,A . The focusing functions f ± 1 ðx, x A , ωÞ, 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 D A equal to those in the actual state B, and independent of the x 3 -coordinate for x 3 ≤ x 3,0 and for x 3 ≥ x 3,A . Hence, the condition for the validity of equation (66) is fulfilled. Analogous to equation (56), the field-normalised focusing functions f ± 1 ðx, x A , ωÞ are related to the full focusing function f 1 ðx, x A , ωÞ, according to

Kirchhoff-Helmholtz Integrals for Backward Propagation
The focusing function f + 1 ðx, x A , ωÞ is incident to the domain D A from the half-space x 3 < x 3,0 (see Figure 5). It propagates and scatters in the inhomogeneous domain D A , focuses at x A on surface ∂D A , and continues as f + 1 ðx, x A , ωÞ in the half-space x 3 > x 3,A . The back-scattered field leaves D A via surface ∂D 0 and continues as f − 1 ðx, x A , ωÞ in the half-space x 3 < x 3,0 . The focusing conditions at the focal plane ∂D A are [18] Here, x L,A denotes the lateral coordinates of x A . The operators ∂ 3 and the factor ð1/2Þiωβðx A Þ 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, 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ðx L Þ 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 ∂D 0,1 replaced by D A and ∂D 0,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, x A , ωÞ. For state B, we consider the decomposed actual field, with a point source at x B in the half-space x 3 > x 3,0 and source spectrum sðωÞ. The parameters α and β in state B are the actual parameters inside as well as outside ∂D 0,A . First, we consider a source in state B which emits waves in the positive x 3 -direction, hence B + B ðx, ωÞ = δðx − x B ÞsðωÞ, B − B ðx, ωÞ = 0, and P ± B ðx, ωÞ = P ±,+ ðx, x B , ωÞ. Substituting all this into equations (70) and (71) (with B ± 0 = B ± ), using equations (98) and (99) in the integrals along ∂D A , gives where χ A is the characteristic function of the domain D A . It is defined by equation (90), with D and ∂D 0,1 replaced by D A and ∂D 0,A , respectively. Next, we replace the source in state B by one which emits waves in the negative x 3 -direction, hence B + B ðx, ωÞ = 0, B − B ðx, ωÞ = δðx − x B ÞsðωÞ, and P ± B ðx, ωÞ = P ±,− ðx, x B , ωÞ. This gives Equations (100), (101), (102), and (103) are single-sided representations for backward propagation of the one-way wave fields P ±,± ðx, x B , ωÞ from ∂D 0 to x A . Similar results have been previously obtained [47,48], but without decomposition at x B . An advantage of these equations over equations (91), (92), (93), and (94) is that the backward propagated fields P ±,± ðx A , x B , ωÞ are expressed entirely in terms of integrals along the surface ∂D 0 .
x 3,0 x 3,A x A Figure 5: Configuration for the derivation of the single-sided Kirchhoff-Helmholtz integrals for backward propagation.