View source for Time reversal

{{metamodelGetDP|time_reversal}}

== Additional information ==

To run the time reversal model, open '''TR.pro''' with Gmsh.

=== The experiment and the code ===

The model presented here was used in the PhD of B. Thierry <ref name=PhDThierry>B. Thierry, ''Analyse et Simulations Numériques du Retournement Temporel et de la Diffraction Multiple'', PhD. Thesis, Université de Nancy, France, 2011, http://tel.archives-ouvertes.fr/tel-00651312/</ref> (who also wrote the code) to study the super resolution phenomenon in a deterministic medium and in the time harmonic regime.

==== Time reversal mirrors ====

The principle of time reversal is to use the reversibility of the wave equation in non dissipative media to back-propagate a signal on the source that first emitted it<ref>M. Fink and C. Prada. ''Acoustic time-reversal mirrors''. Inverse Problems, 17(1) :1761–1773, 2001. </ref>. Time Reversal Mirrors (TRM), designed by Matthias Fink and his team in the 90's, are composed of a large number of transducers which can play alternatively the role of emitters and receivers. Associated to a numerical treatment module, they are able to record a signal, time reverse it, and back propagate it into the medium. The applications are numerous, such as submarine communication or kidney stone destruction.

==== Super resolution phenomenon ====

In a homogeneous medium, the wave back-propagated by the TRM focuses the source with a resolution generally given by
$$
R = \frac{\lambda d}{a},
$$
where
* $\lambda = 2\pi/k$ is the central wavelength
* $d$ is the distance between the TRM and the source
* $a$ is the aperture of the TRM (its size for a linear mirror)
When the medium is heterogeneous, it is well-known that the back propagated wave will still focus on the point source with moreover a more accurate precision than in the homogeneous medium. In the litterature, this is generally imputed to the multi-pathing effect: the multiple scattering produces multiple paths between the sources and the mirror, which then receives more information compared to the same experiment in a homogeneous medium.

==== Time reversal experiment ====

Consider a point source $\mathbf{s}$ and a rectangular time reversal mirror $\mathcal{M}$ placed at finite distance of the source. The experiment can be split into three steps
* First, the source emits a wave that propagates throughout the medium
* Second, the TRM receives and records the wave
* Finally, it back-propagates it into the medium. The back-propagated field should then focus on the point source.
The context of the experiment is the following:
* The time reversal mirror is non intrusive (does not reflect waves), continuous and thick. This last condition is to avoid cracks in the mesh.
* Waves are assumed to be time harmonic of pulsation $\omega$: $\mathscr{U}(\mathbf{x}, t) = u(\mathbf{x})e^{-ikct}$, where $k= \omega/c$ and $c$ is the speed of sound in the vacuum.
** The time-harmonic assumption implies that the time reversal operation is nothing but a phase conjugation.

==== Context of the code ====

The simulations presented here are achieved in the time-harmonic regime and in two dimensions. The equation to be solved is then the Helmholtz equation. The heterogeneous medium is modelled as multiple circular obstacles with an inner wavelength (a contrast). The wavenumber, even if not fully random, is still randomized by the position of the disks and the constrast function. The code allows to:
* Simulate the time reversal experiment in a homogeneous medium (no obstacles)
* Study the super resolution phenomenon by adding circular scaterrers between the source and the TRM
* For both media, monochromatic waves are available (one frequency) or point source waves with a bandwidth, which permits to simulate time-dependent problems (but is highly time consuming, in terms of CPU usage).

In addition, the graphical interface allows to:
* Place the TRM and the point source
* Specify the aperture of the mirror
* Add the obstacles, randomly placed in a customizable box. 
* The contrast function can be ramdonized or be manually specified, either the same for each obstacle or not.

=== Mathematical problem (monochromatic wave) ===

==== Modelling ====

