# [

Imaging non-reciprocal media]
Unified wave field imaging of inhomogeneous non-reciprocal media
Kees Wapenaar and Christian Reinicke]
Kees Wapenaar and Christian Reinicke

Delft University of Technology, P.O. Box 5048, 2600 GA Delft, The Netherlands

Acoustic imaging methods often ignore multiple scattering. This may lead to false images in cases where multiple scattering is strong. Marchenko imaging has recently been introduced as a data-driven way to deal with internal multiple scattering. Promising results have been obtained with geophysical and ultrasonic data. Given the increasing interest in non-reciprocal materials, both for acoustic and electromagnetic applications, we propose to modify the Marchenko method for imaging of such materials.

We start by formulating a unified wave equation for non-reciprocal materials, exploiting the similarity between acoustic and electromagnetic wave phenomena. This unified wave equation forms the basis for deriving reciprocity theorems that interrelate wave fields in a non-reciprocal medium and its adjoint. Next, we reformulate these theorems for downgoing and upgoing wave fields. From these decomposed reciprocity theorems we derive representations of the Green’s function inside the non-reciprocal medium, in terms of the reflection response at the surface and focusing functions inside the medium and its adjoint. These representations form the basis for deriving a modified version of the Marchenko method for imaging of non-reciprocal media. We illustrate the proposed method at the hand of the numerically modeled reflection response of a horizontally layered medium.

The appendices contain a detailed derivation of the unified wave equation for non-reciprocal media and the decomposition of the reciprocity theorems.

## 1 Introduction

Acoustic imaging methods are traditionally based on the single-scattering assumption (Claerbout, 1971; Stolt, 1978; Berkhout & van Wulfften Palthe, 1979; Williams & Maynard, 1980; Devaney, 1982; Bleistein & Cohen, 1982; Maynard et al., 1985; Langenberg et al., 1986; McMechan, 1983; Esmersoy & Oristaglio, 1988; Oristaglio, 1989; Norton, 1992; Wu, 2004; Lindsey & Braun, 2004; Etgen et al., 2009). Multiply scattered waves are not properly handled by these methods and may lead to false images overlaying the desired primary image. Several approaches have been developed that account for multiple scattering. For the sake of the discussion it is important to distinguish between different classes of multiply scattered waves. Waves that have scattered at least once at the surface of the medium are called surface-related multiples. This type of multiple scattering is particularly severe in exploration geophysics. However, because the scattering boundary is known, this class of multiples is relatively easily dealt with. Successful methods have been developed to suppress surface-related multiples prior to imaging (Verschuur et al., 1992; Carvalho et al., 1992; van Borselen et al., 1996; Biersteker, 2001; Pica et al., 2005; Dragoset et al., 2010). Waves that scatter several times inside the medium before being recorded at the surface are called internal multiples. Internal multiple scattering may occur at heterogeneities at many scales. We may distinguish between deterministic scattering at well-separated scatterers, giving rise to long period multiples, and diffuse scattering in stochastic media. Of course this distinction is not always sharp. In this paper we only consider the first type of internal multiple scattering, which typically occurs in layered media (which, in general, may have curved interfaces and varying parameters in the layers). Several imaging approaches that account for deterministic internal multiples are currently under development, such as the inverse scattering series approach (Weglein et al., 1997; Ten Kroode, 2002; Weglein et al., 2003), full wavefield migration (Berkhout, 2014; Davydenko & Verschuur, 2017), and Marchenko imaging. The latter approach builds on a 1D autofocusing procedure (Rose, 2001, 2002; Broggini & Snieder, 2012), which has been generalised for 2D and 3D inhomogeneous media (Wapenaar et al., 2012, 2014; Broggini et al., 2014; Behura et al., 2014; Meles et al., 2015; van der Neut et al., 2015; van der Neut & Wapenaar, 2016; Thorbecke et al., 2017; Van der Neut et al., 2017; Singh et al., 2017; Mildner et al., 2017; Elison et al., 2018). This methodology predicts the internal multiples in a data-driven way and suppresses their contribution to the final image. Promising results have been obtained with geophysical (Ravasi et al., 2016; Ravasi, 2017; Staring et al., 2018) and ultrasonic data (Wapenaar et al., 2018).

