TRAC method in dissipative media—a first analysis in frequency domain and homogeneous media

We propose to explore the time-reversed absorbing condition (TRAC) method in the case of dissipative homogeneous media. In previous work, the TRAC method was derived from the time-reversibility of the (undamped) wave equation and proved to be efficient in both the time-domain and the frequency-domain. Namely, two main utilisations of the TRAC method have been probed: (a) redatuming, i.e. moving virtually the measurements by reconstructing the wavefield and (b) tracking down the location of a possible inclusion inside the domain. In this paper, we focus on the redatuming application and investigate the feasibility of the TRAC method in the case of dissipation. In particular, we will see that performing the classical TRAC method, i.e. ignoring the dissipation, may give satisfactory results, even for larger values of dissipation. An analysis is provided in the frequency-domain and one-space dimension and shows satisfactory updated versions of the TRAC method. Moreover, a systematic error study in two-space dimension is illustrated via numerical examples.


Introduction
Time reversal (TR) was introduced in the 90s by Matthias Fink as a response to imaging problems [1,2]. The principle of TR is to take advantage of the time reversibility of the undamped wave equation to reconstruct the wavefield from measurements recorded on a receivers' array, also called time reversal mirror (TRM) [3,4]. These data are time-reversed and sent back into the domain, creating a wavefield propagating back, and through, any heterogeneous components included in the domain. This has been widely advertised as 'recreating the past' in the TR community.
Many work on TR has been done since the seminal papers. In particular, even though TR is a powerful tool and is robust even with respect to complex media, the generated time-reversed wavefield back-propagating in the ambient medium to its original source has a non-decreasing energy due to the diffraction limit [1,5]. Therefore, a focal spot is observed at the location of the source, which prevents a good resolution of the image. In [6,7], the authors proposed a way to limit the effects of the diffraction limit while working with TRMs. To tackle this issue, a study about the TR operator decomposition was achieved to gain a better understanding of the technique [8,9]. Other techniques using filtering also examined ways to reduce the diffraction limit and improve the resolution of the image [10][11][12] in the context of heterogeneous media, summarised in 2013 in [13]. Some theoretical results were proposed for dissipative media in [14,15] as well.
About ten years ago, an alternative was proposed in order to remove the issue of the diffraction limit, via a numerical method called TRAC for time reversed absorbing condition [16][17][18]. The main idea of this method is to introduce an artificial sink, also called trial domain, in the computational domain that encloses the source or the scatterer acting as a passive source. The choice of the boundary condition to impose on the new artificial boundary is crucial and it has been proven that a time-reversed absorbing boundary condition (ABC) was the right condition [16], hence the name of the method. Then, the type and order of the chosen ABC totally depends on the shape of the artificial boundary and the level of accuracy one wants to reach.
There has been many work on absorbing or non-reflecting boundary conditions, since the first papers [19][20][21], namely for different shapes like in [22][23][24] for an ellipse and [25,26] for arbitrary shapes. A wide literature also consider the question of the accuracy by developing always higher-order of approximation of the ABC, see for instance [27][28][29][30][31][32]. Moreover, the ABC actually depends on the model as well. Therefore, many types of ABCs have been designed for the acoustic waves in the case of heterogeneous media [29,33] and multiple scattering [34][35][36]. Then, for the Maxwell equations, non-reflecting boundary conditions are proposed in [37][38][39], and for the elastic waves see [40][41][42][43].
The TRAC method has two main applications. The first application is the ability to reduce the size of the computational domain by moving virtually the measured data from the receivers' array to the artificial trial domain. This principle is known as redatuming and was first developed for seismic waves [44] and later generalised in [45]. This first application has been explored for the TRAC method in [46,47] and coupled with inverse problem techniques to corroborate the fact that the new gathered virtual data provide a good reconstruction of the expected image in a considerably reduced computational domain. This naturally improved the cost in terms of computational time and storage for the subsequent inverse problem. More recently, full waveform redatuming has been explored via a least-squares version of the TRAC method in [48].
The second application of the TRAC method consists of tracking inclusions or scatterers. The method is using the information provided via the reconstruction of the wavefield to determine the location of the scatterer. It has been indeed shown that a bad position of the trial domain, i.e. not enclosing entirely the scatterer, would lead to an erroneous reconstruction of the wavefield. By considering the associated error, the position of the trial domain can be corrected to locate the inclusion better. This has been studied in [17,18,49], but also for crack and source identification in [50][51][52]. Recently, an algorithm for an automatic detection of the inclusion has been proposed in [53]. This application is out of the scope of this paper, and we refer to [17,18,53] for the feasibility of detecting a scatterer with the TRAC method.
The novelty of this paper is to consider a dissipative version of the TRAC method while performing redatuming, that is, for the reconstruction of the wavefield from measurements inside the computational domain. Until now, the computational domain was always assumed homogeneous, constant and non-dissipative, which limited its utilisation in real applications containing dissipation. Therefore, we propose a study of the method when an anti-dissipative term is added by TR. We will use simple ABCs including the dissipative component, derived in appendix, since it seems that very little work has been proposed on that topic. Only recently, Barucq et al [54] designed ABCs in the context of 3D atmospheric waves, whose model contains a damping term. Our 1D and 2D results are consistent with the ABCs derived in [54].
In the rest of the paper, we consider the wave equation in the frequency domain 1 , which will be explicitly expressed using the Helmholtz equation. By doing so, we are able to compute by hand (or using symbolic calculators such as Maple) solutions to the forward problem and diverse versions of the time reversed counterparts of the forward problem.
This paper is structured as follows: a reminder of the principle of the TRAC method in the frequency domain, followed by research questions studied through this paper, are presented in section 2. In section 3 we propose an analysis in one-space dimension and frequency domain. Then, the two-space dimension case will be illustrated in section 4 by numerical results. Finally, after the conclusion, a derivation of the ABCs in the presence of dissipation is given in appendix.

