Éphane St

Descombes

S T Éphane Lanteri
email: stephane.lanteri@inria.fr

Ludovic Moya
email: ludovic.moya@inria.fr

LOCALLY IMPLICIT DISCONTINUOUS GALERKIN TIME DOMAIN METHOD FOR ELECTROMAGNETIC WAVE PROPAGATION IN DISPERSIVE MEDIA APPLIED TO NUMERICAL DOSIMETRY IN BIOLOGICAL TISSUES

Keywords: time-domain Maxwell's equations, dispersive media, Debye, discontinuous Galerkin method, locally implicit scheme AMS subject classifications. 35Q61, 65L20

HAL is a multi-disciplinary

Locally implicit discontinuous Galerkin time domain method for electromagnetic wave propagation in dispersive media applied to numerical dosimetry in biological tissues

Stéphane Descombes, Stéphane Lanteri, Ludovic Moya

To cite this version:

Stéphane Descombes, Stéphane Lanteri, Ludovic Moya. Locally implicit discontinuous Galerkin time domain method for electromagnetic wave propagation in dispersive media applied to numerical dosimetry in biological tissues. SIAM Journal on Scientific Computing, 2016, 38 [START_REF] Converse | Ultrawide-band microwave space-time beamforming for hyperthermia treatment of breast cancer: a computational feasibility study[END_REF], pp.A2611-A2633. [START_REF] Dolean | Locally implicit Discontinuous Galerkin method for time domain electromagnetics[END_REF].1137/15M1010282. hal-01133694v2

1. Introduction. This article is concerned with the numerical simulation of electromagnetic wave propagation in dispersive media. These are materials in which either or both of the electromagnetic material parameters ε and µ are functions of frequency. Note that the conductivity σ may also be a function of frequency, but its effect can be incorporated in the complex permittivity. In reality, all materials have frequency-dependent ε and µ, but many materials can be approximated as frequency-independent over a frequency band of interest, simplifying their analysis and simulation. Here, we will focus on the much more common case of frequencydependent permittivity. A lot of practical electromagnetic wave propagation problems involve such propagation media, such as those involving the interaction of an electromagnetic wave with biological tissues. The numerical modeling of the propagation of electromagnetic waves through human tissues is at the heart of many biomedical applications such as the microwave imaging of cancer tumors or the treatment of the latter by hyperthermia. For example, microwave imaging for breast cancer detection is expected to be safe for the patient and has the potential to detect very small cancerous tumors in the breast [START_REF] Fear | Microwaves for breast cancer detection[END_REF][START_REF] Klemm | Towards contrast enhanced breast imaging using ultra-wideband microwave radar system[END_REF][START_REF] Schea | Contrast-enhanced microwave imaging of breast tumors: a computational study using 3D realistic numerical phantoms[END_REF]. Beside, the definition of microwavebased hyperthermia as an immunotherapy strategy for cancer can also be cited [START_REF] Converse | Ultrawide-band microwave space-time beamforming for hyperthermia treatment of breast cancer: a computational feasibility study[END_REF][START_REF] Karanasiou | Development and laboratory testing of a noninvasive intracranial focused hyperthermia system[END_REF]. The electroporation technique can also be an application, which consists of applying nanopulses to the tissues, enabling only intracellular membranes to be affected, and then opening the route to therapeutic strategies such as electrochemotherapy or gene transfer [START_REF] Miklavčič | A validated model of in vivo electric field distribution in tissues for electrochemotherapy and for DNA electrotransfer for gene therapy[END_REF][START_REF] Šel | Sequential Finite Element model of tissue electropermeabilization[END_REF][START_REF] Serša | Application of electroporation in electrochemotherapy of tumors, in Electroporation based Technologies and Treatment[END_REF][START_REF] Sukharev | Electroporation and electrophoretic DNA transfer into cells: the effect of DNA interaction with electropores[END_REF][START_REF] Tsong | Electroporation of cell membranes[END_REF]. Because for all these biomedical applications experimental modeling is almost impossible, computer simulation is the approach of choice for understanding the underlying physical phenomena. In order to do so, one needs accurate and efficient numerical modeling techniques, able to deal with the complex issues characterizing the associated propagation problems.

Numerical simulation of wave propagation in dispersive media started in early 1990's in the framework of Finite Difference Time Domain (FDTD) methods, for details and references see e.g. [START_REF] Inan | Numerical Electromagnetics: The FDTD Method[END_REF] or [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]. Finite Element Time Domain (FETD) methods were not explored until 2001 [START_REF] Jiao | Time-domain Finite-Element modeling of dispersive media[END_REF] and Discontinuous Galerkin Time Domain (DGTD) methods for solving Maxwell's equations in dispersive media have been considered more recently. DGTD methods share almost all the advantages of FETD methods (large spectrum of applications, complex geometries, etc) and FVTD methods (ability to capture discontinuous solutions). The DGTD method has other nice properties which explain the renewed interest it gains in various domains in scientific computing as witnessed by books or special issues of journals dedicated to this method [START_REF] Cockburn | Discontinuous Galerkin Methods[END_REF][START_REF]Special issue on Discontinuous Galerkin methods[END_REF][START_REF] Dawson | Special issue on Discontinuous Galerkin methods[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF]. Concerning Maxwell's equations in the non-dispersive case, DGTD discretizations have been introduced in [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws. V. Multidimensional systems[END_REF] and [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF], Concerning the dispersive case, in [START_REF] Lu | Discontinuous Galerkin methods for dispersive and lossy Maxwell's equations and PML boundary conditions[END_REF][START_REF] Lu | Discontinuous Galerkin Time domain methods for GPR simulation in dispersive media[END_REF], a DGTD method, which uses piecewise high-order polynomials for spatial discretization and a Runge-Kutta method for time integration, is derived for linear dispersive media of Debye type, the treatment of the dispersive character relying on an Auxiliary Differential Equation (ADE) approach. Numerical results for two-dimensional problems are given and no proof of convergence is provided. In [START_REF] Huang | Interior penalty Discontinuous Galerkin method for Maxwell's equations in cold plasma[END_REF][START_REF] Huang | Interior penalty Discontinuous Galerkin methods for Maxwell's equations in dispersive media[END_REF], a priori error estimates are proved for the second-order formulation of Maxwell's equations coupled to dispersive models discretized by an interior penalty DGTD formulation. Some two-dimensional numerical tests are included for supporting the analysis. In [START_REF] Wang | Error analysis of a Discontinuous Galerkin method for Maxwell's equations in dispersive media[END_REF], different dispersive media are treated, considering a locally divergence-free DG method. The scheme is written and studied in its semi-discretized version, while the fully discrete scheme is described but not analyzed. Finally, in [START_REF] Lanteri | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF], which deals with the Debye model, the same DG method than the one considered in the present article is used with a second-order leap-frog scheme for time integration. Stability estimates are derived through energy conservation and convergence is proved for both the semi-discrete and the fully discrete scheme. A two-dimensional artificial numerical problem is presented to validate the theoretical findings. All these studies rely on fully explicit or fully implicit time schemes. Here we propose and analyze a DG method coupled with an efficient time integration method for dealing with grid induced stiffness when using non-uniform (locally refined) meshes and use these methods to study the interaction of electromagnetic waves with biological tissues. Stability and convergence analysis of the resulting scheme are studied, and three-dimensional problems are presented.This scheme belongs to the class of implicit-explicit (IMEX) scheme providing a single time step for stability, independent of the fine grid. Other methods also exist based on explicit local time-stepping techniques (see [START_REF] Diaz | Energy conserving explicit local time stepping for second-order wave equations[END_REF][START_REF] Seny | Multirate time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows[END_REF][START_REF] Gödel | Gpu accelerated adamsbashforth multirate discontinuous galerkin fem simulation of high-frequency electromagnetic fields[END_REF] and [START_REF] Grote | Explicit local time-stepping methods for Maxwell's equations[END_REF] specially for Maxwell's equations), which consist to use smaller time steps, given by a local stability criterion, precisely where the smallest elements are located to overcome the bottleneck caused by local mesh refinement in explicit time integrators. This paper is organized as follows. Section 2 presents the formulation of Maxwell's equations for Debye dispersive media. The Debye model is most often used to model electromagnetic wave interactions with water-based substances, such as biological materials. In particular, biological tissues are well represented by multi-pole Debye models. The dispersive character will be taken into account via an ADE which relates the electric polarization to the electric field. In Section 3 Maxwell's equations in dispersive media are discretized according to a DG formulation and we adapt the locally implicit time integration method from [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] with the additional ADE, The stability of the resulting locally implicit DGTD method is analyzed via an energy approach that provides a rigorous stability criterion. We also state a convergence result to prove that the locally implicit DGTD method for Maxwell's equations in dispersive media is second-order accurate. Finally, in section 4, we present some numerical results for three-dimensional problems.

2. The continous problem : Maxwell's equations in dispersive media. In this section we recall the initial and boundary value problem used to model the propagation of an electromagnetic wave in a dispersive medium. The physical behaviour of the medium is assumed to be described by the Debye model modeled by a combination of the first order Maxwell's equations and an ADE representation.