To date, the application of Marchenko imaging has been restricted to reciprocal media. With the increasing interest in non-reciprocal materials, both in electromagnetics (Willis, 2011; He et al., 2011; Ardakani, 2014) and in acoustics (Willis, 2012; Norris et al., 2012; Nassar et al., 2017; Attarzadeh & Nouh, 2018), it is opportune to modify the Marchenko method for imaging of non-reciprocal media. We start with a brief review of the wave equation for non-reciprocal media. By restricting this to scalar waves in a 2D plane, it is possible to capture different wave phenomena by a unified wave equation. Next, we formulate reciprocity theorems for waves in a non-reciprocal medium and its adjoint. From these reciprocity theorems we derive Green’s function representations, which form the basis for the Marchenko method in non-reciprocal media. We illustrate the new method with a numerical example, showing that it has the potential to accurately image a non-reciprocal medium, without false images related to multiple scattered waves.

## 2 Unified wave equation for non-reciprocal media

Consider the following unified equations for 2D wave propagation in the -plane in inhomogeneous, lossless, anisotropic, non-reciprocal media

(1) | |||

(2) |

These equations hold for transverse-electric (TE), transverse-magnetic (TM), horizontally polarised shear (SH) and acoustic (AC) waves. They are formulated in the space-time domain, with . Operator stands for differentiation in the direction. Lower-case subscripts and take the values 1 and 3 only; Einstein’s summation convention applies for repeated subscripts. Operator stands for temporal differentiation. The macroscopic wave field quantities ( and ), the effective medium parameters (, and ) and the macroscopic source functions ( and ) are specified for the different wave phenomena in Table 1 (note that ). For details we refer to Appendix A.

By eliminating from equations (1) and (2) we obtain a scalar wave equation for field quantity , according to

(3) |

see Appendix A for the derivation. Here is the inverse of . Compare equation (3) with the common wave equation for waves in isotropic reciprocal media

(4) |

In equation (3), replaces , with being responsible for the non-reciprocal behaviour. Moreover, replaces , thus accounting for anisotropy of the effective non-reciprocal medium.

## 3 Reciprocity theorems for a non-reciprocal medium and its adjoint

We derive reciprocity theorems in the space-frequency -domain for wave fields in a non-reciprocal medium and its adjoint. To this end, we define the temporal Fourier transform of a space- and time-dependent function as

(5) |

where is the angular frequency and the imaginary unit. For notational convenience we use the same symbol for quantities in the time domain and in the frequency domain. We use equation (5) to transform equations (1) and (2) to the space-frequency domain. The temporal differential operators are thus replaced by , hence

(6) | |||

(7) |

with , , and . A reciprocity theorem formulates a mathematical relation between two independent states (Fokkema & van den Berg, 1993; de Hoop, 1995; Achenbach, 2003). We indicate the sources, medium parameters and wave fields in the two states by subscripts and . Consider the quantity

(8) |

Applying the product rule for differentiation, using equations (6) and (7) for states and , using , integrating the result over domain enclosed by boundary with outward pointing normal vector and applying the theorem of Gauss, we obtain

(9) | |||

This is the general reciprocity theorem of the convolution type. When the medium parameters , and are identical in both states, then the first integral on the right-hand side vanishes, but the second integral, containing , does not vanish. This confirms that is the parameter responsible for the non-reciprocal behaviour. When we choose , then the second integral also vanishes. For this situation we call state , with parameters , and , the actual state, and state , with parameters , and , the adjoint state. We indicate the adjoint state by a superscript . Hence

(10) | |||

This reciprocity theorem will play a role in the derivation of Green’s function representations for the Marchenko method for non-reciprocal media (section 4). Here we use it to derive a relation between Green’s functions in states and . For the adjoint state we choose a unit monopole point source at in , hence , where is the Dirac delta function. The response to this point source is the Green’s function in state , hence . Similarly, for state we choose a unit monopole point source at in , hence and . We substitute these expressions into equation (10) and set the other source quantities, and , to zero. Further, we assume that Neumann or Dirichlet boundary conditions apply at , or that the medium at and outside is homogeneous and reciprocal. In each of these cases the boundary integral vanishes. We thus obtain (Slob & Wapenaar, 2009; Willis, 2012)

(11) |