Principle of the TRAC method
Before describing the TRAC method, we will recall the associated forward problem and state the studied inverse problem.

Forward problem
We consider an inclusion D living in some known medium, which can be assumed homogeneous or not, dissipative or not. A source f generates an incident wavefield u I that impinges on the inclusion, whose scattering response u S is recorded on a receivers' array. The resulting forward total wavefield u T := u I + u S satisfies: with Sommerfeld radiation conditions. In (1), ω denotes the frequency and d is the dimension, likely to be 1, 2 or 3. We will only consider the case of a single frequency in this paper. Parameter κ represents the physical properties of the medium and the inclusion, and can be defined as: Figure 1. Geometry of the domain: (a) the total field u T , as sum of the incident wave u I and the scattering response u S , is recorded on the receivers' array Γ; (b) once recorded and digitally time-reversed, the time-reversed total wavefield u T R is sent back from Γ into the computational domain. In the inverse problem, the inclusion D may be known or (more likely) not.
where c is the wave propagation speed, σ the dissipation and i = √ −1. Both c and σ are assumed to be functions in space with possible discontinuity at the interface between the inclusion and the surrounding medium. Therefore, κ can denote any type of inclusions: sound-soft, sound-hard or penetrable. Moreover, we use the following notation: the subscript 0 is used for the surrounding medium while the subscript D is for the inclusion. The characteristic function χ D defines the area delimited by the boundary of the inclusion D. Thus, Remark 2.1. For some applications using electromagnetic waves, the Helmholtz equation is used as a good approximation [56][57][58]. In this case, the parameters involved are the electric permittivity ε, the magnetic permeability µ and the electric conductivity σ, the latter inducing damping. With these notations, κ becomes:

Inverse problem
In this paper, our goal is to perform redatuming from acquired measurements, that is, we wish to reconstruct the scattered (or total) wavefield from the only knowledge of this wavefield at the recording location. To do so, we introduce a bounded computational domain Ω. We also denote by Γ the receivers' array, which may be part of ∂Ω or not, and which may be continuous or discrete, see figure 1. The studied inverse problem can be stated as follows: Let ω be a given fixed frequency. Given measurements of the total wavefield on Γ , reconstruct the total wavefield everywhere inside the computational domain Ω. TR offers a solution to the above inverse problem: the measurements are digitally timereversed and sent back into the computational domain from Γ. Obviously, in the time-harmonic case there is no time to be reversed, but this corresponds to the complex conjugate of the total wavefield via Fourier Transform. Abusively, we will keep talking about time-reversed field, when the data is processed digitally either in time or in frequency. As a result, the time-reversed total wavefield u T R , where the index R stands for reverse, back-propagates to its original source, in theory. To get a perfect time-reversed total wavefield, we need to assume that the inclusion D is known, which is usually not the case, even though we may assume some approximate location of it. The time-reversed problem reads: with an artificial boundary condition on the remaining ∂Ω\Γ and an appropriate boundary condition on ∂D corresponding to the type of inclusion, if known.
To overcome the fact that D is unknown, the classical approach is to back-propagate the wavefield inside the computational domain assuming there is no inclusion, i.e. Ω\D is taken as the entire Ω in (3), see [5] and references therein. This will obviously lead to an incorrect reconstruction.

Time reversed absorbing conditions
Our aim is to reconstruct the total wavefield from measurements using TR techniques [1,2]. To perform the TRAC method, we introduce an artificial convex domain B, for instance a ball, which is assumed to enclose inclusion D, see figure 2. As already developed in previous papers [17,18], the time-reversed total wavefield v T R can be reconstructed when solving: with appropriate radiation conditions on the boundary ∂Ω, if relevant. Here, the operator TRAC denotes an appropriate time-reversed ABC imposed on ∂B for the purpose of the method. The choice of the ABC is independent of the type and shape of inclusion D that generated the scattered wavefield. Note that TR in the frequency domain corresponds to complex conjugation [17], which is denoted by an overline here. The underlying inverse problem can be updated into: Let ω be a given fixed frequency. Given measurements of the total wavefield on Γ , reconstruct the total wavefield everywhere inside the computational domain Ω \ B. If by any chance, the location of B with respect to D is not satisfied, then the reconstruction will not hold independently of the chosen TRAC operator [17,18,53]. ■