Let $M$ random disks $\Omega^-_p$ be placed between the source and the mirror, with $M$ possibly equal to $0$ (homogeneous medium). Let the three step, emission, reception/time reversal and back propagation, be mathematically modelled as follows:
* '''Emission''': let $u_E$ be the wave emitted by the point source $\mathbf{s}$. Thus, $u_E$ is the solution of
\begin{equation}\label{eq:ue}
\begin{cases}
(\Delta + k^2n^2) u_E = -\delta_{\mathbf{s}} & \text{ in } \mathbb{R}^2,\\
u_E \text{ outgoing}.
\end{cases}
\end{equation}
Operator $\Delta v := \sum_{j=1}^2 \partial^2_{x_j}v$ is the Laplacian operator, $\delta_{\mathbf{s}}$ is the Dirac distribution centered on $\mathbf{s}$ and the first equation is known as the Helmholtz equation. The quantity $k$ is the wavenumber and $n$ is the contrast function such that
$$
n(\mathbf{x})=\begin{cases}
{\rm n}_p & \text{ if }\mathbf{x}\in\Omega^-_p,\\
1 & \text{ otherwise.}
\end{cases}
$$
The outgoing condition stands for the Sommerfeld radiation condition at infinity
$$
\lim_{\|\mathbf{x}\|\to\infty} \frac{1}{\|\mathbf{x}\|^{1/2}}\left(\nabla u \cdot \frac{\mathbf{x}}{\|\mathbf{x}\|} - iku\right) = 0.
$$
The unique solution of problem \ref{eq:ue} is nothing but the Green function $G_k^n(\cdot, \mathbf{s})$ centered on $\mathbf{s}$.
$$
\forall \mathbf{x}\neq\mathbf{s},\qquad u_E(\mathbf{x})= G_k^n(\mathbf{x},\mathbf{s}).
$$
In the homogeneous case ($n=1$), an analytic expression of the emitted wave $u_E^0$ is known since the Green function, denoted by $G_k^0(\cdot, \mathbf{s})$, is given by:
$$
\forall \mathbf{x}\neq\mathbf{s},\qquad G_k^0(\mathbf{x},\mathbf{s}) = \frac{i}{4}H^{(1)}_0(k\|\mathbf{x}-\mathbf{s}\|),
$$
where $H^{(1)}_0$ is the zeroth order Hankel function of the first kind.
In a heterogeneous medium and except for particular case, the Green function is unknown and must be approximated numerically.
* '''Reception''': the TRM $\mathcal{M}$ measures the emitted field ($u_E|_{\mathcal{M}}$) and time-reversed it. As highlighted above, this operation is equivalent to a phase conjugation: $\overline{u_E}|\mathcal{M}$.
* Finally, the mirror '''back-propagates''' the signal. This last step is modelled as new wave $u_{back}$, solution of the following problem
\begin{equation}\label{eq:uback}
\begin{cases}
(\Delta + k^2) u_{back} = -\overline{u_E}|_{\mathcal{M}} = -\overline{G_k^n(\cdot,\mathbf{s})}|_{\mathcal{M}} & \text{ in } \mathbb{R}^2,\\
u_{back} \text{ outgoing}.
\end{cases}
\end{equation}
The expression of $u_{back}$ is hence given by 
$$
\forall \mathbf{y}\in \mathbb{R}^2\setminus\overline{\mathcal{M}}, \qquad u_{back}(\mathbf{y}) = \int_{\mathcal{M}} G_k^0(\mathbf{x},\mathbf{y}) \overline{u_E(\mathbf{x})}\; {\rm d}\mathbf{x}.
$$
Replacing $u_E$ by the Green function leads to a convolution product between Green functions:
$$
\forall \mathbf{y}\in \mathbb{R}^2\setminus\overline{\mathcal{M}}, \qquad u_{back}(\mathbf{y}) = \int_{\mathcal{M}} G_k^n(\mathbf{x},\mathbf{y}) \overline{G_k^n(\mathbf{x},\mathbf{s})}\; {\rm d}\mathbf{x}.
$$
For a homogeneous medium, as the Green function is known explicitely, so is the back propagated wave. This function can be plot directly with, ''e.g.'' a Python script, Matlab, Scilab, etc. With GMSH/GetDP, $u_{back}$ will however be computed by solving problem \ref{eq:uback}, but using the value of the Green function. For the heterogenous case and as the Green function is unknown, the emission and the back propagation steps are both computed numerically.

==== Theoretical results (homogenous medium only) ====