The left-hand side is the response to a source at in the adjoint medium (with parameter ), observed by a receiver at . The right-hand side is the response to a source at in the actual medium (with parameter ), observed by a receiver at . Note the analogy with the flow-reversal theorem for waves in flowing media (Lyamshev, 1961; Godin, 1997; Wapenaar & Fokkema, 2004).

Next, we consider the quantity

(12) |

Superscript denotes complex conjugation. Following the same steps as before, we obtain

(13) | |||

This is the general reciprocity theorem of the correlation type. When the medium parameters , and are identical in both states, then the first and second integral on the right-hand side vanish. Hence

(14) | |||

Also this reciprocity theorem will play a role in the derivation of Green’s function representations for the Marchenko method for non-reciprocal media.

## 4 Green’s function representations for the Marchenko method

We use the reciprocity theorems of the convolution and correlation type (equations (10) and (14)) to derive Green’s function representations for the Marchenko method for non-reciprocal media. The derivation is similar to that for reciprocal media (Wapenaar et al., 2014); here we emphasise the differences. We consider a spatial domain , enclosed by two infinite horizontal boundaries and (with below ), and two finite vertical side boundaries (at ), see Figure 1. The positive -axis points downward. The normal vectors at and are and , respectively. The boundary integrals in equations (10) and (14) along the vertical side boundaries vanish (Wapenaar & Berkhout, 1989). Assuming there are no sources in in both states, the reciprocity theorems thus simplify to

(15) |

and

(16) |

For the derivation of the representations for the Marchenko method it is convenient to decompose the wave field quantities in these theorems into downgoing and upgoing fields in both states. Consider the following relations

(17) |

with wave vectors and defined as

(18) |

Here and are downgoing and upgoing wave fields, respectively. Operator in equation (17) is a pseudo-differential operator that composes the total wave field from its downgoing and upgoing constituents (Corones et al., 1983; Fishman et al., 1987; Wapenaar & Berkhout, 1989; Fishman, 1993; de Hoop, 1992, 1996; Wapenaar, 1996; Haines & de Hoop, 1996; Fishman et al., 2000). Its inverse decomposes the total wave field into downgoing and upgoing fields. For inhomogeneous isotropic reciprocal media, the theory for this operator is well developed. For anisotropic non-reciprocal media we restrict the application of this operator to the laterally invariant situation. In Appendix B we use equations (17) and (18) at boundaries and to recast reciprocity theorems (15) and (16) as follows

(19) |

and

(20) |

Note that the assumption of lateral invariance only applies to boundaries and ; the remainder of the medium (in- and outside ) may be arbitrary inhomogeneous. Equation (19) is exact, whereas in equation (20) evanescent waves are neglected at boundaries and .

In the following we define (at ) as the upper boundary of an inhomogeneous, anisotropic, non-reciprocal, lossless medium. Furthermore, we define (at , with ) as an arbitrary boundary inside the medium. We assume that the medium above is homogeneous. For state we consider a unit source for downgoing waves at , just above (hence, , with ). The response to this unit source at any observation point is given by , where and denote the downgoing and upgoing components of the Green’s function. For at , i.e., just below the source, we have and , with denoting the reflection response at of the medium below . At , we have . For state we consider a focal point at at . The medium in state is a truncated medium, which is identical to the actual medium between and , and homogeneous below . At a downgoing focusing function , with , is incident to the truncated medium. This function focuses at , hence, at we have . The response to this focusing function at is . Because the truncated medium is homogeneous below , we have at . The quantities for both states are summarised in Table 2.

Note that the downgoing focusing function , for at , is the inverse of the transmission response of the truncated medium (Wapenaar et al., 2014), hence

(21) |

for at . To avoid instabilities in the evanescent field, the focusing function is in practice spatially band-limited.

Substituting the quantities of Table 2 into equations (19) and (20) gives

(22) |

and

(23) |

respectively. These are two representations for the upgoing and downgoing parts of the Green’s function between at the acquisition surface and inside the non-reciprocal medium. They are expressed in terms of the reflection response and a number of focusing functions. Unlike similar representations for reciprocal media (Slob et al., 2014; Wapenaar et al., 2014), the focusing functions in equation (22) are defined in the adjoint of the truncated medium. Therefore we cannot use the standard approach to retrieve the focusing functions and Green’s functions from the reflection response . We obtain a second set of representations by replacing all quantities in equations (22) and (23) by the corresponding quantities in the adjoint medium. For the focusing functions in equation (22) this implies they are replaced by their counterparts in the truncated actual medium. We thus obtain