Research questions
The TRAC method has already been extensively studied in the time domain [16-18, 48, 51-53], but only mentioned quickly in the frequency domain [17,59]. TR relies on the reversibility in time of the undamped wave equation. Therefore, the first studies never mentioned dissipation or damping, which deprives the equation of the time reversibility. However, in many applications, dissipation is present to some extent and the surrounding medium may not be homogeneous. That is why, we propose to work on the TRAC method for more general media. It is well-known that a damping term in the wave equation gives an attenuated solution over time, modelled by a decreasing exponential. As a result, TR of such a term will induce an increasing exponential in the solution. Hence, it is legitimate to wonder if a dissipation can be considered for the TRAC method. Since we are working on a finite time frame (in the case of the time-domain wave equation), it is likely that the solution does not have time to blow up and the increasing exponential will not prevent a good reconstruction [60]. In the frequency domain, we do not have an increasing exponential in time anymore, but the related increasing amplitude is still concerning. Therefore, the following research questions will be investigated. Question 2.1. To which extent does the anti-dissipative term in the TRAC method affect the reconstruction of the total (or scattered) wavefield?
A numerical analysis in one-space and two-space dimensions will allow us have a better understanding of the error made using the dissipative version of the TRAC method. While deriving the dissipative TRAC method, multiple choices of time-reversed absorbing condition will be available, raising the next question: Question 2.2. What can an ABC look like in case of dissipation?
In this paper, we will derive ABCs for dissipative media and observe that the resulting differential operator is hardly translatable from the frequency domain to the time domain via an inverse Fourier transform. As a result, a last 2 question comes to mind: Question 2.3. If a non-dissipative TRAC problem is considered while there is clear dissipation in the experiment (physical or numeric), how wrong qualitatively and quantitatively is this approximation for the reconstruction of the total (or scattered) wavefield? That is, which order of approximation should we consider before the reconstruction is not acceptable?
To answer partly to the questions above, we will focus on time-harmonic TRAC for the study.

Analysis in one-space dimension (and frequency domain)
For this section, we assume that the physical characteristic functions, c (or ε, µ) and σ, are piecewise constant, i.e. constant in the surrounding medium and constant in the inclusion D. Working in one-space dimension reduces the Helmholtz equation to a linear constant coefficient second order ordinary differential equation that can be easily solved by hand to access analytical solutions. Therefore, we show that the TRAC method provides exact reconstructions of the total wavefield for the correct location of B with respect to D.
In a second part of this section, we also consider some approximations of the time reversed absorbing condition on the boundary of B. This consideration is meant to explore how wrong the TRAC method can be in presence of dissipation if the time reversed ABC is inexact. In particular, the exact dissipative absorbing condition is hardly translatable in the time domain. Therefore, we are tempted to use the non-dissipative condition, even though we know that the surrounding medium is dissipative, or a simpler expression of the time reversed absorbing condition that will be translatable in time domain. We aim to quantify the error made by these approximations.

Exact TRAC reconstruction in 1D
In one-space dimension, see figure 3, the acoustic total wavefield satisfies: where The derivation of the boundary conditions is given in appendix section 'In the frequency domain and one-space dimension (1D)'. Problem (5) is written for penetrable inclusions, i.e. inclusions that will also let part of the wavefield propagate through itself. For sake of completeness, we also consider fully reflective inclusions known as sound-hard (as if wave speed c D inside the inclusion was very large) and sound-soft (as if c D = 0). In practice, we work with the forward problem with the corresponding boundary conditions at the boundary of the scatterer (Neumann for sound-hard and Dirichlet for sound-soft): Now assume that we record data at is introduced and assumed to enclose inclusion D, see figure 4. Now, the computational domain Ω\B corresponds to the two disconnected segments (x L , x B − ) and (x B + , x R ). The TRAC problem to reconstruct the total wavefield from measurements satisfies: where u I is the incident wavefield satisfying (5) Parameter ω σ0 in the Helmholtz equation is now complex conjugated and only takes the properties of the surrounding medium. The last line of (8) corresponds to the time reversed measurements. Note that in the case of sound-soft or sound-hard inclusions, the time reversed solution (as well as the forward one) is constant equal to zero in too). Therefore, for sound-soft or sound-hard inclusions, we solve (8), and compare with the forward reference, only in (x B + , x R ).

Proposition 3.1. The general solution to (8), independently of the form of the incident and the total wavefield, is
where operators TRAC ± ex stand for Proof. All calculations can be made by hand or using a symbolic calculator such as Maple.
For the sake of space, we do not give the details here.  (8) gives the exact reconstruction of the total field, solution to (5) where the wavefield propagates inside the inclusion with the physical properties of the surrounding medium. ■