First, we recall Maxwell's equations which govern the electric field E and the magnetic field H in matter. Let Ω be a bounded smooth convex domain, let T > 0, Maxwell's equations in [0, T ] are given by

                     ∂ B ∂t = -curl E , ∂ D ∂t = curl H -J c -J s , div B = 0, div D = 0, (2.1)
where B and D are the magnetic and electric flux densities, J c represents the conduction current density and J s a given source current density. Introducing σ, the conductivity, Ohm's law governs the definition of the conduction current density, i.e. J c = σ E. Moreover, we will assume that J s = 0. Finally we will consider dielectric media for which magnetic effects are negligible. Introducing µ 0 the magnetic permittivity in vacuum, this yields the relation B = µ 0 H. For linear dispersive materials, the constitutive relationship between the electric flux density and electric field phasors, valid at each frequency ω and at each point in space x, can be expressed in terms of ε 0 , the electric permittivity in vacuum, and ε r the complex relative permittivity, by the following relation

D ( x, ω) = ε 0 ε r ( x, ω) E ( x, ω) (2.2)
where D and E denote the Fourier transform of D and E. Introducing P the electric polarization phasor, χ e the frequency-domain electric susceptibility of the material and ε ∞ the infinite frequency relative permittivity, equation (2.2) is often written as

D ( x, ω) = ε 0 ε ∞ ( x) E ( x, ω) + P ( x, ω) , = ε 0 ε ∞ ( x) E ( x, ω) + ε 0 χ e ( x, ω) E ( x, ω) , (2.3)
The use of a frequency-dependent electric permittivity is linked to the fact that the electric flux density D in a dielectric resulting from an applied alternating electric field E has a different phase with respect to E. This behavior is due to the inertia of the polarization P which, when frequency becomes high enough, cannot follow the rapid variations of the field, giving rise to a relaxation (i.e., a momentary delay or lag) of the permittivity. Due to the product of two frequency-dependent quantities in (2.3), the time-domain equivalent is given by the following convolution integral

D ( x, t) = ε 0 ε ∞ ( x) E ( x, t) + P ( x, t) , = ε 0 ε ∞ ( x) E ( x, t) + ε 0 t 0 E ( x, t -s) χ e ( x, s) ds.
(2.4)

In order to avoid the direct numerical treatment of the convolution integral in (2.4), approaches based on a recursive calculation of the convolution integral (Recursive Convolution Method -RCM) on one hand, and on the integration of an auxiliary differential equation (ADE) on the other hand, have been developed for implementing (2.4), see e.g. [START_REF] Luebbers | A frequencydependent Finite-difference Time-domain formulation formulation dispersive materials[END_REF][START_REF] Luebbers | Fdtd for nth-order dispersive media[END_REF][START_REF] Kelley | Piecewise linear recursive convolution for dispersive media using fdtd[END_REF] and [START_REF] Kashiwa | A treatment by fdtd method of dispersive characteristics associated with electronic polarization[END_REF][START_REF] Joseph | Direct time integration of Maxwell's equations in linear dispersive media wi th absorption for scattering and propagation of femtosecond electromagnetic pulses[END_REF], respectively. In this article we assume that the medium is a single-pole Debye type dispersive medium and concentrate on the ADE technique. Debye media are characterized by a complex valued, frequencydomain susceptibility function χ e (ω) that has one or more real poles at separate frequencies. For a single-pole Debye medium, using the notation j = √ -1, we have

χ e (ω) = ε s -ε ∞ 1 + jωτ = ∆ε 1 + jωτ ,
with ε s , called the static relative permittivity, the permittivity at zero frequency (ε s > ε ∞ ) and τ the Debye relaxation time constant, characteristic of the material. With this particular functional form of the susceptibility and (2.3) we obtain

P ( x, ω) = ε 0 ∆ε 1 + jωτ E ( x, ω) ,
and finally by inverse Fourier transformation

P + τ ∂ P ∂t = ε 0 (ε s -ε ∞ ) E. (2.5)
We can now state Maxwell's equations in a Debye dispersive medium. We infer from (2.1) and (2.5) that the magnetic field H, the electric field E and the electric polarization P satisfy the following system of equations in [0, T ]

                 µ 0 ∂ H ∂t = -curl E , ε 0 ε ∞ ∂ E ∂t = curl H - ε 0 (ε s -ε ∞ ) τ E -σ E + 1 τ P , ∂ P ∂t = ε 0 (ε s -ε ∞ ) τ E - 1 τ P .
(2.6)

The boundary conditions are usually of the following type: we introduce the decomposition ∂Ω = Γ m ∪ Γ a of the boundary of Ω on which we impose

   n × E = 0 on Γ m , n × E - µ ε n × H × n = n × E inc - µ ε n × H inc × n on Γ a ,
where n denotes the unit outward normal to ∂Ω and E inc , H inc is a given incident field. The first boundary condition is called metallic (referring to a perfectly conductive surface) while the second one is called absorbing and takes the form of the Silver-Müller condition which is a first order approximation of the exact absorbing boundary condition.

To simplify the presentation in the rest of the paper, we assume that Γ a = ∅. However, we note that the analysis carried on in the following can be easily extended to the more general setting that includes both a metallic and an absorbing boundary condition.

Remark 2.1. One can prove that if the fields are initially divergence free, then the fields are divergence free at any time, see [START_REF] Li | Unified analysis of leap-frog methods for solving time-domain Maxwell's equations in dispersive media[END_REF]. That is why we do not need to consider the third and fourth equations of (2.1) in the system (2.6).

The question of the existence and uniqueness of solutions of system (2.6) has been studied in [START_REF] Lanteri | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF]. Let us introduce H (curl, Ω) the classical subspace of L 2 (Ω) fields with curl in L 2 (Ω) and H 0 (curl, Ω) the classical subspace of H (curl, Ω) fields with zero tangential trace. Let us define

H 0 (Ω) = H (curl, Ω) 3 × H 0 (curl, Ω) 3 × L 2 (Ω) 3 ,
it is proven in [START_REF] Lanteri | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] that if we assume that H 0 , E 0 , P 0 belongs to H then system (2.6) has a unique weak solution belonging to C ([0, T ] , H 0 (Ω)). From now on, we consider a normalized form of Maxwell's equations (2.6). We introduce the normalized space, time variables and physical fields through the relations

x = x, t = c 0 t, τ = c 0 τ, σ = Z 0 σ, E = E, H = Z 0 H and P = P /ε 0 ,
where c 0 = 1/ √ ε 0 µ 0 is the speed of light in vacuum and Z 0 = µ 0 /ε 0 is the free space intrinsic impedance. For convenience of presentation we omit in the sequel the "∼" notation except for τ to avoid any ambiguity since τ will denote an element of a mesh, and we denote the relative magnetic permeability, µ r , by µ and use the notation ∆ε = ε s -ε ∞ . We can now write Maxwell's equations as

                 µ ∂ H ∂t = -curl E , ε ∞ ∂ E ∂t = curl H - ∆ε τ E -σ E + 1 τ P , ∂ P ∂t = ∆ε τ E - 1 τ P .
(2.7)

3. Locally implicit DGTD method for Maxwell's equations in dispersive media. In this section Maxwell's equations (2.7) are discretized in space according to a DG formulation then we adapt the locally implicit time integration method from [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin Time-domain method for the 3d heterogeneous Maxwell equations on unstructured meshes[END_REF] with the additional ADE and analyze its stability, via an energy approach which provides a rigorous stability criterion.

3.1. Semi-discretization in space by the DG method. As in [START_REF] Lanteri | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] we write the semi-discrete system of (2.7) using a DG method formulated on simplicial meshes following the approach of [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin Time-domain method for the 3d heterogeneous Maxwell equations on unstructured meshes[END_REF]. From now on we assume that Ω is a bounded convex polygonal domain, let Ω h be a partition of Ω into a set of N h tetrahedra τ i of size h i with boundary ∂τ i . We assume that this mesh is shape regular. By convention h denotes the maximum diameter of the (non-uniform) grid elements. An internal face is denoted by a ik = τ i ∩ τ k , for i = k and n ik is the unitary normal vector oriented from τ i to τ k . For the boundary interface, the index k corresponds to a fictitious element outside the domain. We denote by ν i the set of indices of the elements which have a common interface with τ i with the following decomposition

ν i = ν i i ∪ ν m i , with ν i i = {k ∈ ν i : a ik is an internal interface} and ν m i = {k ∈ ν i : a ik ∈ Γ m }.
As with a finite element method for the given partition Ω h we seek approximate electric and magnetic fields and electric polarization E h , H h , P h of E, H, P solution of (2.7) in a subspace V h . We introduce P pi (τ i ) the space of polynomial functions of degree at most p i inside the element τ i and define the following finite dimensional subspace

V h = v = (v 1 , v 2 , v 3 ) T ∈ L 2 (Ω) 3 , v k τi ∈ P pi (τ i ), ∀k, ∀τ i ∈ Ω h .
Following again the DG approach introduced in [START_REF] Lanteri | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF], inside each finite element τ i , the local electric, magnetic fields and electric polarization ( E h τi , H h τi , P h τi ) = ( E i , H i P i ) are expressed as linear combination of linearly independent basis vector Φ il , 1 ≤ l ≤ 3d i ,