Propositions 3.8 and A.7 (appendix) of <ref name=PhDThierry />, shows the focusing property of the back propagated field for a thin and a thick TRM in 2 dimensions. As only a thick TRM is considered here, then only this result is summaries in this page. First, let be introduced the cone $\mathcal{C}$ of apex $\mathbf{s}$ and basis the mirror (see figure on right). Then, when the wavenumber $k$ tends to infinity, the back-propagate field $u_{back}$ focuses on the point source $\mathbf{s}$ in the sens that
* if $\mathbf{y} \in \mathbb{R}^2\setminus\overline{\mathcal{C}}$, then $|u_{back}(\mathbf{y})| = O(k^{-2})$,
* if $\mathbf{y} \in \mathcal{C}$, then $ |u_{back}(\mathbf{y})| = O(k^{-3/2})$,
* if $\mathbf{y} = \mathbf{s}$, then $ |u_{back}(\mathbf{y})| = O(k^{-1})$.

=== Mathematical problem (Band of frequencies) ===

This example is also the opportunity to show results using a band of frequencies. In that case, $k\in[k_0 - \Delta_k, k_0+\Delta_k]$, with, ''e.g.'' $k_0 =10$ and $\Delta_k = k_0/5$. 

==== Emitted wave ====

The emitted wave $\mathscr{U}_E$ is then obtained by integrating over all the frequencies:
$$
\mathscr{U}_E(\mathbf{x}, t) = \int_{k=k_0-\Delta_k}^{k_0+\Delta_k} u_E^k(\mathbf{x})e^{-i k ct}\;{\rm d}k
= \int_{k=k_0-\Delta_k}^{k_0+\Delta_k} G_k^n(\mathbf{x})e^{-i k ct}\;{\rm d}k,
$$
where $c$ is the velocity in the medium.

==== Back propagated wave $\mathscr{U}_{back}$ by the TRM ====

Recall that the expression of the back propagated field $u_{back}(\cdot\,,k)$ of frequency $k$, in the time harmonic regime, is given by,
$$
\forall k\in[k_0 - \Delta_{k},k_0 + \Delta_{k}],\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad  u_{back}(\mathbf{y},k) = \int_{\mathcal{M}}\overline{G_{k}^n(\mathbf{x},\mathbf{s})}G_{k}^n(\mathbf{x},\mathbf{y})\;{
\rm d}\mathbf{x}.
$$
The back propaged wave $\mathscr{U}_{back}$ is hence given by
$$
\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad  \mathscr{U}_{back}(\mathbf{y},t) = \int_{k=k_0-\Delta_k}^{k_0+\Delta_k}\int_{\mathcal{M}}\overline{G_{k}^n(\mathbf{x},\mathbf{s})}G_{k}^n(\mathbf{x},\mathbf{y})\;{
\rm d}\mathbf{x} e^{-ikct}{\rm d}k.
$$

==== $3$D and homogeneous medium: theoretical results ====

In the $3$D and homogeneous case, the expression of the emitted field reduces to
$$
\forall t \in\mathbb{R}, \forall\mathbf{x}\in\mathbb{R}^{3}\setminus\{\mathbf{s}\}, \qquad \mathscr{U}_E(\mathbf{x},t) = \frac{\Delta_k{\rm sinc}\left[\Delta_k\left(\|\mathbf{x}-\mathbf{s}\|-ct\right)\right]  e^{ik_0\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)}}{2\pi\|\mathbf{x}-\mathbf{s}\|},
$$
where ${\rm sinc}(x) = \sin(x)/x$. Remark first that, due to the quantity $e^{ik_0\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)}$, the wave $\mathscr{U}_E$ is propagating in a radial direction with respect to $\mathbf{s}$. The modulus of $\mathscr{U}_E(\mathbf{x},t)$ especially depends on the cardinal sine term, as:
$$
\forall t \in\mathbb{R}, \forall\mathbf{x}\in\mathbb{R}^{3}\setminus\{\mathbf{s}\}, \qquad \left|\mathscr{U}_E(\mathbf{x},t)\right| = \left|\frac{\Delta_k{\rm sinc}\left[\Delta_k\left(\|\mathbf{x}-\mathbf{s}\|-ct\right)\right] }{2\pi\|\mathbf{x}-\mathbf{s}\|}\right|.
$$
In particular, for a point $\mathbf{x}$ close to $\mathbf{s}$, the modulus of $\mathscr{U}_E(\mathbf{x},t)$ is reached at $t = 0$, when the value of the cardinal sine is $1$. This shows the presence of an impulsion at time $t=0$ at the source point. On the other hand, on a point $\mathbf{x}$ located at a distance $D$ from the source point, the modulus of $\mathscr{U}_E(\mathbf{x},t)$ is given by
$$
\forall t \in\mathbb{R}, \qquad \left|\mathscr{U}_E(\mathbf{x},t)\right| = \left|\frac{\Delta_k{\rm sinc}\left[\Delta_k\left(D-ct\right)\right]}{2\pi D}\right|.
$$
This quantity reaches its maximal value when $t = -D/c$. In other words, the perturbation created at time $t=0$ reaches the distance $D$ and thus the MRT, at time $t = D/c$.