Approximate TRAC reconstruction in 1D
Working in one-space dimension and frequency domain allows for analytical solutions and exact calculations. In this section, the underlying idea is to study the effect of an approximate non-dissipative time reversed ABC when dissipation is involved in the forward problem. We can already see that the exact ABC in 1D and frequency domain is hardly translatable to the time domain, hence the need for an approximation of the ABC and a quantification of the inherent error. We will compare three levels of approximation to the exact TRAC problem (8): • The 1st-order dissipative TRAC method, using a first-order approximation of the timereversed ABC on ∂B, for very small values of σ 0 .
• The zeroth-order dissipative TRAC method, using a zeroth-order approximation of the timereversed ABC on ∂B, (namely the dissipation is omitted in the ABC) • The classical TRAC method, where the dissipation is omitted in the reconstruction, Let's take the following notation for the ABCs: and Using the proposed approximations above, we obtained the following three results:

. The first-order dissipative TRAC method:
A first-order approximation to the TRAC problem can be written as: The general solution is: with Similarly, the zeroth-order dissipative TRAC method can be expressed using the zerothorder ABC and we have the following result:

Proposition 3.3. The zeroth-order dissipative TRAC method:
A zeroth-order approximation (omitting dissipation on the ABCs) to the TRAC problem can be written as: The general solution is: with And finally, we propose to compute the reconstruction using the classical TRAC method, which does not contain any dissipation term.

Proposition 3.4 (The classical TRAC method).
In the eventuality that the dissipation is unknown altogether, the erroneous TRAC problem becomes: The general solution is: For the sake of illustration, we use the following parameters: 5, 6), and c 0 = 1, σ 0 = 0.7, c D = 1.5, σ D = 2. In figure 5, we plot the complex conjugated total wavefield for reference, together with the exact TRAC reconstruction (9), the first-order dissipative TRAC reconstruction (13), the zeroth-order dissipative TRAC reconstruction (15) and the classical TRAC reconstruction (17). We observe that all three approximations adding some kind of dissipation are superimposed in (x L , x B − ) ∪ (x B + , x R ) for this choice of parameters, which also reveals that (13) and (15) may give a good approximation, while the classical TRAC solution is totally wrong.
The TRAC method in the presence of dissipation requires an anti-dissipative term in the Helmholtz equation to be valid. We observe that, despite the anti-dissipative component, the solutions to the different variations of the TRAC method do not blow up or increase exponentially in a way preventing the reconstruction. We will now quantify the error made using the proposed approximations.
Let's consider the behaviour of the error as a function of the ambient dissipation σ 0 , so that we can determine the values of σ 0 for which the error is negligible or, on the contrary, that will give a very poor reconstruction. We consider the relative L 2 -error for this study where variant refers to the different possible approximations of the TRAC method listed above and we recall that u T is the complex conjugated of the forward total field for reference. Note that the L 2 -norm is taken either in , since the behaviour differs on each interval. In particular, for sound-soft or sound-hard inclusions (in the current settings), the reconstruction only occurs in For the numerical test, we consider the same parameters as for figure 5, except that σ 0 varies in [0, 5]. In figure 6, we display the relative L 2 -error calculated on the interval [x B + , x R ] We compare the error for the different types of inclusions. In the top row, the relative L 2 -error shows a maximum and then seems to converge towards zero at the infinity for the first and zeroth order dissipate TRAC methods. However, for the classical TRAC method, we see a strictly increasing error, which confirms that the classical TRAC method (omitting any dissipation) is not reliable. For the other two methods, we seek for the value of σ 0 giving the maximal error and the range of values for σ 0 so that the relative L 2 -error is below 1%. Note that the error is having its maximum for the same value of σ 0 independently of the type of inclusion. We summarise the results in table 1.
The observations on these results are the following: • For the first and zeroth order dissipative TRAC methods, the relative L 2 -error only present one maximum; For larger values of σ 0 , the error tends to zero, as the solution in the vicinity of x B + has an decaying amplitude due to the value of the dissipation. The reconstruction is indeed of good quality around the receiver x R , i.e. the wavefield has not propagated far enough yet to show major discrepancies. On the other hand, in the vicinity of the trial domain x B + , the reconstruction depends more on the choice of the time-reversed absorbing condition, leading to our second point. • In the vicinity of the trial domain x B + , the reconstructed wavefield may display a change of phase compared with the reference wavefield. This is clearly observed for the values of σ 0 where we have the peaks in the error curves. The phase may be off for larger values of σ 0 as well, but is not visible because of the decaying amplitude. • The error in the case of the classical TRAC method does not show the same pattern and is only strictly increasing: the main error comes from the amplitude that remains constant instead of decaying in [x B + , x R ] rapidly to zero as the dissipation is more important in the medium.
The bottom graph displayed in figure 6 is showing the relative L 2 -error in the vincinity of the other side of the trial domain, i.e. in a neighbourhood of x B − . We also focussed on computing the error in a small interval around the boundary, with ϵ = 0.5, to ensure that we would capture the discrepancies due to the time-reversed absorbing condition. This error is only computed for the penetrable case and all three approximations show a peak for very different values of σ 0 , before the amplitude becomes to small to account of any error. Unsurprisingly, the peak of the error is attained for smaller values of σ 0 as the accuracy of the method gets lower. By zooming in the solution, we indeed observe the same phenomena as described before: phase shift inducing the peak, and constant amplitude for the classical TRAC method. Remark 3.2. The observations made above concerned this specific test with the chosen parameters. However, it seems that the global behaviour remains even with other sets of parameters. Unfortunately, this has not been proven yet, as the formulae for the relative L 2 -error with respect to σ 0 could not be explicit and analysed.
Finally, we used the one-space dimension example to illustrate the effect of approximate TRAC methods on the reconstruction, while the exact dissipative TRAC method and all analytical solutions were actually available. This was essentially to give an account of the error in known and explicit grounds. In the next section, we investigate the two-space dimension case, still in the frequency domain, where explicit solutions are not available anymore for a clear error analysis.