E i ( x, t) = 3di l=1 E il (t)Φ il ( x), H i ( x, t) = 3di l=1 H il (t)Φ il ( x), P i ( x, t) = 3di l=1 P il (t)Φ il ( x), (3.1)
where d i denotes the local number of degrees of freedom associated to the interpolation degree p i in τ i , i.e.

d i = (p i + 1) (p i + 2) (p i + 3) /6
                         τi µ i Φ•∂ t H h d x + τi curlΦ• E h d x - ∂τi Φ•( E h × n) ds = 0, τi ε i∞ Φ•∂ t E h d x - τi curlΦ• H h d x + ∂τi Φ•( H h × n) ds + τi σ i Φ• E h d x - τi ∆ε i τ i Φ• E h d x + τi Φ τ i • P h d x = 0, τi Φ•∂ t P h d x - τi ∆ε i τ i Φ• E h d x + τi Φ τ i • P h d x = 0.
(3.2) One of the main features of the DG approach is that the approximate fields are allowed to be discontinuous across element boundaries. Then, for such discontinuous fields, we must define an approximate trace (i.e. a numerical trace) to evaluate the integrals over ∂τ i . In this study, we choose to use a centered approximation

∀i, ∀k ∈ ν i , E h a ik = E i a ik + E k a ik 2 and H h a ik = H i a ik + H k a ik 2 . (3.3)
Now using ( E h τi , H h τi , P h τi ) = ( E i , H i , P i ) for volume integrals, evaluating the surface integrals in (3.2) with the centered numerical flux (3.3) and re-integrating by parts, yields

                                           τi µ i Φ•∂ t H i d x = - 1 2 τi curlΦ• E i + curl E i •Φ d x + 1 2 k∈νi a ik Φ• E k × n ik ds, τi ε i∞ Φ•∂ t E i d x = 1 2 τi curlΦ• H i + curl H i •Φ d x - 1 2 k∈νi a ik Φ• H k × n ik ds - τi σ i Φ• E i d x - τi ∆ε i τ i Φ• E i d x + τi Φ τ i • P i d x, τi Φ•∂ t P i d x = τi ∆ε i τ i Φ• E i d x - τi Φ τ i • P i d x.
(3.4)

Then for each face on the boundary of Ω h , a ik , the trace of a fictitious neighboring element is needed for the computation of the numerical flux, we impose E k a ik = -E i a ik and H k a ik = H i a ik . The system of equations (3.4) can be written in terms of scalar unknowns. According to the decomposition (3.1) of the local electric, magnetic fields and electric polarization ( E i , H i , P i ), denoting the column vectors (E il ) 1≤l≤3di , (H il ) 1≤l≤3di and (P il ) 1≤l≤3di by E i , H i and P i , respectively, and replacing Φ by the basis vectors (Φ ij ) 1≤j≤3di , we obtain the equivalent system

                   M µ i ∂ t H i = -K i E i + k∈νi S ik E k , M ε∞ i ∂ t E i = K i H i - k∈νi S ik H k -D i E i -M ∆ε/ τ i E i + M 1/ τ i P i , M 1 i ∂ t P i = M ∆ε/ τ i E i -M 1/ τ i P i (3.5)
with the following matrices : M ς i (ς stands for µ, ε, ∆ε/ τ or 1/ τ ) are the symmetric, positive definite mass matrices

(M ς i ) jl = τi Φ T ij ς i Φ il d x (1 ≤ j, l ≤ 3d i ), K i is the symmetric stiffness matrix (K i ) jl = 1 2 τi Φ T ij curlΦ il + Φ T il curlΦ ij d x (1 ≤ j, l ≤ 3d i ),
S ik are the rectangular interface matrices

(S ik ) jl = 1 2 a ik Φ T ij (Φ kl × n ik ) ds (1 ≤ j ≤ 3d i , 1 ≤ l ≤ 3d k ).
D i is the symmetric, positive semi-definite conduction matrix

(D i ) jl = τi Φ T ij σ i Φ il d x (1 ≤ j, l ≤ 3d i ).
Taking into account the trace of fictitious neighboring elements and using the previous decomposition ν i = ν i i ∪ ν m i , we can rewrite the local system of ODEs (3.5) as

             M µ i ∂ t H i = -K i E i + k∈ν i i S ik E k - k∈ν m i S m ik E i , M ε∞ i ∂ t E i = K i H i - k∈ν i i S ik H k + k∈ν m i S m ik H i -D i E i -M ∆ε/ τ i E i + M 1/ τ i P i , M 1 i ∂ t P i = M ∆ε/ τ i E i -M 1/ τ i P i , (3.6)
with (S m ik

) jl = 1 2 a ik Φ T ij (Φ il × n ik ) ds (1 ≤ j, l ≤ 3d i ).
These set of local semidiscrete systems (3.6) for each finite element τ i can be transformed into a global system. Gathering all electric, magnetic and electric polarization unknowns in column vectors of size d = N h i=1 d i , denoted E, H and P , respectively, we write

       M µ ∂ t H = -KE + AE -BE, M ε∞ ∂ t E = KH -AH + BH -DE -M ∆ε/ τ E + M 1/ τ P, M 1 ∂ t P = M ∆ε/ τ E -M 1/ τ P,
where M ς and K are 3d × 3d block diagonal mass and stiffness matrices with diagonal blocks equal to M ς i and K i , respectively. Then the matrices M ς are symmetric, positive definite and K is symmetric. A is a 3d × 3d block sparse matrix, whose nonzero blocks are equal to S ik when a ik is an internal interface. From (3.1) and the equality n ki = -n ik we can check that S ki = S T ik and then A is symmetric

(S ki ) lj = 1 2 a ik -Φ T kl (Φ ij × n ik ) ds = 1 2 a ik Φ T ij (Φ kl × n ik ) ds = (S ik ) jl .
B is a 3d×3d block diagonal matrix, whose nonzero blocks equal to S m ik . In the present case (S m ik ) jl = -(S m ik ) lj , then S m ik = -(S m ik ) T and B is skew-symmetric. To simplify the notations and for the rest of this section we assume the different coefficients are constant and introduce the notation S = K -A -B, allowing now to drop the ∼ on τ we then obtain the following system

             µM ∂H ∂t = -S T E, ε ∞ M ∂E ∂t = SH - (ε s -ε ∞ ) τ M E -σM E + 1 τ M P, M ∂P ∂t = (ε s -ε ∞ ) τ M E - 1 τ M P. (3.7)
Now we use the Cholesky factorization of the mass matrix M = L M L T M , where L M is a triangular matrix, and we introduce the change of variables

E = L T M E, H = L T M H and P = L T M P in (3.7), hence                  µ ∂ H ∂t = -S T E, ε ∞ ∂ E ∂t = S H - (ε s -ε ∞ ) τ E -σ E + 1 τ P , ∂ P ∂t = (ε s -ε ∞ ) τ E - 1 τ P , (3.8)
where

S = L -1 M S (L -1 M ) T .
For convenience of presentation we again omit from now on the notation tilde in (3.8) and we focus now on the stability of this semi-discrete scheme. In the rest of the paper we denote • , • the discrete L 2 inner product on R d and • 2 the corresponding norm. Let (H h , E h , P h ) be the solution of the semi-discrete system (3.8), defining the semi-discrete energy by

E h = 1 2 µ H h 2 2 + ε ∞ E h 2 2 + 1 (ε s -ε ∞ ) P h 2 2 ,
as in [START_REF] Lanteri | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF], we can prove that E h is a decreasing function in time so that E h (t) ≤ E h (0), showing stability in the discrete L 2 sense on a fixed mesh.

3.2. Time integration methods. Adapting the locally implicit time integration method from [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] with an additional ADE we introduce the implicit-explicit time integration scheme for the semi-discrete system (3.8)

                                             µ H n+1/2 -H n ∆t/2 = -S T E n , ε ∞ E n+1 -E n ∆t = S 0 H n+1/2 + 1 2 S 1 H n+1 + H n - (ε s -ε ∞ ) 2τ E n+1 + E n - 1 2 σ E n+1 + E n + 1 2τ P n+1 + P n , P n+1 -P n ∆t = (ε s -ε ∞ ) 2τ E n+1 + E n - 1 2τ P n+1 + P n , µ H n+1 -H n+1/2 ∆t/2 = -S T E n+1 , (3.9)
where S = S 0 + S 1 is a matrix splitting. The method is implicit in S 1 and explicit in S 0 . For S 0 = 0 we recover the second order leap-frog scheme and for S 1 = 0 the second order Crank-Nicolson scheme both written in a three-stage form. We adopt the splitting defined in [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF], i.e. S 1 = SS H , where S H is a diagonal matrix of dimension the length of H defined by (S H ) jj = 0, component H j of H to be treated explicitly, 1, component H j of H to be treated implicitly.

From the third and the fourth equations of (3.9) we express the electric polarization P n+1 and the magnetic field H n+1 as

P n+1 = 2τ -∆t 2τ + ∆t P n + ∆t 2τ + ∆t (ε s -ε ∞ ) E n+1 + E n ,
and

H n+1 = H n+1/2 - 1 µ ∆t 2 S T E n+1 .
Now we substitute these expressions into the second stage of (3.9)