(24) |

and

(25) |

respectively. Because in practical situations we do not have access to the reflection response in the adjoint medium, we derive a relation analogous to equation (11) for this reflection response. To this end, consider the quantities in Table 3, with and just above , and with denoting a boundary below all inhomogeneities, so that there are no upgoing waves at . Substituting the quantities of Table 3 into equation (19) (with replaced by ) gives

(26) |

Equations (22) (25), with replaced by , form the basis for the Marchenko method, discussed in the next section.

Table 3: Quantities to derive equation (26). at at 0

## 5 The Marchenko method for non-reciprocal media

The standard multidimensional Marchenko method for reciprocal media (Slob et al., 2014; Wapenaar et al., 2014) uses the representations of equations (22) and (23), but without the superscript , to retrieve the focusing functions from the reflection response. Here we discuss how to modify this method for non-reciprocal media. We separate the representations of equations (22) (25) into two sets, each set containing focusing functions in one and the same truncated medium. These sets are equations (23) and (24), with the focusing functions in the truncated actual medium, and equations (22) and (25), with the focusing functions in the truncated adjoint medium. We start with the set of equations (23) and (24), which read in the time domain (using equation (26))

(27) |

and

(28) |

respectively. We introduce time windows to remove the Green’s functions from these representations. Similar as in the reciprocal situation, we assume that the Green’s function and the time-reversed focusing function on the left-hand side of equation (27) are separated in time, except for the direct arrivals (Wapenaar et al., 2014). This is a reasonable assumption for media with smooth lateral variations, and for limited horizontal source-receiver distances. Let denote the travel time of the direct arrival of . We define a time window , where is the Heaviside function and a small positive time constant. Under the above-mentioned assumption, we have . For the focusing function on the left-hand side of equation (27) we write (Wapenaar et al., 2014)

(29) | |||||

where is the inverse of the direct arrival of the transmission response of the truncated medium and the scattering coda. The travel time of is and the scattering coda obeys for . Hence, . Applying the time window to both sides of equation (27) thus yields

(30) |

We assume that the Green’s function and the focusing function in the left-hand side of equation (28) are separated in time (without overlap). Unlike for reciprocal media, we need a different time window to suppress the Green’s function, because the latter is defined in the adjoint medium. To this end we define a time window , where denotes the travel time of the direct arrival in the adjoint medium. Applying this window to both sides of equation (28) yields

(31) |

Equations (30) and (31), with given by equation (29), form a set of two equations for the two unknown functions and (with at ). These functions can be resolved from equations (30) and (31), assuming , , , and are known for all and at . The reflection responses and are obtained from measurements at the upper boundary of the medium. This involves deconvolution for the source function, decomposition and, when the upper boundary is a reflecting boundary, elimination of the surface-related multiple reflections (Verschuur et al., 1992). The travel times and , and the inverse of the direct arrival of the transmission response, , can be derived from a background model of the medium and its adjoint. A smooth background model is sufficient to derive these quantities, hence, no information about the scattering interfaces inside the medium is required. The iterative Marchenko scheme to solve for and reads

(32) | |||||

(33) |

with

(34) |

starting with . Once and are found, is obtained from equation (29) and, subsequently, the Green’s functions and are obtained from equations (27) and (28). Note that only is defined in the actual medium. To obtain in the actual medium we consider the set of equations (22) and (25), which read in the time domain (using equation (26))

(35) |

and

(36) |

respectively. The same reasoning as above leads to the following iterative Marchenko scheme for the focusing functions in the truncated adjoint medium

(37) | |||||

(38) |

with

(39) |

starting with . Here can be derived from the adjoint background model. Once the focusing functions and are found, the Green’s functions and are obtained from equations (35) and (36).

We conclude this section by showing how and can be used to image the interior of the non-reciprocal medium. First we derive a mutual relation between these Green’s functions. To this end, consider the quantities in Table 4. Here in state is the reflection response at of the adjoint medium below , with defined just above and the medium in state being homogeneous above . Substituting the quantities of Table 4 into equation (19) (with and replaced by and , respectively) gives