Numerical experiments in two-space dimensions (and frequency domain)
As it is well known, 2D is not the best candidate for exact formulations. Therefore, we do not present any proof for Claim 2.1. However, a good approximation for the ABC may induce a good enough reconstruction of the total (or scattered) wavefield. For the first studies on the TRAC method, the first order Bayliss-Turkel boundary condition was used successfully for elliptic trial domain [16][17][18]. More recently, the second order Enquist-Madja boundary condition was also proposed in [53] with excellent results, as well.
In this section, we follow [17,53], and consider both a first-order Bayliss-Turkel-like ABC and a second-order Engquist-Madja-like ABC to perform the TRAC method in 2D in dissipative (and homogeneous) media.

A numerical experiment in 2D-forward problem
To illustrate the feasibility of the TRAC method with dissipation in two-space dimension and frequency domain, we propose a benchmark example. For this, we assume that the surrounding medium is constant and homogeneous, but contains a dissipative component. Also, we use the Helmholtz equation with parameters from electromagnetism: electric permittivity ε, magnetic permeability µ and electric conductivity σ, see remark 2.1. Furthermore, we denote by ε 0 and µ 0 the electric permittivity and the magnetic permeability, respectively, in the vacuum.
Consider a square computational domainΩ, which contains one or several inclusions. We assume that the surrounding medium is homogeneous and has a relative permittivity ε r = 5, a relative permeability µ r = 1 and an electric conductivity σ 0 = exp(−3) ∼ 0.05. These values correspond to the permittivity, permeability and conductivity of the soil [61]. The size ofΩ is taken as 2 m in width and 2 m in depth. The inclusions are located somewhere around the centre ofΩ.
To solve the forward problem, we use the finite element method through the software FreeFem++ [62], namely P 2 elements, and the total-field/scattered-field (TF/SF) formulation [56,57] with second-order Engquist-Madja boundary conditions, see derivation in appendix section 'In the frequency domain and two-space dimension (2D)' of the appendix. A point source is modelled by where H (2) 0 is the zeroth-order Hankel function of Type 2, (x s , y s ) are the coordinates of the point source and k σ := ωµ 0 µ r (ωε 0 ε r − iσ 0 ).
We work with different types of inclusions, still following the parameters from [61]: • sound-hard, with ε r = 1, µ r = 10 4 and σ = 10 7 . These would for instance correspond to iron mines buried in the land; • penetrable, with ε r = 3, µ r = 1 and σ = 0, to model mines with plastic wrapping.
We also consider two sizes for the inclusions: a small one with 10 × 4 cm 2 and a large one with 60 × 10 cm 2 .
In figure 7, we show the numerical results for the forward problem in the case of three inclusions. We take a computational domainΩ = [−1, 1] × [−1, 1] m 2 and the point source is located at (0.3, 2), deliberately slightly uncentred to avoid a perfect alignment of all the components, figure 7(b). The ambient (or surrounding) medium considered is soil and, of the three inclusions, we look at two small iron mines (centered at (0.2, 0) and (0, −0.2), respectively)  and a larger plastic mine enclosing one of the small ones, see figure 7(a). The large inclusion is centered in the computational domain. We display the total wavefield computed using the TF/SF formulation, with a layer wide of 0.2, in figure 7(c), where the scattering due to the inclusions is clearly identifiable, and the associated scattered field, obtained by subtracting the incident wave to the total field, is shown in figure 7(d). We observe all contributions from the three inclusions.
Let's now assume that we record the scattered wavefield on a square receivers' array (Γ = boundary of Ω = [−0.78, 0.78] × [−0.78, 0.78], i.e. slightly inside the region where the total field was computed). Note that we consider a receivers' array with full aperture and continuous. From these measurements, we will perform TR using the TRAC method.