On the other hand, the back-propagated field reduces to:
\begin{equation}\label{eqRT:Ureyt}
\mathscr{U}_{back}(\mathbf{y},t) =   2{\rm d}k\int_{\mathcal{M}}  
\frac{{\rm sinc}\left[{\rm d}k\left(\left[\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right] - ct\right)\right] 
e^{ik_0\left(\left[\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right]-ct\right)}}{\|\mathbf{x}-\mathbf{s}\|\|\mathbf{x}-\mathbf{y}\|}
\;{\rm d}\mathbf{x}.
\end{equation}
It should be pointed out that on the point $\mathbf{y}=\mathbf{s}$, the cardinal sine term vanishes and the above expression reduces to
$$
\mathscr{U}_{back}(\mathbf{s},t) =   2{\rm d}k \int_{\mathcal{M}}  
\frac{{\rm sinc}\left(-{\rm d}k ct\right)}{\|\mathbf{x}-\mathbf{s}\|^{2}}
 \; {\rm d}\mathbf{x} = 2{\rm d}k{\rm sinc}\left({\rm d}k ct\right) \int_{\mathcal{M}}  
\frac{{\rm d}\mathbf{x}}{\|\mathbf{x}-\mathbf{s}\|^{2}}. 
$$
Hence, when $\mathbf{y} = \mathbf{s}$, the modulus of $\mathscr{U}_{back}(\mathbf{s},t)$ reaches its maximal value at time $t = 0$. On the other hand and on a point $\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}}$, when $t=0$ the equality \ref{eqRT:Ureyt} becomes
$$
\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad \mathscr{U}_{back}(\mathbf{y},t = 0) =   2{\rm d}k\int_{\mathcal{M}}  
\frac{{\rm sinc}\left[{\rm d}k\left(\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right)\right] 
e^{ik_0\left(\left[\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right]-ct\right)}}{\|\mathbf{x}-\mathbf{s}\|\|\mathbf{x}-\mathbf{y}\|}
\;{\rm d}\mathbf{x}.
$$
On $\mathbf{y} = \mathbf{s}$, the quantity $\mathscr{U}_{back}(\mathbf{y},t = 0)$ does not depend on $k$ anymore and reads as
$$
\mathscr{U}_{back}(\mathbf{s},t = 0) =   2{\rm d}k\int_{\mathcal{M}}  
\frac{{\rm d}\mathbf{x}}{\|\mathbf{x}-\mathbf{s}\|^{2}}.
$$
For every other point $\mathbf{y}$ (up to symetric points of $\mathbf{s}$ with respect to the TRM), the stationary phase theorem can be applied, as in the time harmonic regime. Hence, a decrease of $\left|\mathscr{U}_{back}(\mathbf{y},t = 0)\right|$ with respect to $k_0$ can be proven. This decrease is moreover stronger du to the presence of the cardinal sine term. As a consequence, at time $t = 0$, the wave $\mathscr{U}_{back}$ focuses on the point source $\mathbf{s}$.

Finally, remark that at $t=0$, the general expression of $\mathscr{U}_{back}$, becomes
$$
\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad \mathscr{U}_{back}(\mathbf{y},t = 0) = \int_{k_0-{\rm d}k}^{k_0+{\rm d}k} \mathscr{U}_{back}(\mathbf{y},k) \;{\rm d} k.
$$
In other words, the wave $\mathscr{U}_{back}$ on $\mathbf{y}$ and at time $t=0$ is obtained by the sum of every time-harmonic fields $\mathscr{U}_{back}(\cdot\,,k)$, for $k\in[k_0-{\rm d}k,k_0+{\rm d}k]$. The quantity $\mathscr{U}_{back}(\mathbf{y},t = 0)$ can then be seen as an average over the frequencies $k$ of the different fields $u_{back}(\cdot\,,k)$.