ε ∞ E n+1 -E n ∆t = S 0 H n+1/2 + 1 2 S 1 (H n + H n+1 ) - 1 2τ (ε s -ε ∞ ) E n - 1 2 σE n + 1 2τ ∆t 2τ + ∆t (ε s -ε ∞ ) E n - 1 2τ (ε s -ε ∞ ) E n+1 - 1 2 σE n+1 + 1 2τ ∆t 2τ + ∆t (ε s -ε ∞ ) E n+1 - 1 µ ∆t 4 S 1 S T E n+1 + 1 2τ P n + 1 2τ 2τ -∆t 2τ + ∆t P n .
(3.10) Note that S 1 S T = S 1 S T 1 , the second stage (3.10) is then equivalent to the following linear system

ME n+1 = b n , (3.11)
where

M = 1 - 1 ε ∞ ∆t 2τ ∆t 2τ + ∆t (ε s -ε ∞ ) + 1 ε ∞ ∆t 2τ (ε s -ε ∞ ) + 1 ε ∞ ∆t 2 σ I + 1 ε ∞ µ ∆t 2 4 S 1 S T 1 , (3.12
) and

b n = 1 + 1 ε ∞ ∆t 2τ ∆t 2τ + ∆t (ε s -ε ∞ ) - 1 ε ∞ ∆t 2τ (ε s -ε ∞ ) - 1 ε ∞ ∆t 2 σ E n + 1 ε ∞ ∆tS 0 H n+1/2 + 1 ε ∞ ∆t 2 S 1 H n+1/2 + H n + 1 ε ∞ ∆t 2τ 1 + 2τ -∆t 2τ + ∆t P n .
(3.13) Hence we can write the locally implicit scheme (3.9) as

                       H n+1/2 = H n - 1 µ ∆t 2 S T E n , ME n+1 = b n , P n+1 = 2τ -∆t 2τ + ∆t P n + ∆t 2τ + ∆t (ε s -ε ∞ ) E n+1 + E n , H n+1 = H n+1/2 - 1 µ ∆t 2 S T E n+1 , (3.14)
where the matrix M and the right hand side b n are defined by (3.12) and (3.13), respectively.

3.3. Stability of the fully discrete scheme on a fixed mesh. In this section we are interested in the stability of the fully discrete locally implicit scheme (3.9). The derivations in the remainder of this section follow an energy approach which provides a rigorous criterion for stability. We first exhibit a discrete energy, which is a quadratic form of the numerical unknowns, we then show that this energy is a positive definite quadratic form and finally we prove that this energy is decreasing, which achieves the stability analysis. We define the discrete electromagnetic energy, denoted E n , as

E n = 1 2 µ H n 2 2 + ε ∞ E n 2 2 + 1 (ε s -ε ∞ ) P n 2 2 - ∆t 2 4µ S 0 S T E n , E n , (3.15)
We first state a condition on the time step ∆t such that E n is a positive definite quadratic form. Lemma 3.1. The quadratic form E n given by (3.15) is a positive definite quadratic form of the numerical unknowns H n , E n and P n if

∆t < 2 √ ε ∞ µ ρ S 0 S T 0 , (3.16)
where ρ(S 0 S T 0 ) denotes the spectral radius of S 0 S T 0 . Proof. By definition of S 0 = S(I -S H ) the matrix S 0 S T is symmetric, since S 0 S T = S(I -S H )S T = S(I -S H )(I -S H )S T = S 0 S T 0 . Then, from (3.15) we write the electromagnetic energy as

E n = 1 2 µ H n 2 2 + ε ∞ E n 2 2 + 1 (ε s -ε ∞ ) P n 2 2 - ∆t 2 4µ S 0 S T 0 E n , E n , = 1 2 µ H n 2 2 + ε ∞ E n 2 2 + 1 (ε s -ε ∞ ) P n 2 2 - ∆t 2 4µ S 0 E n 2 2 .
Furthermore, we have

S T 0 E n 2 ≤ S T 0 2 E n 2 = ρ S 0 S T 0 E n 2 , hence E n ≥ 1 2 µ H n 2 2 + ε ∞ - ∆t 2 4µ ρ S 0 S T 0 E n 2 2 + 1 (ε s -ε ∞ ) P n 2 2 ,
from which it follows that under the condition (3.16), the discrete energy E n is a positive definite quadratic form of the numerical unknowns H n , E n and P n . We now prove the following result. Lemma 3.2. The discrete energy (3.15) is decreasing so that E n ≤ E 0 . Proof. From the first and fourth equation of (3.9) we have

H n+1/2 = H n - ∆t 2µ S T E n , H n+1/2 = H n+1 + ∆t 2µ S T E n+1 .
(3.17)

By substituting half of each expression of (3.17) for H n+1/2 into the second stage of (3.9) and the first expression of (3.17) into the fourth stage of (3.9), together with the final expression of P n+1 we obtain

ε ∞ E n+1 -E n = ∆t 2 S H n+1 + H n + ∆t 2τ P n+1 + P n - (ε s -ε ∞ ) ∆t 2τ E n+1 + E n - ∆t 2 σ E n+1 + E n + ∆t 2 4µ S 0 S T E n+1 -E n , 1 (ε s -ε ∞ ) P n+1 -P n = ∆t 2τ E n+1 + E n - ∆t 2τ (ε s -ε ∞ ) P n+1 + P n , µ H n+1 -H n = - ∆t 2 S T E n+1 + E n .
Taking inner product with E n+1 +E n , P n+1 +P n and H n+1 +H n yields, respectively

ε ∞ E n+1 2 2 -E n 2 2 = ∆t 2 S H n+1 + H n , E n+1 + E n + ∆t 2τ P n+1 + P n , E n+1 + E n - (ε s -ε ∞ ) ∆t 2τ E n+1 + E n 2 2 - ∆t 2 σ E n+1 + E n 2 2 + ∆t 2 4µ S 0 S T E n+1 -E n , E n+1 + E n , 1 (ε s -ε ∞ ) P n+1 2 2 -P n 2 2 = ∆t 2τ E n+1 + E n , P n+1 + P n - ∆t 2τ (ε s -ε ∞ ) P n+1 + P n 2 2 , µ H n+1 2 2 -H n 2 2 = - ∆t 2 S T E n+1 + E n , H n+1 + H n .
(3.18) We recall that S 0 S T = S 0 S T 0 , hence

S 0 S T E n+1 -E n , E n+1 + E n = S 0 S T 0 E n+1 -E n , E n+1 + E n , = S 0 S T 0 E n+1 , E n+1 -S 0 S T 0 E n , E n , = S 0 S T E n+1 , E n+1 -S 0 S T E n , E n .
Substituting this expression into the first equation of (3.18) and adding the three equations of (3.18) yields

E n+1 -E n ∆t = - 1 4τ (ε s -ε ∞ ) P n+1 -P n 2 2 + 1 2τ P n+1 + P n , E n+1 + E n - (ε s -ε ∞ ) 4τ E n+1 -E n 2 2 - 1 4 σ E n+1 -E n 2 2 , then E n+1 -E n ∆t = - 1 4τ (ε s -ε ∞ ) P n+1 + P n -(ε s -ε ∞ ) E n+1 + E n 2 2 - 1 4 σ E n+1 -E n 2 2 , ≤ 0.
This concludes the proof of Lemma 3.2.

Convergence.

In this section we are interested in the PDE convergence of the locally implicit method (3.9). As in [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] we wonder whether the method retains its second-order ODE convergence under stable simultaneous space-time grid refinement ∆t h, h → 0 towards the exact PDE solution. This is not a priori clear due to the component splitting which can introduce order reduction through error constants which grow with h -1 , for h → 0, see for example [START_REF] Moya | Temporal convergence of a locally implicit Discontinuous Galerkin method for Maxwell's equations[END_REF]. By applying the method described in [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] and [START_REF] Moya | Temporal convergence of a locally implicit Discontinuous Galerkin method for Maxwell's equations[END_REF], it is possible to prove that the subdivision into coarse and fine elements is not detrimental to the second-order ODE convergence of the method (3.9), under stable simultaneous space-time grid refinement towards the exact underlying PDE solution. This convergence result is given in the following theorem :

Theorem 3.3. Let H h (t), E h (t) and P h (t) denote the exact solutions of the Maxwell problem in dispersive media under consideration, restricted to the space grid, i.e. the exact solutions of the system of ODEs

                     µ d dt H h (t) = -S T E h (t) + ζ H h (t) , ε ∞ d dt E h (t) = SH h (t) - (ε s -ε ∞ ) τ E h (t) -σE h (t) + 1 τ P h (t) + ζ E h (t) , d dt P h (t) = (ε s -ε ∞ ) τ E h (t) - 1 τ P h (t) + ζ P h (t) , (3.19)
where ζ H h , ζ E h and ζ P h denote the spatial truncation errors. Assume a Lax-Richtmyer stable space-time grid refinement ∆t ∼ h, h → 0. On the interval [0, T ] the approximations H n , E n and P n of method (3.9) then converge with order two to H h (t), E h (t) and P h (t).