A quick analogy to the 1D case
Before we present the numerical results of the TRAC method with the dissipative component, we propose a quick analog analysis in 1D to determine, as a first estimate, the potential problematic values of σ 0 in terms of accuracy of the reconstruction. To do so, we consider the case of one small sound-hard inclusion located at the centre of the computational domain. To perform the TRAC method, we introduce a trial ball enclosing the inclusion, as in figure 8.
In figure 9, we plot the obtained relative L 2 -error for this 1D analog version of the 2D problem. We propose the three possible approximations from section 3, but really consider only the two on the left and centered columns. The right column presents the case of the classical TRAC and we know that the error is mostly due to an amplitude discrepancy. For the two dissipative versions of the TRAC method, we observe one peak each occurring at σ 0 = 0.034 for both. The fact that the peaks coincide is related to this specific example, as we saw in the section 3.2 that they may happen for different values of σ 0 . The main conclusion is that the dissipation of the soil σ 0 = e −3 ∼ 0.05 is beyond that peak and reconstructions may be affected.

2D dissipative TRAC method-numerical results
To analyse the effects of dissipation when performing the TRAC method, we explore different versions of the time-reversed absorbing condition, for different orders of accuracy and different shapes of the trial domain. A set of dissipation values chosen within the range used in section 4.2 is also used and we compute the relative L 2 -error for each numerical experiment.
As mentioned in section 2.3, the choice of the shape of the trial domain does not depend on the scatterer. Therefore, we can choose several configurations and compare the quality of the construction for each of them. Moreover, depending on the chosen shape, we can opt for different orders of accuracy of the ABC to be imposed on the trial domain. In the following, we consider a circular [17], an elliptic [18], a square [53] and an octagonal trial domain. The latter is a hybrid version between the circle and the square to avoid sharp corners but still uses a higher order of boundary condition.
For a chosen trial domain B, we are solving the TRAC problem to reconstruct the scattered wavefield where the operator TRAC refers to a chosen appropriate time-reversed ABC. An energy conservation only in the very specific 3D case of concentric spheres without any dissipation has been proposed in [16]. Moreover, by construction of the TRAC problem, see [17], the time-reversed boundary condition is designed to allow the in-coming wavefield from the measurements to be absorb within the trial domain, as if the scattered response of the scatterer were going back to its source. It has been observed in previous papers, namely [17,18], that TR is surprisingly robust to noisy data, particularly in the time-domain. Therefore, solving the same boundary value problem with slightly corrupted data gives a very close solution. In that sense we find some aspects of the definition of well-posed in the sense of Hadamar. In the dissipative case, we also refer to [63] for stability, as we work with bounded domains (and finite times).
For this study, we use three TRAC operators of different orders of accuracy: (i) a zeroth-order Engquist-Majda-like complex conjugated boundary condition [19], available for any shape of the trial domain, (ii) a first-order Bayliss-Turkel-like complex conjugated boundary condition [20,22], available for a circle or an ellipse, where r denotes the distance of a point of ∂B to the centre of the circle or the ellipse; (iii) a second-order Engquist-Majda-like complex conjugated boundary condition, for a square or an octagonal boundary, with ∂/∂n the normal derivative to the inner ∂B and ∂/∂τ the corresponding tangential derivative.
We give a derivation of the ABCs in 2D in appendix section 'Derivation for a planar boundary'.

A first numerical result.
We solve (20) while considering a square trial domain with a second order Engquist-Majda complex conjugated boundary condition with or without dissipation. As in 1D, we can consider three variants of the dissipative TRAC method: the first order dissipative TRAC method with operator TRAC 2 on ∂B, the zeroth order dissipative TRAC method with the same operator where σ = 0 and the classical TRAC method, where σ = 0 everywhere in (20).  figure 11. These two values are chosen on each side of the peak estimated in section 4.2, which was happening for σ 0 = 0.034. For each reconstruction, we also compute the relative L 2 -error according to: where variant refers to the different possible approximations of the TRAC method. The exponent recalls the chosen trial domain, here a square, and the order of the ABC. Figure 10 displays the reconstruction in the case of an ambient dissipation σ 0 = 0.01. We compare the three variants of the TRAC method with the expected reconstruction, which is a cropped image of the forward scattered field after complex conjugation. All images are to the scale for a better eye comparison, and the relative L 2 -errors are also recalled in the captions. From these pictures, we can make the following observations: • We do not see much difference between the first and zeroth order approximations of the TRAC reconstructions and even the corresponding L 2 -errors, 30.97% and 33.17% respectively, are close. The reconstructions cover most of the features of the expected scattered wavefield, but also looks more oscillatory where the amplitude is low, like around the corners. • The classical TRAC method seems to perform better in view of the L 2 -error, e □,2 classical (σ 0 ) = 20.55%: on the graph, we see that the scattered field is nicely reconstructed but the amplitude from the external boundary (measurements) to the trial domain does not increase exponentially, due to the missing anti-dissipative term in the Helmholtz equation. In figure 11, we show the reconstruction in the case of an ambient dissipation σ 0 = e −3 ∼ 0.05 for the same three variants of the TRAC method: • This time, the first order dissipative TRAC method gives a better reconstruction than the zeroth order dissipative TRAC method. Clearly, for larger values of σ 0 , the role of the dissipation within the time reversed boundary condition is essential. Despite having the best L 2 -error, the 1st order dissipative TRAC still shows oscillations in low amplitude areas. • As before, the classical TRAC method seems to reconstruct the main features of the scattered wavefield, but the amplitude is not amplifies, leading to a large relative L 2 -error of 70.50%. By using a logarithmic scale, the features are better displayed and the reconstruction actually seems even better than the two others.