=== Numerical resolution ===

==== Approximation with a PML ====
To be solved using finite element method, problem \ref{eq:uback} is truncated using a Perfectly Matched Layer (PML), as in [[Acoustic_scattering|the scattering problem with a PML]]. 
<!--First, let us build the domain $\Omega_{e}$ on which an approximation of the back propagated field will be computed. A rectangular shape for $\Omega_{e}$ is natural in that problem:
$$
\begin{array}{rcl}
\Omega_e &=& \left(\left]\alpha_1,\beta_{1}\right[\times\left]\alpha_{2},\beta_{2}\right[\right) \ \setminus\  \overline{\Omega^-} \\
                 &=& \left(\left]\alpha_1,\beta_{1}\right[\times\left]\alpha_{2},\beta_{2}\right[\right) \ \bigcap\  {\Omega^+},
\end{array}
$$
where $\alpha_1,b_{1},\alpha_{2}$ and $\beta_{2}$ are real constant such that $\alpha_{j}<\beta_{j}$ for $j=1,2$. Obviously, this set is sufficiently large to contain the source and the time reversal mirror.

Around the rectangle $\Omega_e$ is built the absorbing layer $\Omega_{PML}$ defined by
$$
\Omega_{PML} = \left(\;\left]\alpha_1 - \delta^{PML}_{1}, \beta_{1}+\delta^{PML}_{1}\right[\times\left]\alpha_{2} - \delta^{PML}_{2},\beta_{2} + \delta^{PML}_{2}\right[ \;\right)\setminus\overline{\Omega_e},
$$
where $\delta^{PML}_{1}$ and $\delta^{PML}_{2}$ are the thickness of the PML in the $x$ et $y$ direction, respectively. Let $\Sigma_{int}$ be the boundary separating $\Omega_e$ from $\Omega_{PML}$ and $\Sigma_{PML} = \partial\Omega_{PML}\setminus\Sigma_{int}$, the boundary which truncated the PML.
Finally, let us introduce $\Omega_{T} = \overline{\Omega_e}\bigcup\Omega_{PML}$ the total computation domain and.  In practice, the thickness $\delta^{PML}_{1}$ and $\delta^{PML}_{2}$ are very small, e.g. : 10 elements.
[[File:TRFreeSpace_PML.png|center|600px|frameless|Domains and parameters.|alt=The different domains and parameters.]]

The approximation $u^{PML}_{back}$ of $u_{back}$ in $\Omega_e$ is solution of a partial differential equation in $\Omega_T$ and satisfies a boundary condition on $\Sigma_{PML}$. The field $u^{PML}_{back}$ is dissipated through it passage in $\Omega_{PML}$. This absorption is obtained by modifying the Helmholtz equation. However, let us first focus on the boundary condition. If the outgoing waves are sufficiently damped by the PML then the amplitude will be close to zero in a neighborhood of $\Sigma_{PML}$. Thus, the role of the boundary condition will be unimportant <ref>P. Petropoulos, ''On the Termination of the Perfectly Matched Layer with Local Absorbing Boundary Conditions'', Journal of Computational Physics, 1998</ref>. That is why is proposed to set a homogeneous Dirichlet boundary condition on $\Sigma_{PML}$. However, an absorbing boundary condition could be employed to slightly enhance the accuracy of the result.