Remark 3.1. The proof of this theorem can be found in [START_REF]Locally implicit Discontinuous Galerkin Time-domain methods for electromagnetic wave propagation in biological tissues[END_REF] and more recently in [START_REF] Descombes | Temporal convergence analysis of a locally implicit discontinuous galerkin time domain method for electromagnetic wave propagation in dispersive media[END_REF]. Futhermore, following the previous results and the techniques developed in the recent prepint [START_REF] Hochbruck | Error analysis of a second order locally implicit method for linear Maxwell's equations[END_REF] in which an error analysis of the fully discrete locally implicit DGDT scheme is presented in the case of non-dispersive media, it is also possible to obtain an error analysis in our dispersive case.

4. Numerical results. The simulations discussed in this section have been performed on a workstation equipped with an Intel Xeon 2.40 GHz processor and 16 GB of RAM memory. The linear system associated to the locally implicit method (3.9) (see (3.11)) is solved using the MUMPS (MUltifrontal Massively Parallel sparse direct Solver) optimized sparse direct solver [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF]. We consider the three-dimensional (normalized) Maxwell's equations in Debye dispersive media (2.7) where we add a given source current density. We denote by Ω h the computational domain and Ω exp h the set of tetrahedra that belong to the region where the explicit method is used into the implicit-explicit splitting. The critical step size, denoted ∆t c , used in the numerical tests is then given by

∆t c = CF L × min τi∈Ω exp h 4 × V i A i , (4.1)
where V i and A i are the volume and the surface area of tetrahedron τ i , respectively. The values of the CF L number corresponds to the numerical stability, i.e. the limit beyond which we observe a growth of the discrete energy.

4.

1. An artificial validation test. We establish an artificial three-dimensional problem to validate the locally implicit DGTD method proposed here for solving the time domain Maxwell's equations in dispersive media with the Debye model. For that purpose, the computational domain is the cube [0, 1] 3 and we first substitute the solution of the propagation of a standing wave in a cubic PEC cavity given by, for a real ω Then we substitute the electric field and the polarization into the second equation of (2.7) to obtain the following current density to impose to have an exact solution

H = π ω sin (ωt)
                                 J s,x ( x, t) = cos (πx) sin (πy) sin (πz) × ε s -ε ∞ 1 + ω 2 τ 2 1 τ cos (ωt) + ω sin (ωt) - ε s -ε ∞ τ + σ cos (ωt) , J s,y ( x, t) = 0, J s,z ( x, t) = ε s -ε ∞ τ sin (πx) sin (πy) sin (πz) cos (ωt) × ε s -ε ∞ 1 + ω 2 τ 2 1 τ cos (ωt) + ω sin (ωt) - ε s -ε ∞ τ + σ cos (ωt) , (4.3
) the angular frequency is given by ω = 2πf (rad•s -1 ) with f = 260 MHz and wavelength λ = 1.15 m. We impose a metallic boundary condition n × E = 0 on ∂Ω. In our simulations we choose µ = ε ∞ = 1, ε s = 5, σ = 0 and τ = 9.4 picoseconds (ps). The total simulation time is set to T = 1.67 × 10 -8 s which corresponds to T = 5 m (normalized unit). We investigate the space-time convergence order (i.e. for a stable simultaneous space-time grid refinement ∆t h, h → 0) of the fully explicit method and the locally implicit methods. We measure the maximal L 2 -norm of the error for a sequence of four successively locally refined tetrahedral meshes, see Figure 4.1 for examples of meshes with the identification of the region for implicit treatment. Propagation of a standing wave in a cubic PEC cavity: cross sections of 3D meshes of two locally refined meshes used for the numerical convergence study (implicit treatment: red regions), on left: 2, 968 tetraedra, 536 in the implicit region, so approximately 18% of the elements of the mesh and 635 vertices; on right: 40, 616 tetraedra, 536 in the implicit region, so approximately 1, 3% of the elements of mesh, 7, 759 vertices.