Remark 4.2.
It is known that ABCs applied on corners may be troublesome to implement [64]. In the current state, FreeFem++ allows to implement the problem using finite elements and no special treatment has been made on the corners of the trial domain. The results do not show artefacts due to the corners, as in figure 10(d). Moreover, we performed the same tests with a circular trial domain (and first order ABC) to check our assumption about the corners and the same features appear, confirming that the corner issue has been addressed by the software directly in the variational formulation.

Comparison between different time-reversed ABCs for a given scatterer.
For the same settings of the forward experiment, we computed multiple TRAC reconstructions with different shapes of the trial domain and different orders of accuracy of the time-reversed ABC. A total of eight simulations have run for the range of dissipation values σ 0 ∈ {0, 0.0001, 0.001, 0.01, 0.014, 0.02, 0.034, 0.042, 0.05, 0.067, 0.08, 0.1}.
The relative L 2 -errors are plotted in figure 12 and a clear statement is that using a square trial domain with a second order time-reversed ABC gives the best reconstruction for this model. Table 2. 2D numerical results: list of the different scatterers involved in the numerical analysis of the TRAC method with dissipation. Small and large scatterers are rectangles of 4 cm × 10 cm and 10 cm × 60 cm, respectively, with their longer side along the x-axis.

Scatterer
Number Centre location Size Type

Comparison between different scatterers for the second order dissipative TRAC method
on a square trial domain. Since performing the TRAC method with a square trial domain and a second order ABC seems to give the best results so far, we considered different kinds of scatterers within the same medium and the same settings for the forward problem. The type of scatterers is described in table 2. Note that Scatterer (g) is the model we worked with in the previous sections, see figure 7(a). We display in the bottom row of figure 13 the curves of the L 2 -error for the seven types of scatterer listed in table 2. On the top row, we recall the error obtained as an estimation for the 1D analogy in section 4.2 for a sound-hard inclusion. Clearly, the scales are different, with much larger errors in the 2D case. However, we notice a similar behaviour of the curves with a peak occurring for σ 0 between 0.01 and 0.05 (not always quite at 0.034, though). The second main observation is that the classical TRAC method in the 2D case seems to perform better than the other with an L 2 -error always under 100%, while the other can become unlimitedly large. These curves show that the dissipative TRAC methods look less stable and give more oscillatory solution than the classical TRAC method, whose major default is the incorrectly recovered amplitude, while the other features are clearly recognisable.
As a result, the 1D analogy allows to get a rough estimate of the problematic values of σ 0 for an acceptable reconstruction. However, the classical TRAC method reveals itself as the most accurate method in term of information recovery. A more comprehensive list of tests has been performed and the observations always lead to the same result, that is, the 1st order dissipative TRAC method with a high order ABC (either on a square or octagonal trial domain) performs best overall, but the classical TRAC method recovers the features better despite an incorrect amplitude.

Conclusion, discussion and perspective
In this paper, we explored the feasibility and the limit of the TRAC method in the presence of dissipation in the model. Dissipation indeed annihilates the time reversibility of the wave equation, and an anti-dissipative term is added in the time-reversed problem, which jeopardises the stability of the solution. We then conducted an analysis in one-space dimension and frequency domain to prove that the TRAC method using an exact ABC gives an exact reconstruction of the targeted wavefield. Furthermore, since the exact ABC in the frequency domain does not lead to a differential operator in the time domain after inverse Fourier Transform, we investigated three possible variants to the exact TRAC method with different levels of approximation. Ranges of values for the dissipation have been demonstrated for which the approximated TRAC methods can still give satisfactory reconstructions. A similar analysis is conducted in two-space dimension and reveals that the reconstruction is sensitive to the order of accuracy of the chosen time-reversed absorbing condition. Moreover, if we disregard the discrepancy in the amplitude, due to the missing anti-dissipative term in the classical TRAC method, the latter proves to give the most satisfactory recovery of the main features of the reconstructed scattered wavefield.
As it stands, the current dissipative version of the TRAC method is not yet transposable to the time domain. Another level of approximation is needed to get a differential operator, nicer to implement, after inverse Fourier Transform. Another question to be explored is about heterogeneous media. In a subsequent study, we plan to extend our findings to the time domain in two-space dimensions with heterogeneous media.

Data availability statement
No new data were created or analysed in this study.

Acknowledgment
M C would like to thank the University of Auckland for granting her a Summer Research Scholarship from December 2021 to February 2022, which allowed her to discover this branch of Applied Mathematics and contribute to this paper.
which can be only achieved if C 2 = 0, i.e. only the dissipative wavefield e λσ x propagates through {x L }.
The same calculation can be done on {x R } to show that C 1 must vanish to keep only the dissipative solution e −λσ x as x goes to +∞.