Let us now focus on the absorption. As said previously, this absorption is obtained by modifying the Helmholtz equation inside $\Omega_{PML}$. More precisely, the problem solved by $u^{PML}_{back}$ is the following
\begin{equation}\label{eqRT:HelmholtzPML}
\begin{array}{r c l l}
\displaystyle{\partial_{x_{1}}\left( \frac{S_{x_{2}}}{S_{x_{1}}}\partial_{x_{1}}u^{PML}_{back} \right) +
\partial_{x_{2}}\left( \frac{S_{x_{1}}}{S_{x_{2}}}\partial_{x_{2}}u^{PML}_{back} \right)  + 
k^{2} S_{x_{1}}S_{x_{2}}u^{PML}_{back}} &=& f & (\Omega_{T})\\
u^{PML}_{back} &=& 0 & (\Sigma_{PML})
\end{array}
\end{equation}
with
$$
\forall\mathbf{x}=(x_1,x_2)\in\overline{\Omega_{T}},\qquad\begin{cases}
\displaystyle{S_{x_{1}}(\mathbf{x}) = k + i\sigma_{x_{1}}(x_1)}, & \\
\displaystyle{S_{x_{2}}(\mathbf{x}) = k + i\sigma_{x_{2}}(x_2)}. &
\end{cases}
$$
The functions $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are called damping functions. When $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are equal to zero, then equation \ref{eqRT:HelmholtzPML} reduce to (classical) Helmholtz equation, that is why in general, these functions vanish on $\Omega_{e}$. On the other hand, when the damping functions are positive, they act as a dissipative term in the equation. Remark that, by introducing the following matrix
$$
D = \left[ 
\begin{array}{c c}
\frac{S_{x_{2}}}{S_{x_{1}}} & 0\\
0 & \frac{S_{x_{2}}}{S_{x_{1}}}
\end{array}
 \right],