We plot this error as a function of 1/h, in logarithmic scale, the use of the logarithmic scale allows to visualize the convergence rates as the slopes of the curves. We use the DGTD-P k method, with k = 1 or 2. The obtained results, for both methods, are summarized in Figure 4.2 (left). The results are similar for both time integration methods, which confirm that the subdivision into coarse and fine elements is not detrimental for the convergence order of the locally implicit method. The order of convergence is slightly stronger than the theoretical behavior for k = 1 (about 1.6 instead of 1.0), while for k = 2 we obtain the expected order, i.e. the order two. Indeed our DG approach is based on a centered numerical flux for the approximation of the boundary integral term at the interface between neighboring elements, then the convergence rate is O(h k ) for a k-th interpolation order [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin Time-domain method for the 3d heterogeneous Maxwell equations on unstructured meshes[END_REF]. The convergence result is slightly weaker than available results for upwind fluxes with O(h k+1 ) convergence rate [START_REF] Cockburn | Discontinuous Galerkin Methods[END_REF][START_REF] Sármány | Dispersion and dissipation error in highorder runge-kutta Discontinuous Galerkin discretisations of the Maxwell equations[END_REF]. Consequently a suboptimal convergence rate is obtained in the L 2 -norm; nevertheless, this setting yields a non-dissipative numerical scheme which preserves a discrete form of the electromagnetic energy in the case on non-dispersive and nonconducting media, assuming metallic boundary conditions only link in [START_REF] Dolean | Locally implicit Discontinuous Galerkin method for time domain electromagnetics[END_REF] and [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin Time-domain method for the 3d heterogeneous Maxwell equations on unstructured meshes[END_REF]. We also plot in Figure 4.2 (right) the error as a function of the CPU time. For a given error we can observe the gain in CPU time with the locally implicit approach compared to the fully explicit method. Now we focus on the linear systems to be solved and the efficiency of the locally implicit method. We consider the locally refined tetrahedral mesh composed of 40, 616 tetrahedra and 7, 756 vertices (Figure 4.1, on right). First, we are interested in the sparsity of the matrix to be inverted and the cost of the factorization step. In Table 4.1 for each matrix of the linear system to solve, we indicate the number of nonzero elements and the fill-in ratio percentage. We also indicate the total size (memory requirement) of all internal data used during nu- merical factorization and the total CPU time for analysis and factorization. Note that for the fully implicit scheme based on DGTD-P 1,2 methods the number of nonzero is 75, 081, 406 and 408, 726, 916 respectively, and will require too much memory to be a reasonable alternative to the fully explicit method. Finally we present some numerical results for the locally implicit method and the fully explicit method. In Table 4.2 we observe that the locally implicit method allows to overcome the step size limitations caused by the local refinement. With implicit-explicit approaches the sizes of the time step are about 34 times larger which yields significant gains of final CPU time, about 8 times lower than the fully explicit case. Regarding the time evolution of the error, Figure 4.3, we observe that the IMEX scheme and the fully explicit scheme based on DGTD-P 4.2. Exposure of head tissues to a localized source radiation. We now consider a more realistic test problem in the field of bioelectromagnetics. This problem is of particular importance to the assessment of potential adverse effects resulting from human exposure to electromagnetic fields radiated from wireless systems. More precisely, the goal is to simulate the propagation of an electromagnetic wave in the head tissues when the wave is emitted from a mobile phone antenna. Here, we do not take into account in the numerical modeling a detailed phone structure and simply consider a dipolar source localized near an ear, and concentrate on demonstrating the usefulness of the proposed locally implicit DGTD method (3.9). The geometrical model of the head is heterogeneous and consists of four tissues namely, the skin, the skull, the Cerebro Spinal Fluid (CFS) and the brain. The computational domain is artificially bounded by a sphere on which a Silver-Müller condition is imposed and we use an unstructured tetrahedral mesh consisting of 61, 358 vertices and 366, 208 tetrahedra. The surface meshes of the different tissues are shown on Figure 4.4 and the unstructured mesh on Figure 4.5. The Debye model parameters that we have used for the tissues are given in Table 4.3. The numbers have been obtained by applying an appropriate fitting procedure from a set of values for the permittivity of tissues at given frequencies, that have been obtained with a base available at: http://niremf.ifac.cnr.it/tissprop/. In the present case, there is no incident field and the source term is given by

Tissue

J s,z ( x, t) = f 0 e -(t-4tp) 2 /tp sin (2πf c (t -t p )) g ( x) g ( x) , (4.4)
where f 0 is the amplitude of the signal, t p = 400 ps and f c = 1.8 GHz, and g is a three-dimensional gaussian function with (x 0 , y 0 , z 0 ) = (-0.100, 0.025, -0.015) the center of the gaussian spatial support

g( x) = e -α((x-x0) 2 +(y-y0) 2 +(z-z0) 2 ) . (4.5)
The total physical simulation time is set to T = 10 × T c where T c = 1/f c , and a discrete Fourier transform of the electromagnetic field at the frequency f c is computed on the fly. For the simulations we only consider a DGTD-P 1 method with the implicit region consisting of 5092 tetrahedra, which represents 1.4% of the total number of elements. The parameter α in (4.5) has been chosen such that the source term J s,z is strongly localized, which justifies the local refinement around the center of the gaussian, the support of the gaussian g is then strictly included in the implicit region.

The localization and the distribution of the elements for the implicit treatment are given in Figure 4.5 and Table 4.4, respectively. Finally, Table 4.5 gives the number of nonzero elements, the fill-in ratio percentage for the matrix of the linear system to be solved, the total size (memory requirement) of all internal data used during numerical factorization and the total CPU time for analysis and factorization. 4.5 Exposure of head tissues to a localized source radiation: data and factorization of the matrix of the linear system to be solved for the locally implicit DGTD-P 1 method.

A first simulation is performed with the locally implicit DGTD-P 1 method (3.9) and a second simulation is realized using the original locally implicit DGTD-P 1 method for non-dispersive media with the same localized source term and adopting the appropriate values of the electric permittivity and the electric conductivity for each tissue at the central frequency f c . Time evolution of the electric component E z at two selected points is shown on Figure 4.6 for the dispersive and the nondispersive calculations. As expected, the curves obtained for the first point, localized at the center of the gaussian in vacuum are similar. For the second point, localized in the head, we can observe the dispersion and a slight attenuation of the signal for the dispersive calculation compared to the non-dispersive calculation.

A third simulation is performed with the fully explicit DGTD-P 1 method. We recall that the critical time step size for the stability of the locally implicit method is then about 15.5 times larger than the critical time step size for the stability of the fully explicit method. We recall that the SAR is a measure of the rate at which electric energy is absorbed by the tissues when exposed to a radio-frequency electromagnetic field. It represents the power absorbed per mass of tissue and has units of watts per kilogram (W•kg -1 ). The SAR is then defined as σ| E f our | 2 /ρ, where E f our denotes the electric field in the frequency-domain, resulting from the discrete Fourier transform of the temporal field, and ρ is the density which depends on the tissues. We do not observe noticeable differences in the different tissues for the SAR patterns, the obtained results for the locally implicit method (3.9) and the original locally implicit method are globally similar.

On Figures 4.7 we show the contour lines of the local SAR normalized by the maximal local SAR, in logarithmic scale, for the calculations with the locally implicit DGTD-P 1 method (3.9). For convenience of presentation we do not show the results obtained with the fully explicit method because they are indistinguishable from that shown in Figure 4.7. Finally, the total simulation times are equal to 6 h 38 min for the locally implicit DGTD-P 1 method and 24 h 56 min for the fully explicit DGTD-P 1 method. Hence, the locally implicit DGTD-P 1 method allows a reduction of the computing time by a factor of 3.8.

Conclusion and future works.

We have presented a locally implicit DGTD method for simulating the propagation of electromagnetic waves in Debye-type dispersive media in view of biomedical applications. The latter most often involve irregularly shaped structures corresponding to biological tissues. Modeling realistically the interfaces between tissues is particularly important if one is interested in evaluating accurately the impact of field discontinuities on the local SAR distribution. DGTD methods formulated on unstructured tetrahedral meshes are appealing in this context but can be very expensive in terms of computing time when based on a fully explicit time integration scheme. The locally implicit solution strategy considered in this study allows to overcome the restrictive stability condition of global time step explicit schemes at the expense of the inversion of large sparse linear system at each time step. Despite this computational overhead, numerical and performance results obtained here for three-dimensional problems clearly illustrate the benefits of the proposed locally implicit DGTD method. Future works will be concerned with the algorithmic adaptation of the method to parallel computing architectures and the study of high order in time approximations; indeed, composition methods allow to reach higher order as seen in [START_REF] Descombes | High-order locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell's equations[END_REF] for a non-dispersive case. Noteworthy, the method can be easily extended to other dispersion models such as the Drude and Drude-Lorentz models for the interaction of electromagnetic waves with metals at optical frequencies.

-- 1 + ω 2 τ 2 1 τ

 11 sin (πx) cos (πy) cos (πz) 2 cos (πx) sin (πy) cos (πz) -cos (πx) cos (πy) sin (πz) cos (πx) sin (πy) sin (πz) 0 sin (πx) sin (πy) cos (πz)   , into the third equation of (2.7) to obtain the analytic expression of the electric polarization and find P = -(ε s -ε ∞ ) τ cos (ωt) + ω sin (ωt)   cos (πx) sin (πy) sin (πz) 0 sin (πx) sin (πy) cos (πz)

Fig. 4 .

 4 Fig.4.1. Propagation of a standing wave in a cubic PEC cavity: cross sections of 3D meshes of two locally refined meshes used for the numerical convergence study (implicit treatment: red regions), on left: 2, 968 tetraedra, 536 in the implicit region, so approximately 18% of the elements of the mesh and 635 vertices; on right: 40, 616 tetraedra, 536 in the implicit region, so approximately 1, 3% of the elements of mesh, 7, 759 vertices.

Fig. 4 . 2 .

 42 Fig. 4.2. Numerical convergence and maximum error (L 2 -norm) in function of final CPU time for the locally implicit and fully explicit DGTD-P k methods (left -right, respectively).

Fig. 4 . 3 .

 43 Fig. 4.3. Time evolution of the error (L 2 -norm) for the locally implicit and fully explicit DGTD-P k methods.

Fig. 4 . 4 .

 44 Fig. 4.4. Surface meshes of the skin, the skull and the CSF.

Fig. 4 . 5 .

 45 Fig. 4.5. Cross sections of the 3D mesh (61, 358 vertices, 366, 208 tetrahedra, implicit treatment: 5092 tetrahedra, dark red region).

Fig. 4 . 6 .

 46 Fig. 4.6. Exposure of head tissues to a localized source radiation: time evolution of the E z component of the electric field at selected spatial locations (center of the gaussian and in the head).

Fig. 4 . 7 .

 47 Fig. 4.7. Exposure of head tissues to a localized source radiation: contour lines of the local SAR normalized by the maximal local SAR (logarithmic scale).

Table 4 . 1

 41 1,2 methods give similar results. Data and factorization of the matrix of the linear system to be solved for the locally implicit methods, with DGTD-P k methods (k = 1, 2).

	P k	Matrix		# nonzero	Fill-in ratio	RAM size CPU time
		order		(nz)		(percentage) (MB)	(s)
	P 1	487,392	1,247,266	0.000525	210	8
	P 2	1,218,480	5,437,770	0.000366	937	72
		Fully exp. method		Loc. imp. method
	P k	∆t (m)	CPU time (s)	∆t (m)	CPU time (s)	CPU (LF2) CPU (IMEX)
	P 1	1.1247e-4	7285		2.6040e-3	922	7.9
	P 2	7.8267e-5	31568		4.1502e-3	3786	8.3

Table 4 . 2

 42 Critical time step size and CPU time for the fully explicit method and the locally implicit method, with DGTD-P k method (k = 1, 2).

Table 4 . 3

 43 Debye model parameters for the different tissues.

Table 4 . 4

 44 Exposure of head tissues to a localized source radiation: distribution of elements in the implicit region.

	Vacuum	Skin	Skull	CSF	Brain	Total
	4209		102	21	720	40	5092
	P k	Matrix order	# nonzero Fill-in ratio	RAM size	CPU time
				(nz)	(percentage) (MB)	(s)
	P 1	4394496	15048090	0.000078	1298	95

Table


Emmanuel Vincent
email: emmanuel.vincent@inria.fr

Aghilas Sini

François Charpillet

AUDIO SOURCE LOCALIZATION BY OPTIMAL CONTROL OF A MOBILE ROBOT

Keywords: Source localization, occupancy grid, active sensing, mobile robot control

We consider the task of audio source localization using a microphone array on a mobile robot. Active localization algorithms have been proposed in the literature that can estimate the 3D position of a source by fusing the measurements taken for different poses of the robot. The robot movements are typically fixed, however, or they obey heuristic strategies, such as turning the head and moving towards the source, which may be suboptimal. In this paper, we propose to control the robot movements so as to locate the source as quickly as possible. We represent the belief about the source position by a discrete grid and we introduce a dynamic programming algorithm to find the optimal robot motion minimizing the entropy of the grid. We report initial results in a real environment.

INTRODUCTION

Robot audition is an emerging research field at the interface of audio signal processing, artificial intelligence, and control theory [START_REF] Nakadai | Design and implementation of robot audition system 'HARK' -open source software for listening to three simultaneous speakers[END_REF]. Today, assistive robots typically carry several microphones enabling them to locate and to recognize speech and other sound events. This enables the detection of visually hidden sound sources and efficient interaction with humans.

Source localization techniques fall into three classes [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF]. One approach is to compute the time delay of arrival (TDOA) between every two microphones using generalized crosscorrelation with phase transform (GCC-PHAT) [START_REF] Knapp | The generalized crosscorrelation method for estimation of time delay[END_REF] and to derive the source position by triangulation. It is typically outperformed by steered response power (SRP) [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF] or multiple signal classification (MUSIC) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] techniques that compute the pseudo-likelihood of each candidate position on a grid and pick the maxima on that grid. See [START_REF] Dibiase | Robust localisation in reverberant rooms[END_REF][START_REF] Badali | Evaluating real-time audio localization algorithms for artificial audition in robotics[END_REF][START_REF] Blandin | Multi-source TDOA estimation in reverberant audio using angular spectra and clustering[END_REF] for experimental comparisons. Binaural extensions of the above techniques have been designed for situations when the array is mounted on a head [START_REF] Nakadai | Active audition for humanoid[END_REF]. Particle filtering-based tracking algorithms have also been studied for moving sources [START_REF] Lehmann | Particle filter design using importance sampling for acoustic source localization and tracking in reverberant environment[END_REF][START_REF] Valin | Robust localization and tracking of simultaneous moving sound sources using beamforming and particle filtering[END_REF][START_REF] Fallon | Acoustic source localization and tracking of a time-varying number of speakers[END_REF].

All three categories of techniques have been implemented on robots [START_REF] Nakadai | Active audition for humanoid[END_REF][START_REF] Martinson | Auditory evidence grids[END_REF][START_REF] Kim | Speaker localization using the TDOA-based feature matrix for a humanoid robot[END_REF][START_REF] Nakamura | Real-time super-resolution sound source localization for robots[END_REF]. They are most often used to estimate the source angle of arrival (AoA), as with a static far-field microphone array. Mobile robots are not restricted to this, however, and they can estimate the full 3D position of the source by taking measurements at different poses. So-called active source localization algorithms [START_REF] Nakadai | Active audition for humanoid[END_REF] have been proposed to integrate these measurements into a single location estimate by means of triangulation [START_REF] Sasaki | Multiple sound source mapping for a mobile robot by self-motion triangulation[END_REF][START_REF] Kneip | Binaural model for artificial spatial sound localization based on interaural time delays and movements of the interaural axis[END_REF], occupancy grids [START_REF] Martinson | Auditory evidence grids[END_REF][START_REF] Dejong | Auditory occupancy grids with a mobile robot[END_REF], nonlinear extensions of Kalman filtering [START_REF] Portello | Acoustic models and Kalman filtering strategies for active binaural sound localization[END_REF], or particle filtering [START_REF] Marković | Active speaker localization with circular likelihoods and bootstrap filtering[END_REF]. Occupancy grids are one of the most successful frameworks for environment modeling in robotics [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF]. The algorithms in [START_REF] Martinson | Auditory evidence grids[END_REF][START_REF] Dejong | Auditory occupancy grids with a mobile robot[END_REF] estimate the posterior probability of each candidate position on a grid by rescaling the GCC-PHAT pseudolikelihood. Such an inverse sensor model treats each cell of the grid separately and it results in probability values which do not integrate to 1, which leads to inconsistent maps [START_REF] Thrun | Learning occupancy grid maps with forward sensor models[END_REF].

The above algorithms operate by making the robot follow a fixed patrolling path [START_REF] Martinson | Discovery of sound sources by an autonomous mobile robot[END_REF] or heuristic motion strategies inspired from humans [START_REF] Lu | Motion strategies for binaural localisation of speech sources in azimuth and distance by artificial listeners[END_REF], such as turning the head towards the source [START_REF] Nakadai | Active audition for humanoid[END_REF][START_REF] Berglund | Sound source localisation through active audition[END_REF] and getting closer to it [START_REF] Martinson | Discovery of sound sources by an autonomous mobile robot[END_REF][START_REF] Song | Olfaction and hearing based mobile robot navigation for odor/sound source search[END_REF]. Arm movements resulting in dynamic change of the array aperture have also been studied [START_REF] Barfuss | An adaptive microphone array topology for target signal extraction with humanoid robots[END_REF]. These strategies have been shown to improve localization performance experimentally but they may be suboptimal, especially for arrays of three or more microphones which do not fit the geometry of binaural hearing.

In this paper, we provide two contributions. Firstly, we represent the belief about the source position by a discrete grid and we update it over time using a rigorous forward sensor model [START_REF] Thrun | Learning occupancy grid maps with forward sensor models[END_REF]. Secondly, we quantify the information carried by the grid by its entropy and we introduce a dynamic programming algorithm to find the optimal robot motion minimizing the expected future entropy. The structure of the rest of the paper is as follows. The position mapping algorithm and the motion control algorithm are described in Sections 2 and 3, respectively. We report initial results in a real environment in Section 4 and we conclude in Section 5.

SOURCE POSITION MAPPING

We represent the environment as a discrete 3D grid whose cells are binary random variables encoding the presence or absence of an audio source at each position. The goal of mapping is to estimate the posterior probability of these variables given a set of measurements. The number of parameters of the posterior grows exponentially with the size of the grid, hence approximations are required.

A popular approximation is to assume independence and to estimate the occupancy of each cell separately [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF]. This approach used in [START_REF] Martinson | Auditory evidence grids[END_REF][START_REF] Dejong | Auditory occupancy grids with a mobile robot[END_REF] suffers from two drawbacks. First, the resulting maps are inconsistent: when the robot does not move, the probabilities will converge to 1 for all cells along the line from the microphone array to the source, while the source is present in one of those cells only. Secondly, the probabilities do not integrate to 1, which makes it difficult to quantify the amount of information carried by the grid.

Discrete grid

To address these issues, we adopt a rigorous approach based on a probabilistic forward sensor model [START_REF] Thrun | Learning occupancy grid maps with forward sensor models[END_REF]. We assume that there is a single active sound source located at absolute position s. We partition the time axis into time frames indexed by t and we denote by p t = [x t , y t , θ t ] the pose of the microphone array at time t, i.e., its absolute position [x t , y t , z] and its orientation θ t w.r.t. the y-axis. We further denote by m t the acoustic measurement at time t. The posterior probability of the source position to be estimated is then equal to

P (s = i|m 1:T , p 1:T ) (1)
where m 1:T = {m 1 , . . . , m T }, p 1:T = {p 1 , . . . , p T }, and i = [x i , y i , z i ] is any cell of the grid.

Forward sensor model

The sensor model p(m t |s = i, p t ) defines the likelihood that a given localization technique applied to one signal frame provides a measurement m t given the source position i and the robot pose p t . In the case of a linear far-field microphone array, the measurements m t are AoA estimates relative to the array and the sensor model can be expressed as

p(m t |s = i, p t ) = p(m t |d it , α it ) (2)
where

d it = [(x i -x t ) 2 + (y i -y t ) 2 + (z i -z) 2 ] 1/2
(3)

α it = arccos (y i -y t ) cos θ t -(x i -x t ) sin θ t d it (4)
are the distance and the AoA of cell i relative to the array at time t, respectively. More general parameterizations can be found for nonlinear or nonplanar arrays.

As an example, let us consider the scenario of a Turtlebot1 equipped with a Kinect. We estimate the source AoA from the Kinect's 4-microphone linear array output using MUSIC with generalized singular value decomposition (GSVD) [START_REF] Nakamura | Real-time super-resolution sound source localization for robots[END_REF] as implemented in HARK [START_REF] Nakadai | Design and implementation of robot audition system 'HARK' -open source software for listening to three simultaneous speakers[END_REF]. MUSIC-GSVD is a variant of MUSIC that is robust to spatially correlated noise and that operates in real time. The input covariance matrix R xx (t, f ) and the noise covariance matrix R nn (f ) are classically estimated in each time frame t and each frequency bin f by averaging the short time Fourier transform (STFT) of the input signal and a noise-only signal. The SVD of R nn (f ) -1 R xx (t, f ) is then computed and the left singular vectors are used to compute the MUSIC angular spectrum, whose maximum yields the estimated source AoA m t .

In order to learn the sensor model ( 2), we simulated speech recording in the presence of spatially isotropic Gaussian noise via the image method [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF] using Roomsimove2 . The room size, the reverberation time (250 ms), the intensity of speech and noise, and the noise spectrum roughly match those of the real environment in Section 4. The microphones were supposed to be omnidirectional since we did not have access to the head-related transfer functions (HRTFs) of the Kinect coupled with the Turtlebot. For each of 360 true AoAs (from 0 • to 359 • ) and 5 distances (from 0.5 to 3 m), we built the histogram of AoAs estimated by MUSIC-GSVD over 50 time frames and 100 random robot poses.

The resulting probability density is illustrated in Fig. 1 for two AoAs. At small distance, it concentrates around the true AoA and its symmetric w.r.t. the microphone axis, a phenomenon known as front-back confusion for humans. As distance increases, it becomes smeared and spurious peaks appear at 0 • and 180 • due to lower signal-to-noise ratio.

Fig. 1 also shows that the likelihood of correct measurement is higher for broadside AoAs (close to 90 • ) in the auditory fovea [START_REF] Nakadai | Active audition for humanoid[END_REF][START_REF] Berglund | Sound source localisation through active audition[END_REF], and for endfire AoAs (close to 0 • or 180 • ) which do not suffer from front-back confusion. Nevertheless the posterior probability of correct localization (not shown in the figure) decreases from broadside to endfire.

Updating the grid

The grid (1) is initialized with uniform probability. It is then recursively updated after each new measurement using Bayes law as

P (s = i|m 1:T , p 1:T ) = P (m T |s = i, p T )P (s = i|m 1:T -1 , p 1:T -1 ) i P (m T |s = i , p T )P (s = i |m 1:T -1 , p 1:T -1 ) (5)

ROBOT MOTION CONTROL

The mapping algorithm in Section 2 makes it possible to estimate the 3D coordinates of the source to a certain extent. We now turn to the question of controlling the robot motion so as to estimate it to the best extent possible in a given time. Bottom: likelihood of correct measurement up to tolerance .

p(-ε ≤ m t -α it ≤ ε | d it =2 m)

Cost function

For a given set of measurements, we quantify the amount of information carried by the grid by its entropy:

H(s|m 1:T , p 1:T ) = - i P (s = i|m 1:T , p 1:T ) log P (s = i|m 1:T , p 1:T ). ( 6
)
The lower the entropy, the greater the amount of information.

Let us assume that the robot has taken measurements up to a certain time T and that we wish to move it to a new pose at time T + 1. In order to find the optimal pose at time T + 1, we need to compute the entropy conditionally to the motion sequence p T +1:T +K for all possible motion sequences up to a fixed horizon T + K. This entropy cannot be deterministically computed, since future measurements m T +1:T +K are unavailable, but its expectation can be expressed as E m T +1:T +K [H(s|m 1:T +K , p 1:T +K )] = m T +1:T +K P (m T +1:T +K |m 1:T , p 1:T +K )H(s|m 1:T +K , p 1:T +K ) [START_REF] Nakadai | Active audition for humanoid[END_REF] with

P (m T +1:T +K |m 1:T , p 1:T +K ) = i P (s = i|m 1:T , p 1:T ) K k=1 P (m T +k |s = i, p T +k ). (8)
The expression in [START_REF] Nakadai | Active audition for humanoid[END_REF] is intractable due to combinatorial explosion of the set of future measurements with increasing K.

In order to address this issue, we compute the expectation of the entropy separately for each future pose p T +k instead:

c p T +k = E m T +k [H(s|m 1:T , m T +k , p 1:T , p T +k ) (9)
= m T +k P (m T +k |m 1:T , p 1:T , p T +k ) H(s|m 1:T , m T +k , p 1:T , p T +k ) (10)
where

P (m T +k |m 1:T , p 1:T , p T +k ) = i P (s = i|m 1:T , p 1:T )P (m T +k |s = i, p T +k ) (11)
and H(s|m 1:T , m T +k , p 1:T , p T +k ) is obtained via ( 5) and [START_REF] Blandin | Multi-source TDOA estimation in reverberant audio using angular spectra and clustering[END_REF] for each m T +k and p 1:T . We interpret the quantity (9) as the cost of a future pose and we denote it as c p T +k .

Dynamic programming algorithm

We assume that the cost of moving from one pose p T +k to the next p T +k+1 is 0 when this motion is feasible and +∞ otherwise. The optimal motion sequence is given by:

p T +1:T +K = min p T +1:T +K K k=1 c p T +k . (12)
This is a conventional dynamic programming problem. In the terminology of hidden Markov models, the cost of each pose is the log-observation probability and the cost of each move is the log-transition probability. The optimal sequence is obtained via the Viterbi algorithm (a.k.a. the Bellman algorithm in robotics). Once it has been found, the robot moves to the optimal pose p T +1 , takes a new measurement m T +1 , reestimates the optimal sequence up to T + K + 1, and so on.

EXPERIMENTAL EVALUATION

Protocol

We evaluated our position mapping and motion control algorithms for the localization of a speech source in the smart home at Inria Nancy. A Turtlebot equipped with a Kinect is placed at a fixed initial pose. A speech signal is emitted by a small loudspeaker at the same height and at a given distance and angle from the robot. Two different distances (1.2 and 2.4 m) and three angles (45 • , 90 • and 135 • ) are tested.

The room is discretized into a grid with 5 cm resolution. At each iteration, the robot moves to a new position at 30 cm distance and obtains one AoA estimate using HARK [START_REF] Nakadai | Design and implementation of robot audition system 'HARK' -open source software for listening to three simultaneous speakers[END_REF]. The robot orientation is constrained by the direction of movement from one position to the next. The actual robot pose after movement is measured with a Sokuiki laser. The source position is eventually estimated as the point with maximum probability in the grid. We implemented the proposed algorithms in C++ and interfaced them with HARK and the robot actuators using ROS.

Results

Fig. 2 illustrates one test case. After a first measurement, the optimal strategy is not to move towards the estimated source position, which addresses front-back confusion only, but to move in a position slightly aside of it, which makes it possible to estimate the source distance too. The robot eventually passes by the source and moves around it. A similar trajectory was followed in the other test cases.

Table 1 compares the proposed motion control strategy with a random motion strategy, where the next position is chosen randomly among all positions at 30 cm distance. The proposed strategy achieves consistently smaller localization error, down to 0.07 m after 8 to 12 measurements, compared to 0. 1. Average localization error (m) as a function of the initial distance to the source and the number of robot poses.

CONCLUSION

We proposed a motion control strategy for audio source localization by a mobile robot based on discrete grid mapping with a forward sensor model and on the estimation of the grid entropy resulting from each possible movement. We showed that the optimal motion reduces the average localization error up to a factor of 4 compared to random motion. Future work will focus on the design of a simultaneous localization and mapping (SLAM) algorithm applicable when the robot pose is unknown, on the handling of moving or intermittent sources, and on audiovisual integration.

ACKNOWLEDGMENTS

This work was supported by the AEN Inria PAL project.

Fig. 1 .

 1 Fig. 1. Top and middle: example sensor model distributions. Bottom: likelihood of correct measurement up to tolerance .

P(s = i | m 1 Fig. 2 .

 12 Fig. 2. Top: Source location probability after a first measurement at position p 1 (θ 1 = 0). Middle: Expected entropy at the next position p 2 (θ 2 = 0). Bottom: Complete robot trajectory. Initial and successive robot positions are shown as red squares, and the true source position as a green circle.

  [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF] to 0.29 m for the random strategy.

	Distance (m)	1.2		2.4
	Motion	random proposed random proposed
	2 poses	0.45	0.33	0.91	0.71
	4 poses	0.36	0.25	0.77	0.49
	6 poses	0.27	0.13	0.62	0.32
	8 poses	0.19	0.07	0.48	0.19
	10 poses	N/A	N/A	0.36	0.13
	12 poses	N/A	N/A	0.29	0.07
	Table

http://www.turtlebot.com/

http://www.loria.fr/~evincent/Roomsimove.zip


Sébastien Ménigot
email: sebastien.menigot@univ-tours.fr

Jean-Marc Girault
email: jean-marc.girault@univ-tours.fr

Harmonic Magnification by Time Reversal based on a Hammerstein Decomposition

Introduction

The medical ultrasound imaging systems have been improved by taking into account the nonlinear wave propagation. These improvements leads to increase the signal-to-noise ratio and the focusing in tissue harmonic imaging (THI). Since the backscattered nonlinearities are a function of the transmitted signal, enhancing these nonlinearities means to design the best wave [START_REF] Ménigot | Optimization of Contrast Resolution by Genetic Algorithm in Ultrasound Tissue Harmonic Imaging[END_REF]. Time reversal process is well-known to optimize the SNR and free itself to phase aberration, by combining a waveform design in space and in time thanks a physical matched filter [START_REF] Fink | Time-reversal of Ultrasonic Fields. 1. Basic Principles[END_REF]. However, it is not well-adjusted for THI, because the propagation of the time reversed signal destroys the harmonics.

How can we optimize the SNR in tissue harmonic imaging ?

As we want to guarantee a SNR optimization with a good focusing, the method has to include time reversal The solution firstly consists in extracting the harmonic component at 2 f 0 by a Hammerstein filter. Then this time reversed harmonic component at 2 f 0 is frequency shifted to the fundamental component at f 0 and retropropagated in the medium.

Methods

1. Sending a first standard excitation x standard (n) 2. Harmonic extraction based on a Hammerstein model where the nonlinear function is a frequency shifting by modulation (since it is a bijective function on R):

ŷ(n) = 2 ∑ p=1 M ∑ m=1 h p (m)x standard (n -m) • C p (n), (1)
n the discrete time, M the memory of the Hammerstein model and

C p (n) = cos 2π (p-1) f 0 n F s
with F s the sample frequency. The model can be solved by a pseudo-inversion:

h = X T X -1 X T y. ( 2
) h = [h 1 (1), . . . , h 1 (M), h 2 (1), . . . , h 2 (M)], ŷ = [y(M + 1), . . . , y(N)] T with N the sample quantity, the matrix of input signals X = [x 1 x 2 ] with x p =       v p (M) v p (M + 1) . . . v p (N) v p (M + 1) v p (M + 2) 0 . . . . . . v p (N) 0 0       , (3)
and v p (n) = x(n)C p (n). Finally, the second harmonic signal is:

y NL = x T 2 [h 2 (1), . . . , h 2 (M)] T .

Time reversing y NL

4. Annihilation of the second harmonic effects of the second harmonics:

x opt (n) = A • y NL (n) • cos 2π f 0 n F s + φ , (4)
A preserve the transmit power to x standard (n) and φ = arctan ϕ(y NL ) ϕ(cos) , with ϕ(y NL ) the phase of the the second harmonic signal y NL and the ϕ(cos) the phase of a cosinus at the frequency f 0.

Discussions and Conclusion

• Magnify the second harmonics and improve the SNR

• Hammerstein model using frequency shifting with the good phase → Annihilation of the harmonic effects before sending the time-reversed signal • Matched filter for harmonic components → Extension of the time reversal principle to second harmonics

  Program example of a Hammertein decomposition on Runmycode.orgSimulation• 2D nonlinear wave propagation in a cavity[START_REF] Anderson | A 2D Nonlinear Wave Propagation Solver Written in Open-Source MATLAB Code[END_REF] • 8-element probe centred at f c = 4 MHz • Initial Gaussian pulse x standard (n) centred at f 0 = 22 MHz and with a bandwidth of 50%