In the frequency domain and two-space dimension (2D)
Let's first remind three well-known ABCs in 2D without dissipation: (1) First order Engquist-Majda boundary condition [19] on linear boundaries (also zero-th order Bayliss-Turkel boundary condition [20] on circular boundaries): where ∂/∂n denotes the external normal derivative; (2) First order Bayliss-Turkel boundary condition [20] on circular boundaries: where ∂/∂n denotes the external normal derivative and r is the radius of the circular boundary; (3) Second order Engquist-Majda boundary condition [19] on linear boundaries: where ∂/∂n denotes the external normal derivative and ∂ 2 /∂τ 2 the second tangential derivative.
Derivation for a planar boundary. For the purpose of this paper, we derive analog ABCs including the dissipation term. We proceed as in [65] for planar boundaries, since our goal is to work with second order boundary conditions on square computational domains. The dissipative Helmholtz equation reads: the dissipation coefficient being in κ ∈ C as in (2) (or in Remark 2.1). To derive the boundary condition, we first consider the solution u to (25) in the right half space x ⩾ 0 independently of the boundary condition satisfied at x = 0, as a first example. The other half-planes can be derived similarly.
which translates, after inverse Fourier Transform, into: A second order approximation is also conceivable and yields: which translates, after inverse Fourier Transform, into: Similar derivation can be done for the other three half-planes and we obtained the following absorbing conditions for straight boundaries, where n and τ denote respectively the normal and the tangential coordinates: • zeroth-order: • second-order: Remark A.1. Note that these ABCs are easily implemented in the frequency domain. However, the inverse Fourier Transform in time will again not lead to a partial differential equation and may be too costly to implement. This is the main motivation of this study. In the current form, √ κ prevents a direct translation, though. Approximations will need investigations, which we aim to study in a part II. ■ Numerical results. To conclude this section, we propose to show some numerical examples that validate the ABCs (30) and (31). The numerical experiments will be performed using the software FreeFem++ [62]. For this experiment, we choose a frequency ν equal to 1 GHz, and ω = 2πν denotes the angular frequency in the Helmholtz equation. The medium is assumed to be homogeneous and dissipative, and its physical characteristics are described using electromagnetic notations: where µ 0 and ε 0 denote the magnetic permeability and the electric permittivity in the vacuum, respectively. Electric conductivity σ will be varying for the purpose of the numerical experiment. We also introduce the wavelength, λ, computed from the speed of light c 0 in the vacuum and the frequency: We choose a squared domain with a point source at the centre: f(x, y) := e − x 2 +y 2 δ 2 , with δ = 10 −3 . At the boundary, we use either the zeroth-order or the second-order ABC (30) and (31). For the discretisation, we use P 2 -finite elements.
• Experiment 1-illustration In figure 14, we display the real part of the solution to the dissipative Helmholtz equation. We actually only show the north-east corner, and the rest can be deduced by symmetry. We performed four tests using the dissipation σ = 0.05: (a) using (30) on a squared domain of length 20λ, i.e. distance to the source of 10λ, (b) using (31) with a distance to the source of 10λ, (c) using (30) with a distance to the source of 20λ (length of squared domain of 40λ), and (d) using (31) with a distance to the source of 20λ.
As expected, the solution is decaying faster than the non-dissipative counterpart (tests performed but not displayed) and we also notice that the ABC becomes more accurate as the boundary is farther from the source term. This is especially visible at the corners.

• Experiment 2-quantification of the error
Let's take the solution performed with the second-order ABC, P 2 -finite elements, on a squared mesh of edge 40λ (distance to the source of 20λ) as our reference solution to evaluate the error due to the approximated ABC.
To compute the error, we consider the solution recorded on a vertical line well inside the computational domain x ∈ [−5λ, 5λ], at some distance from the point source. We perform the following three tests for σ = 0 (no dissipation for reference), σ = 0.05 and σ = 0.2: • the recordings are at 5 2 3 λ from the source in a computational domain of half-length 6 2 3 λ; • the recordings are at 5 2 3 λ from the source in a computational domain of half-length 10λ; • the recordings are at 8λ from the source in a computational domain of half-length 10λ; The relative L 2 -error is reported in table 3 and we can clearly see that the ABC is improved as the boundary is farther from the source. We also observe that the errors are of similar range of values as in the non-dissipative case with a better performance from the second-order ABC than the zeroth-order ABC, which indicates that the ABC is doing the right job.
We finally display the recordings (real part, imaginary part and modulus) in figure 15. Again, the second-order ABC performs better than the zeroth-order ABC (less oscillatory), and the ABCs are more accurate on larger computational domains.
In figure 16, we take the same settings and compare the results of figure 15 with the solution obtained thanks to the dissipative first-order Bayliss-Turkel ABC (32). The solution looks smoother (less oscillatory) and the amplitude is closer to the reference. This confirms that (32) may be a good update in the case of dissipation.