$$
then problem \ref{eqRT:HelmholtzPML} can be rewritten as
\begin{equation}\label{eqRT:eqPML}
\left\{\begin{array}{r c l l}
\displaystyle{{\rm div}(D\nabla u^{PML}_{back}) + S_{x_{1}}S_{x_{2}}u^{PML}_{back}}& =& f & (\Omega_{T})\\
\displaystyle{u^{PML}_{back}} &=& 0 & (\Sigma_{PML})
\end{array}\right.
\end{equation}
The choice of the damping functions $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ obviously influces the dissipative power of the PML. Here, these functions are simply linear and, moreover, continuous on $\Omega_{T}$, in order to avoid any non desirated reflection. As they vanish in $\Omega_e$, these functions stays equal to zero on $\Sigma_{int}$ and reach their maximum on $\Sigma_{PML}$, which will be set to $100$. This choice is purely given by the experiment and is not a general coefficient. In fact, for another problem, this parameter must certainly be tuned.

To summarize, the two damping functions $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are given by:
$$
\forall j=1,2, \forall \mathbf{x}=(x_1,x_2) \in\overline{\Omega_{T}}, \qquad \sigma_{x_{j}}(\mathbf{x}) = 
\begin{cases}
0 & \text{ if } \alpha_{j}< x_{j}< \beta_{j},\\
\displaystyle{\frac{100\left|x_{j}-\beta_{j}\right|}{\delta^{PML}_{j}}}&\text{ if } \beta_{j} \leq x_{j} \leq \beta_{j} + \delta^{PML}_{j}, \\
\displaystyle{\frac{100\left|x_{j}-\alpha_{j}\right|}{\delta^{PML}_{j}}}&\text{ if } \alpha_{j}-\delta^{PML}_{j} \leq x_{j} \leq \alpha_{j}.
\end{cases}
$$
For $j=1,2$, function $\sigma_{x_{j}}$ has the following appearance, with respect to $x_{j}$ for a fix $x_{i}$ with $\alpha_{i} - \delta^{PML}_{i}\leq x_{i} \leq \beta_{i} +\delta^{PML}_{i}$, for $i,j = 1,2$ with $i\neq j$. 
[[File:pml_illustration_helmholtz_sigma.png|center|600px|frameless|Damping function.|alt=Damping function.]]


Inside $\Omega_e$, the damping functions vanish whereas inside $\Omega_{PML}$, at least one of them do not vanish. Note that in the "corners" of $\Omega_{PML}$, both damping function are note equal to zero.

Finally, let us just remark that the efficiency of the PML can be improved by using another type of damping functions, as quadratic functions or even functions with an infinite integral on $\Omega_{PML}$, as the one proposed by Bermudez et. al <ref>A. Bermúdez, L. Hervella-Nieto, A. Prieto and R. Rodríguez, ''An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems'', J. Comput. Phys., 2007</ref>.-->

==== Weak formulation ====

To simplify, the approximation of the back propagated field in $\Omega_e$ is still denoted by $u_{back}$ instead of $u^{PML}_{back}$. Function $u_{back}$ will be sough in the Sobolev space $H^1(\Omega_T)$and the weak formulation of problem \ref{eq:uback} reads as
$$
\left\{
\begin{array}{l}
\displaystyle{\text{Find } u \in H^1(\Omega_T) \text{ such that}}\\
\displaystyle{\forall u'\in H^1(\Omega_T),\qquad 
\int_{\Omega}D \nabla u \cdot\nabla u'\;{\rm d}\Omega 
- \int_{\Omega}k^2 S_{x_1} S_{x_2} uu'\;{\rm d}\Omega 
+ \int_{\mathcal{M}} \overline{u_E} u'\;{\rm d}\mathcal{M} 
=0,}
\end{array}\right.
$$
where $u_E$ is here equal to the Green function: $u_E = G(\cdot,\mathbf{s})$.

=== Results ===

==== Homogeneous medium ====

The three following pictures depict respectively the geometry, the real part of $u_{back}$ and the absolute value of $u_{back}$, with the following parameters ($\lambda_{geo} = 2\pi/k_{geo}$):
* $k_{geo} = 10$
* Distance between source and mirror: $20\lambda_{geo}$
* Aperture of the TRM: $15\lambda_{geo}$
* Thickness of the TRM: $0.2\lambda_{geo}$
* $k_{dis} = 15$
* $k = k_{geo} = 10$

For a monochromatic wave of frequency $k_0=k_{geo}$, which corresponds to a harmonic source, the results are:

[[File:TRFreeSpace_geo.png|250px|frameless|alt=Geometry|Geometry and ONELAB interface for a homogeneous medium.]]
[[File:TRFreeSpace_Uback.jpg|250px|frameless|alt=Real part of the back propagated field in a heterohomogeneous medium (monochromatic wave|Real part of the back propagated field in a heterohomogeneous medium (monochromatic wave.]]
[[File:TRFreeSpace_Uback_abs.jpg|250px|frameless|alt=Absolute value of the back propagated field in a heterohomogeneous medium (monochromatic wave)|Absolute value of the back propagated field in a heterohomogeneous medium (monochromatic wave.]]

Now, using the same geometry, a polychromatic wave is used as emitted wave instead of a pulse. The frequency now lies between $[k_0-\Delta_k, k_0+\Delta_k]$ with $k_0 = k_{geo} = 10$ and $\Delta_k=k_0/5 = 2$. The segment $[8,12]$ is sampled by $64$ frequencies and here is the result (on left the real part of the back propagated field, and on the right, its absolute value):

[[File:TR_bandwidthRE.png|250px|frameless|alt=Real part of the back propagated field in a homogeneous medium (polychromatic wave)|Real part of the back propagated field in a homogeneous medium (polychromatic wave).]][[File:TR_bandwidth.png|250px|frameless|alt=Absolute value of the back propagated field in a homogeneous medium (polychromatic wave)|Absolute value of the back propagated field in a homogeneous medium (polychromatic wave).]]

==== Heterogeneous medium ====

The same geometry is used again but with arround 350 circular obstacles placed between the source and the mirror. They are all of radius $\lambda_{geo}/5$ and placed in a box between the source and the TRM. The results are the following, first for a monochromatic wave:

[[File:TRCLUTTER_onelab.png|250px|frameless|alt=Geometry and ONELAB interface for heterogeneous medium|Geometry and ONELAB interface for heterogeneous medium.]]
[[File:TRCLUTTER_monoRE.png|250px|frameless|alt=Real part of the back propagated field|Real part of the back propagated field in a heterogeneous medium (monochromatic wave).]]
[[File:TRCLUTTER_mono.png|250px|frameless|alt=Absolute value of the back propagated field|Absolute value of the back propagated field in a heterogeneous medium (monochromatic wave).]]

And second, for a polychromatic wave (same parameters as for the homogeneous medium). The focal spot is clearly sharper for the heterogeneous medium  (on the left the real part of the back propagated field, and on the right, its absolute value):

[[File:TRCLUTTER_bandwidthRE.png|250px|frameless|alt=Real part of the back propagated field in a heterogeneous medium (polychromatic wave)|Real part of the back propagated field in a heterogeneous medium (polychromatic wave).]]
[[File:TRCLUTTER_bandwidth.png|250px|frameless|alt=Real part of the back propagated field in a heterogeneous medium (polychromatic wave)|Real part of the back propagated field in a heterogeneous medium (polychromatic wave).]]

== References ==

<references />

{{metamodelfooter|time_reversal}}