Maxence Cassier

Trent Degiovanni

S Ébastien Guenneau

Fernando Guevara Vasquez

ACTIVE EXTERIOR CLOAKING FOR THE 2D HELMHOLTZ EQUATION WITH COMPLEX WAVENUMBERS AND APPLICATION TO THERMAL CLOAKING

Keywords: 2010 Mathematics Subject Classification. 35J05, 31B10, 35K05, 65M80 Helmholtz equation, Heat equation, Active cloaking, Potential theory, Green identities

We design sources for the two-dimensional Helmholtz equation that can cloak an object by cancelling out the incident field in a region, without the sources completely surrounding the object to hide. As in previous work for real positive wavenumbers, the sources are also determined by the Green identities. The novelty is that we prove that the same approach works for complex wavenumbers which makes it applicable to a variety of media, including media with dispersion, loss and gain. Furthermore, by deriving bounds on Graf's addition formulas with complex arguments, we obtain new estimates that allow to quantify the quality of the cloaking effect. We illustrate our results by applying them to achieve active exterior cloaking for the heat equation.

Introduction

Our goal is to use specially designed sources to cloak or hide a bounded object from a probing field u i satisfying the two dimensional Helmholtz equation [START_REF] Martínez | Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements[END_REF] ∆u i + k 2 u i = 0, in a region containing the object. Here ∆ denotes the Laplacian. This is called active cloaking since to build the cloak, we use sources rather than passive materials that may be hard to manufacture [START_REF] Kadic | Metamaterials beyond electromagnetism[END_REF]. Moreover, a great advantage of the cloaking strategy we present is that it does not require completely surrounding the object to hide it, hence the exterior cloaking name. This idea was introduced in [START_REF] Vasquez | Active exterior cloaking for the 2d Laplace and Helmholtz equations[END_REF] for the 2D Helmholtz equation with k real (lossless propagative media). Here we allow k to take any values in the complex plane, except for the negative real axis. Thus using a frequency decomposition of the transient regime via a Fourier-Laplace transform on the time variable, our approach applies to cloaking objects for acoustic waves propagating in passive, dissipative, active, or dispersive media, and similarly for diffusive media. Interestingly, complex wavenumbers open a path to exterior cloaking for problems modelled by partial differential equations with second order derivatives in space and with time derivatives of an arbitrary order. Moreover, we derive new error estimates on the convergence of active exterior cloaking and apply our results to cloaking for the heat equation in the transient regime.

∂Ω 1 ∂Ω 2 ∂Ω 3 ∂Ω 4 x 1 x 2 x 3 x 4 Figure 1.
Active exterior cloak for the Helmholtz equation starting from the Green identities applied on the surface ∂Ω of a domain Ω. To hide the object (kite) inside Ω, the sources in the portions ∂Ω j of ∂Ω are moved to new locations x j , j = 1, . . . , N dev . The gray disks are the domains of divergence R j for the new fields, see [START_REF] Bellis | Simulating transient wave phenomena in acoustic metamaterials using auxiliary fields[END_REF]. For illustration purposes we took N dev = 4, but 3 sources would suffice to achieve exterior cloaking.

1.1. Active exterior cloaking. From potential theory [START_REF] Kellogg | Foundations of potential theory[END_REF] or using the Green identities (see e.g. [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]), it is possible to reproduce a solution to the homogeneous Helmholtz equation inside of a bounded region Ω and getting, simultaneously, a zero field outside of Ω. This is achieved by a distribution of monopole and dipole sources on the boundary ∂Ω that can be expressed in terms of the value of the field and its normal derivative on ∂Ω. As observed by Miller [START_REF] Miller | On perfect cloaking[END_REF], this principle can be used for cloaking. Indeed the monopole and dipole distribution can be chosen to generate the cloak field

(2) u c = -u i in Ω 0 outside Ω,
where X denotes the closure of a set X. We see that by linearity, u c + u i cancels out inside Ω without affecting u i outside of Ω. The end result is that objects inside Ω will not scatter and it is impossible to detect the cloaked field outside of Ω. We call this approach the Green identity cloak. A first drawback of this approach is that the probing field u i needs to be known ahead of time. A second drawback is that the sources completely surround the object that we wish to hide. The exterior cloaking approach lifts this second limitation.

To achieve exterior cloaking we follow the approach in [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] for the 2D Helmholtz equation with real wavenumbers, see also [START_REF] Vasquez | Transformation elastodynamics and active exterior cloaking[END_REF] for the 3D Helmholtz equation and [START_REF] Norris | Active elastodynamic cloaking[END_REF][START_REF] O'neill | Active cloaking of inclusions for flexural waves in thin elastic plates[END_REF] for elasticity. See also [START_REF] Onofrei | Active manipulation of fields modeled by the Helmholtz equation[END_REF] for a general analysis. The key observation is that Graf's addition formulas (see e.g. §10.23 in [21] or [START_REF] Martin | Multiple scattering, volume 107 of Encyclopedia of Mathematics and its Applications[END_REF]) can be used to move a monopole (or dipole) located at y to a new location x j . However the price to pay is that the new source is obtained by an infinite superposition of multipolar sources that diverges in the disk Assuming the portions ∂Ω j cover ∂Ω and their intersection is reduced to points, we can achieve the same cloaking effect as the Green identity cloak if we are outside of the region R 1 ∪ R 2 ∪ . . . ∪ R N dev . Here we extend this approach to complex wavenumbers which enables new applications. Moreover, by obtaining bounds for the Graf's addition formula with complex arguments, we derive a simple geometric series ansatz to predict the truncation error for the field generated by the multipolar sources. This extends previous work on the truncation error of Graf's addition formulae [START_REF] Cassier | Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the Foldy-Lax model[END_REF][START_REF] Meng | Bounds for truncation errors of graf's and neumann's addition theorems[END_REF] for real arguments. The truncation estimates allow us to quantitatively predict the quality of the cloaking effect.

1.2. Extending active exterior cloaking to complex k. Many partial differential equations in the frequency domain lead to the Helmholtz equation with a complex wavenumber. To name a few: the telegraph equation, the diffusion equation, Schrödinger equation and the Klein-Gordon equation (see e.g. [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF][START_REF] Ockendon | Applied Partial Differential Equations[END_REF] or [36, §1.1.2]). More generally, consider partial differential equations (PDEs) in the time domain of the form ( 5)

P (∂ t )u = ∆u + f,
where P is a polynomial of degree n and f(x, t) is a source term. Since ( 5) is a constant coefficient PDE, it admits a solution in the distributional sense, for example for any compactly supported source term f. This can be seen from the Malgrange-Ehrenpreis theorem, though without uniqueness or causality guarantees, see e.g. [START_REF] Friedlander | Introduction to the theory of distributions[END_REF].

Many classic equations are of the form [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF]. For example, the wave equation can be obtained with P (z) = z 2 , and the heat equation with P (z) = z. For a causal source f (i.e. f(x, t) = 0 for t < 0), we can analyze equations of the form (5) in the frequency domain by means of the Fourier-Laplace transform

(6) u(x, ω) = ∞ 0 dt[e iωt u(x, t)],
where ω is in general complex and u is assumed to grow sufficiently slowly. For example, we may assume that u(x, t) satisfies for t > 0

(7) u(•, t) ≤ Ce αt (1 + t p ),
where C > 0, p ∈ N, and α ∈ R are constants, see e.g. [START_REF] Arendt | Vector-valued Laplace transforms and Cauchy problems[END_REF][START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF]. Under this assumption, the Fourier-Laplace transform is defined on the half plane

(8) C + α := {ω ∈ C | Im(ω) > α}.
The choice of norm depends on the spatial differential operator. Since we focus on the Laplacian, we use the H1 norm on any bounded open set of interest (α, p, C could depend on the choice of the set). To summarize, in our situation we may assume that u satisfies the growth condition [START_REF] Cassier | Active thermal cloaking and mimicking[END_REF]

and u ∈ L 1 loc ((0, ∞), H 1 loc (R 2 ))
Im(k) R e ( k ) = I m ( k ) R e ( k ) = - I m ( k ) Figure 2.
Diagram of cases for the wavenumber k defined by k = ±i P (-iω), where we used the principal square root and the sign is chosen to match the sign of Im(k). and all its time derivatives up to order n (the degree of the polynomial P ) and the source term

f are in L 1 loc ((0, ∞), L 2 loc (R 2
)) and satisfy also [START_REF] Cassier | Active thermal cloaking and mimicking[END_REF]. This growth condition allows to make sense of the Fourier-Laplace transform for solutions that may grow exponentially in time 2 , as is the case of active media 3 .

Assuming u and all its time derivatives up to order n -1 vanish at t = 0, we get that u(x, ω) satisfies the Helmholtz equation: [START_REF] Cassier | Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the Foldy-Lax model[END_REF] ∆u + k 2 u = -f, with the relation k 2 = -P (-iω). This formalism shows that k may be complex e.g. in the case of the heat equation. Another example is the modified wave equation with P (z) = z 2 + αz, α ∈ R. If α > 0, this is the dissipative wave equation which corresponds to wave propagation in lossy media. If ω > 0, we can choose the root k such that Im(k) > 0 to get spatially decreasing solutions to (9) 4 . Whereas with α < 0, ω > 0 we choose k such that Im(k) < 0 to get spatially increasing solutions corresponding to an amplifying medium (medium with gain), see e.g. [START_REF] Dyatlov | Mathematical theory of scattering resonances[END_REF]. We summarize the possibilities in fig. 2. We have highlighted the region with Re(k 2 ) < 0 (or equivalently |Re(k)| < |Im(k)|), since a form of the maximum principle holds for the Helmholtz equation with such k [START_REF] Kresin | Criteria for validity of the maximum modulus principle for solutions of linear strongly elliptic second order systems[END_REF][START_REF] Krutitskii | The Dirichlet problem for the dissipative Helmholtz equation in a plane domain bounded by closed and open curves[END_REF]. Later in section 3.2, we see that the maximum principle gives a form of stability for the accuracy of approximations to our cloaking approach.

Other situations where complex wavenumber k arises are in passive, dispersive media, where the index of refraction is a complex valued function of frequency ω 2 The growth condition [START_REF] Bretin | On the green function in visco-elastic media obeying a frequency power-law[END_REF] and the regularity assumption u ∈ L 1 loc ((0, ∞), L 2 loc (R 2 )) ensure in particular the existence of the Laplace transform (6) as a Bochner integral with respect to t valued in L 2 [START_REF] Arendt | Vector-valued Laplace transforms and Cauchy problems[END_REF], and thus also pointwise for a. e. x ∈ R 2 and all ω ∈ C + α . We point out that in a more general context, the Fourier-Laplace transform can be extended to spaces of distributions, see e.g. [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF][START_REF] Sayas | Retarded potentials and time domain boundary integral equations[END_REF][START_REF] Zemanian | Realizability theory for continuous linear systems[END_REF].

3

By "active media" we mean there is energy input that may lead to increase of the magnitude of the fields in time, see for instance [START_REF] Skaar | Fresnel equations and the refractive index of active media[END_REF]. By "gain media" we mean that there are spatially growing outgoing solutions of the Helmholtz equation as for e.g. resonant states, see [START_REF] Dyatlov | Mathematical theory of scattering resonances[END_REF]. This is not a universal nomenclature. 4 We shall see in section 2 that the decay is consistent with the choice of Green function [START_REF] Cassier | Mathematical models for dispersive electromagnetic waves: an overview[END_REF]. [START_REF] Tip | Linear absorptive dielectrics[END_REF][START_REF] Figotin | Spectral theory of time dispersive and dissipative systems[END_REF][START_REF] Gralak | Macroscopic maxwell's equations and negative index materials[END_REF][START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF][START_REF] Cassier | Mathematical models for dispersive electromagnetic waves: an overview[END_REF][START_REF] Bellis | Simulating transient wave phenomena in acoustic metamaterials using auxiliary fields[END_REF]. A typical example is the dispersion law given by the Drude-Lorentz model in electromagnetism, which can model both metals and metamaterials with negative index of refraction [START_REF] Ramakrishna | Physics of negative refractive index materials[END_REF]. In acoustics (see for e.g. section 1.1.2 of [START_REF] Gumerov | Fast multipole methods for the Helmholtz equation in three dimensions[END_REF]), complex wavenumbers also arise when studying acoustic (pressure) waves propagating within complex (but homogeneous and isotropic) fluids (i.e. not barotropic) with a relaxation time (due to the presence of solid particles or bubbles) which can be modelled in the time-harmonic regime with

k 2 = ω 2 /(c 2 (1 -iωτ p )).
Here c is the speed of sound (m/s), ω the pulsation frequency (rad/s) and τ p the density relaxation time (s). We note that the frequency domain formulation gives a strategy for active exterior cloaking in the time domain for equations of the form [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF], and even when the powers are fractional or negative (which corresponds to integro-differential equations). One caveat of our approach is that the sources need to make sense physically in [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF]. However, for active thermal cloaking (section 4) the sources we obtain can be thought of as Peltier devices [START_REF] Disalvo | Thermoelectric cooling and power generation[END_REF].

1.3. Structure of the paper. We derive convergence estimates for the Graf addition formula applied to Green functions in section 2, showing that the truncation error of the series can be dominated by that of a geometric series with ratio that depends only on the position of the evaluation point relative to the positions of the original and new sources. This result can be applied to get truncation estimates for the multipolar source expansions that appear in active exterior cloaking for possibly complex wavenumbers (section 3). In section 3, we also use a form of the maximum principle for the Helmholtz equation, which guarantees the truncation errors in a region are maximum on the boundary of the region (this only holds for a class of dissipative or diffusive media). The time domain problem for the heat equation is then considered in section 4. We conclude with some future work and perspectives in section 5.

Moving sources

The field evaluated at x corresponding to a point source located at y is given by the appropriate Green function [START_REF] Cassier | Mathematical models for dispersive electromagnetic waves: an overview[END_REF] G(xy; k) = i 4 H

(1)

0 (k|x -y|),
where H

(1) 0

is the zero-th order Hankel function of the first kind 5 , see e.g. [21, eq. 10.4.3]. Moreover G(x; k) → 0 as |x| → ∞ whenever Im(k) ≥ 0, as can be seen from the large argument asymptotic for Bessel functions [21, eq. 10.2.5]. When Im(k) < 0, the same asymptotic shows that |G(x; k)| → ∞. Thus the choice of Green function is consistent with the loss and gain conventions in the diagram appearing in fig. 2. 5 In the context of the wave equation with constant propagation speed c, we have k = ω/c > 0 and the choice of Green function [START_REF] Cassier | Mathematical models for dispersive electromagnetic waves: an overview[END_REF] corresponds to outgoing waves. This is consistent with the Fourier-Laplace transform convention [START_REF] Bretin | On the green function in visco-elastic media obeying a frequency power-law[END_REF] and with the convention that the corresponding time harmonic field is Re(exp[-iωt]G(x -y; k)). In fact Re(k) > 0 also gives outgoing waves, as can be seen e.g. from adapting the discussion [36, eq. 1.2.12] from 3D to 2D using [21, eq. 10.4.3].

Thanks to the Graf addition formulas [21, eq. 10.23.7], we can move (with three significant caveats) the source from location y to another location x j , indeed:

(11) G(x -y; k) = i 4 ∞ m=-∞ H (1) m (k|x -x j |)J m (k|y -x j |) exp[imθ],
where θ = arg(x-x j )-arg(y -x j ) and arg x is the counter-clockwise angle between the vectors x and (1, 0). Here J m is the m-th order Bessel function of the first kind, see e.g. [21, eq. 10.2.2]. The first caveat is that the new source in [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] is no longer a monopole point source like [START_REF] Cassier | Mathematical models for dispersive electromagnetic waves: an overview[END_REF], but a linear combination of Helmholtz equation solutions that diverge as |xx j | → 0, of the form V m (xx j ) where [START_REF] Cohen | Numerical methods for Laplace transform inversion[END_REF] V m (x) = exp[im arg(x)]H (1) m (k|x|), and that are known as multipolar sources (or cylindrical outgoing waves when k is real). The second caveat is that the Graf addition formula is only valid for k ∈ C \ (-∞, 0]. The third caveat is that the series converges only outside of the disk D xj ,y , as defined in [START_REF] Arendt | Vector-valued Laplace transforms and Cauchy problems[END_REF].

The same method and caveats apply if we desire to move a dipole located at y and oriented in the direction ν(y) normal to the boundary ∂Ω at y, or more precisely [START_REF] Colton | Integral equation methods in scattering theory[END_REF] ∂G ∂ν(y)

(x -y; k) = i 4 ∞ m=-∞ H (1) m (k|x -x j |) ∂ ∂ν(y) J m (k|y -x j |) exp[imθ] ,
where θ is the same as in [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF]. Formally speaking, equation ( 13) can be obtained by taking the gradient term by term (with respect to y) in [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] and then taking the dot product with ν(y). The differentiation term by term can be easily justified by using Lemma 1 in order to prove that the involved series of gradients is locally normally convergent (and thus locally uniformly convergent) with respect to y when x, x j are fixed. We start in section 2.1 by proving convergence estimates for [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] and [START_REF] Colton | Integral equation methods in scattering theory[END_REF]. The convergence errors are illustrated numerically in section 2.2.

Truncation error estimates.

To study the convergence rate of [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] and [START_REF] Colton | Integral equation methods in scattering theory[END_REF] we define the truncation to 2M + 1 terms of the formula [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] for moving the point source at x j to location y by ( 14)

G j,M (x -y; k) = i 4 M m=-M H (1) m (k|x -x j |)J m (k|y -x j |) exp[imθ],
where M is an integer. The truncation error for a monopole and dipole are given respectively by

(15) R j,M (x; k) = |G(x -y; k) -G j,M (x -y; k)| and R j,M (x; k) = ∂ ∂ν(y) [G(x -y; k) -G j,M (x -y; k)] .
In the next theorem, we show that these truncation errors are dominated by the truncation errors of well-known series such as geometric series. Our convergence estimates account for moving sources from different original positions y to a single new position x j . The case of different original positions y is useful in the context of active exterior cloaking (section 3). We point out that the monopole truncation error was derived for real wavenumbers by [START_REF] Cassier | Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the Foldy-Lax model[END_REF]Lemma 9] and [START_REF] Meng | Bounds for truncation errors of graf's and neumann's addition theorems[END_REF], using techniques that are similar to the ones we use here. Theorem 1 applies to the monopole and the dipole truncation errors and allows for estimates that are uniform with respect to original source location y, evaluation point x and complex wavenumbers k.

Theorem 1. Let x j ∈ R 2 , M ≥ 2, Y ⊂ R 2 be a compact set and define the disk

D max j = {x ∈ R 2 |x -x j | ≤ max y∈Y |y -x j |}. Let X ⊂ R 2 \ D max j
and K ⊂ C \ (-∞, 0] be compact sets. Then for any (x, k) ∈ X × K we have the following bounds for the monopole and dipole truncation errors

R j,M (x; k) ≤ C 1 -ln(1 -a x ) - M m=1 a m x m R j,M (x; k) ≤ C 2 a M +1 x 1 -a x , (16)
where C 1 and C 2 may depend on X, Y , K and

a x := max y∈Y |y -x j | |x -x j | . (17)
To prove theorem 1 we need the following asymptotic formulas for Bessel functions that are uniform on the order and are valid on appropriate compact sets of C excluding the negative real axis (-∞, 0]. This is because we use the power series definitions for Bessel functions in [21, §10.8] and the power series expansion for H [START_REF] Martínez | Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements[END_REF] n (z) is not valid for z ∈ (-∞, 0] (as the expansion contains the term (2/π) ln(z/2)J n (z) which has a discontinuity for such z). Lemma 1. Let K 1 be a compact set of C, K 2 be a compact set of C \ (-∞, 0] and n ∈ N. Then there exists constant C K1 , CK1 and CK2 (independent of n) such that:

J n (z) - 1 n! z 2 n ≤ C K1 (n + 1)! |z| 2 n+2 , ∀z ∈ K 1 , n ≥ 0, (18)
J n (z) - 1 2(n -1)! z 2 n-1 ≤ CK1 n! |z| 2 n+1 , ∀z ∈ K 1 , n ≥ 1, (19)
H (1) n (z) - i (n -1)! π 2 z n ≤ CK2 (n -2)! 2 |z| n-2 , ∀z ∈ K 2 , n ≥ 2. ( 20
)
The proof of Lemma 1 is included in appendix A. From Lemma 1, we can deduce the following inequalities for J n , J n and H [START_REF] Martínez | Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements[END_REF] n that are useful in the proof of theorem 1. Let K 1 be a compact set of C and K 2 be a compact set of C \ (-∞, 0]. Applying the inequality [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF], we get that for all z ∈ K 1 and n ≥ 0

|J n (z)| ≤ J n (z) - 1 n! z 2 n | + 1 n! |z| 2 n ≤ C K1 (n + 1)! |z| 2 n+2 + 1 n! |z| 2 n ≤ B K1 n! |z| 2 n with B K1 = max 1, C K1 max z∈K1 |z| 2 2 > 0. (21)
Similarly, one deduces from formula [START_REF] Diatta | Non-singular cloaks allow mimesis[END_REF] and ( 20) that there exists two constants BK1 > 0 and BK2 > 0 such that:

J n (z) ≤ BK1 (n -1)! |z| 2 n-1 , ∀z ∈ K 1 , n ≥ 1, (22)
H (1) n (z) ≤ BK2 (n -1)! 2 |z| n , ∀z ∈ K 2 , n ≥ 2. ( 23
)
We remark that Lemma 1 shows that the inequalities ( 21), ( 22) and ( 23) are optimal in the sense that they bound the functions by their leading order term. We are now ready to prove theorem 1.

Proof.

Step 1: inequality on the monopole truncation error R j,M . We want to apply the Lemma 1 to bound the terms J m (k|y-x j |) and H

(1)

m (k|x-x j |) for m ≥ M + 1 appearing in the expression of R j,M = |G(x -y; k) -G M (x -y; k)|. Noting first that H (1) -m = (-1) m H (1)
m and J -m = (-1) m J m , we obtain that:

R j,M (x; k) ≤ |m|≥M +1
|H (1) m (k|x -

x j |)J m (k|y -x j |)| = ∞ m=M +1 2|H (1) m (k|x -x j |)J m (k|y -x j |)|.
Thus, applying inequalities ( 21) and [START_REF] Dyatlov | Mathematical theory of scattering resonances[END_REF] gives that there exists C 1 > 0 (depending on the compacts X and K but not on the truncation index M ) such that

R j,M (x; k) ≤ C 1 ∞ m=M +1 1 m a m x = C 1 -ln(1 -a x ) - M m=1 a m x m .
Step 2: inequality on the dipole truncation error R j,M . For R j,M (x; k), we have the following (note that θ and n depend on y)

R j,M (x; k) = |m|≥M +1 H (1) m (k|x -x j |) ∂ ∂ν J m (k|y -x j |) exp[imθ] = |m|≥M +1 H (1) m (k|x -x j |)exp[imθ]k (y -x j ) • ν |y -x j | 3 J m (k|y -x j |) + H (1) m (k|x -x j |)exp[imθ]J m (k|y -x j |)(-im) ∂ ∂ν arg(y -x j ) ≤ |m|≥M +1 |H (1) m (k|x -x j |)k (y -x j ) • ν |y -x j | J m (k|y -x j |)| + |H (1) m (k|x -x j |)J m (k|y -x j |)m ∂ ∂ν arg(y -x j )|.
Since H

(1)

-m = (-1) m H (1)
m and J -m = (-1) m J m we can reduce the sum to

R j,M (x; k) ≤ 2 ∞ m=M +1 |H (1) m (k|x -x j |)k (y -x j ) • ν |y -x j | J m (k|y -x j |)| + 2 ∞ m=M +1 |H (1) m (k|x -x j |)J m (k|y -x j |)m ∂ ∂ν arg(y -x j )|. ( 24
)
We deal with the two latter sums separately. We start with the second sum due to its similarity to the monopole error. Noting that

(25) ∂ ∂ν arg(y -x j ) = (y -x j ) ⊥ |y -x j | 2 • ν ≤ 1 |y -x j | ≤ C
where for a vector u = (u 1 , u 2 ) ∈ R 2 , u ⊥ = (-u 2 , u 1 ) and the positive constant C is defined by C = max y∈Y |yx j | -1 . Thus, combining ( 21), ( 23) and [START_REF] Farhat | Cloaking through cancellation of diffusive wave scattering[END_REF] gives that ( 26)

∞ m=M +1 m|H (1) m (k|x -x j |)J m (k|y -x j |) ∂ ∂ν arg(y -x j )| ≤ C 3 ∞ m=M +1 a m x ,
where the positive constant C 3 depends only on Y , X and K. Now, we estimate the first sum of [START_REF] Farhat | Scattering theory and cancellation of gravity-flexural waves of floating plates[END_REF]. By virtue of the estimates ( 22) and ( 23), one gets that there exists a constant

C 4 > 0 such that ∞ m=M +1 |H (1) m (k|x -x j |)k (y -x j ) • ν |y -x j | J m (k|y -x j |)| ≤ C 4 |x -x j | ∞ m=M +1 a m x . Setting C 5 = C 4 (max x∈X |x -x j |) -1 > 0, one obtains (27)
∞ m=M +1 |H (1) m (k|x -x j |)k (y -x j ) • ν |y -x j | J m (k|y -x j |)| ≤ C 5 ∞ m=M +1 a m x .
Combining ( 24), ( 26) and ( 27) yields the second inequality of (16).

2.2.

Numerical experiments for truncation error estimates. We illustrate our bounds numerically in the case where there is only one source to move, i.e. Y = {y}. The bounds in theorem 1 involve a quantity a x that can be estimated from the relative positions of x, y and x j (repectively, the evaluation point and the original and new source positions). The bounds also involve constants C 1 , C 2 that may depend in non-obvious ways on the different choices of compact sets in space and wavenumber. To estimate C 1 (resp. C 2 ) for a particular choice of (x, k) ∈ X × K, we assume the truncation error has the form predicted by the respective upper bound in ( 16), and we find the C 1 (resp. C 2 ) that matches the actual error explicitly for one small value of M . We repeat this estimate on a grid for X × K and then take the maxima of the estimates for C 1 (resp. C 2 ) over the grid.

In fig. 4 we show these bounds for {x} × K, for different choices of wavenumber sets in complex plane and for the fixed evaluation point x = (0, 0.43). The original source location is y = (0, 0) and it is moved to the new location x j = (0, 0.2). We took M = 4 to approximate the constants C 1 and C 2 over {x} × K and then use our estimated C 1 and C 2 to predict the truncations errors with M = 20 terms. The different wavenumber ranges in the complex plane that we considered are summarized in fig. 3. Then in fig. 5, we estimated the constants C 1 and C 2 on X × {k} for four different wavenumbers k ∈ C. Here X is the region

X = {x ∈ R 2 | 1/2 ≤ |y -x j |/|x -x j | ≤ 1}, i.

Condition Range

Color real 15) respectively, and the bounds from theorem 1 for a single point in space which was relatively close to the source. The color represents the set for k as given in fig. 3. For the dipoles we cannot differentiate the different bounds as they all lie on top of each other. We observe that the dipole errors and bounds are larger. This is in line with theorem 1: the dipole error bound decays slower than the monopole error bound.

K 1 = [0.5, 20] green imaginary K 2 = {iθ | θ ∈ K 1 } red dissipative K 3 = {(1 + √ 2i)θ/ √ 3 | θ ∈ K 1 } blue amplifying K 4 = {(99 -i √ 199θ)/100 | θ ∈ K 1 } pink

Active exterior cloaking at fixed frequency

One can achieve active cloaking by observing [START_REF] Miller | On perfect cloaking[END_REF] that a distribution of monopoles and dipoles on the boundary ∂Ω of a bounded open Ω ⊂ R 2 can create a field u c that cancels out the incident or probing field u i inside a region Ω, while vanishing outside, or in other words satisfying [START_REF] Alwakil | Inverse heat mimicking of given objects[END_REF]. By applying the Green identities (see e.g. [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]) or potential theory (see e.g. [START_REF] Kellogg | Foundations of potential theory[END_REF]) the function u c = -u i is given for x / ∈ ∂Ω

u c (x; k) = ∂Ω dS(y)[ -∂u i ∂ν(y) (y; k)G(x -y; k) + u i (y; k) ∂G ∂ν(y) (x -y; k)],
where G(x; k) is the Green function [START_REF] Cassier | Mathematical models for dispersive electromagnetic waves: an overview[END_REF].

∆u i + k 2 u i = 0 in O for k ∈ C \ {0}.
Then as -∆u i = k 2 u i on O, one easily proves by interior elliptic regularity of the minus Laplacian operator (applying iteratively Theorem 2 page 314 of [START_REF] Evans | Partial differential equations[END_REF])

that u i ∈ C ∞ (Ω).
We point out that u c = -u i ∈ C ∞ (Ω) and the outward normal vector ν(y) ∈ L ∞ (∂Ω) since ∂Ω is a Lipschitz boundary. Thus, it is clear that the Dirichlet trace u i is smooth on ∂Ω and that the Neumann trace

[∂u i /∂ν(y)](y; k) is in L ∞ (∂Ω).
Hence, the integrand in (28) is integrable as a sum of two products of L ∞ (∂Ω) functions. Indeed since we have x / ∈ ∂Ω, the Green function G(xy; k) is smooth for y ∈ ∂Ω and its normal derivative [∂G/∂ν(y)](xy; k) is in L ∞ (∂Ω) as a function of y.

To get exterior cloaking, the idea is to move the monopoles and dipoles on the portions ∂Ω j of the boundary ∂Ω to the new source locations x j . Formally, this can be done by replacing the Green function and its normal derivative in the representation formula [START_REF] Friedlander | Introduction to the theory of distributions[END_REF] by their series expansions [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] and [START_REF] Colton | Integral equation methods in scattering theory[END_REF]. Theorem 2 and Remark 1 allow to permute the order of the series and the integral over ∂Ω (since the series is normally convergent with respect to y). Thus we can express the new cloaking field as ( 29)

u e (x; k) = N dev j=1 ∞ m=-∞ b j,m V m (x -x j ; k),
where V m are multipolar sources [START_REF] Cohen | Numerical methods for Laplace transform inversion[END_REF] and the coefficients b j,m are given by [START_REF] Greenleaf | Anisotropic conductivities that cannot be detected by eit[END_REF] in terms of integrals over the ∂Ω j , identical to those obtained in [START_REF] Vasquez | Active exterior cloaking for the 2d Laplace and Helmholtz equations[END_REF]. We emphasize that theorem 2 is valid for complex k with the exception of the negative real axis, whereas the result in [START_REF] Vasquez | Active exterior cloaking for the 2d Laplace and Helmholtz equations[END_REF] is only proven for real k positive. Moreover, theorem 2 leverages on the Graf addition formula truncation error estimates in theorem 1, to give the truncation error when we consider instead the truncated fields:

(30) u (M ) e (x; k) = N dev j=1 M m=-M b j,m V m (x -x j ; k).
This error estimate is novel and applies to the results in [START_REF] Vasquez | Active exterior cloaking for the 2d Laplace and Helmholtz equations[END_REF].

Theorem 2. Let Ω ⊂ R 2 be
region R = R 1 ∪ • • • ∪ R N dev , i.e. the union of the disks R j in (4). Let K be a compact subset of C \ (-∞, 0] and X a compact subset of R 2 \ R. Define the coefficients b j,m in (29) and (30) by (31) b j,m = ∂Ωj dS(y) - ∂u i ∂ν(y) (y; k)U m (y -x j ; k) + u i (y; k) ∂U m (y -x j ; k) ∂ν (y)
,

where j = 1, .., N dev , m ∈ Z, and U m (x; k) = J m (k|x|) exp[-imarg(x)].
Then there exists a constant C > 0 (which may depend on K, X, u i and ∂u i /∂ν(y)) such that for any

(x, k) ∈ X × K, u c (x; k) -u (M ) e (x; k) ≤ C a M +1 1 -a where a = max x∈X max j=1,...,N dev max y∈∂Ωj |y -x j | |x -x j | < 1.
In particular for any x / ∈ R and k / ∈ (-∞, 0], we have u c (x; k) = u e (x; k).

Remark 2. The integral appearing in the definition of the b j,m in (31) can be expressed as a series when the incident field u i is given in terms of its cylindrical wave expansion [59, Theorem 2]. Although the series expansion is proven for the 2D Helmholtz equation with k > 0, we conjecture it is valid for complex k. Proof. We rewrite the boundary integral representation [START_REF] Friedlander | Introduction to the theory of distributions[END_REF] of the cloaking field u c as a sum of integrals over the portions ∂Ω j of the boundary. We then apply [START_REF] Cassier | Bounds on herglotz functions and fundamental limits of broadband passive quasistatic cloaking[END_REF] to yield

u c (x) = N dev j=1 ∂Ωj dS(y) -∂u i ∂ν(y) (y) i 4 ∞ m=-∞ V m (x -x j )U m (y -x j ) + u i (y) i 4 ∂ ∂ν(y) ∞ m=-∞ V m (x -x j )U m (y -x j ) (32)
which holds for x / ∈ R. We approximate the cloak field by u

(M ) e
with coefficients b j,m chosen as in [START_REF] Greenleaf | Anisotropic conductivities that cannot be detected by eit[END_REF] to match the |m| ≤ M terms in the series in [START_REF] Guenneau | Transformation thermodynamics: cloaking and concentrating heat flux[END_REF]. Thus the error we make by approximating u c (x) by u (M ) e at some x / ∈ R can be bounded by

|u c (x) -u (M ) e (x)| = N dev j=1 ∂Ωj dS(y) -∂u i ∂ν(y) (y) i 4 |m|≥M +1 V m (x -x j )U m (y -x j ) + u i (y) i 4 ∂ ∂ν(y) |m|≥M +1 V m (x -x j )U m (y -x j ) ≤ N dev j=1 ∂Ωj dS(y) -∂u i ∂ν(y) (y) i 4 |m|≥M +1 V m (x -x j )U m (y -x j ) + u i (y) i 4 ∂ ∂ν(y) |m|≥M +1 V m (x -x j )U m (y -x j ) .
We notice that

R j,M (x; k) = |m|≥M +1 V m (x -x j )U m (y -x j ) R j,M (x; k) = ∂ ∂ν(y) |m|≥M +1 V m (x -x j )U m (y -x j ) ,
allowing us to apply theorem 1 to bound the truncation by remainders of a geometric series. It follows that there exists two positive constants C 1,j and C 2,j that depend on the compact sets X, Y j = ∂Ω j and K such that

|u c (x) -u (M ) e (x)| ≤ N dev j=1 ∂Ωj dS(y) ∂u i ∂ν(y) (y) 1 4 C 1,j + u i (y) 1 4 C 2,j a M +1 x,j 1 -a x,j , with a x,j = max y∈∂Ω j |y -x j | |x -x j | . Setting C 1 = max j=1,...,N dev C 1,j and C 2 = max j=1,...,N dev C 2,j , it follows that |u c (x) -u (M ) e (x)| ≤ n max j |∂Ω j | C a M +1 1 -a ,
where from Remark 1 we have that u i ∈ C 0 (∂Ω) and ∂u i /∂ν(y) is L ∞ (∂Ω): and on the other hand that the maximum of |u| is only reached on the boundary ∂B of B. We note that the strong maximum principle is not valid for k outside the region [START_REF] Vasquez | Transformation elastodynamics and active exterior cloaking[END_REF] as one can find examples of solutions violating it [START_REF] Kresin | Criteria for validity of the maximum modulus principle for solutions of linear strongly elliptic second order systems[END_REF].

C = sup y∈∂Ω ∂u i ∂ν(y) (y) 1 4 C 1 + max y∈∂Ω u i (y) 1 4 C 2 , and (33)
In particular if u and v are smooth solutions to the homogeneous Helmholtz equation in B with k satisfying [START_REF] Vasquez | Transformation elastodynamics and active exterior cloaking[END_REF], the maximum of the error |u(x)v(x)| is attained only at the boundary ∂B. In other words, the error within the domain is controlled by the error on the boundary (the Dirichlet data). This can be viewed as a form of stability for the boundary integral representation [START_REF] Friedlander | Introduction to the theory of distributions[END_REF]. Moreover, u e and u Finally we point out that when numerically evaluating the boundary representation formula [START_REF] Friedlander | Introduction to the theory of distributions[END_REF], we use finitely many monopole and dipole sources on the domain ∂Ω. Following the same argument, the error we make with this discretization is also maximum on the boundary of any bounded domain B such that B ⊂ R 2 \ ∂Ω. We numerically illustrate in fig. 9 that the maximum principle predicts that the maximum cloaking errors occurs on the boundary of a region and not inside. 3.3.1. Evaluation of the cloak field. We illustrate theorem 2 numerically using a disk region Ω and n = 4 sources, as shown in fig. 6. While we chose to illustrate exterior cloaking with four multipolar sources, only three are necessary in two dimensions to give a non-empty region cloaking [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF]. Cloaking fields u (M ) e with M = 22 are shown in fig. 7 for several representative wavenumbers on the square [0, 10] 2 using a 200 × 200 uniform grid. The disk Ω is centered at [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF][START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF] and with radius δ C = 10/6. The x j are placed on a circle of radius δ D (see fig. 6), where δ D is chosen to maximize the area of the cloaking region for δ C fixed. The optimization is done via a simple geometric argument similar to [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] and gives δ D = 5 √ 2/3. The incident field is generated by a point source at y = (2, 5). To evaluate the truncated exterior cloaking field, u (M ) e , we use an equi-spaced discretization of ∂Ω, into points y i with i = 1, ..., n int . We split ∂Ω into n = 4 regions each associated with a new source location x j . We choose n int so that there is an equal number of discretization points of ∂Ω for each ∂Ω j and choose the x j such that max yi∈∂Ωj |y ix j | is equal for all j in order to keep the size of the theoretical divergence regions R j of our devices equal. The integrals over ∂Ω j that determine the coefficients b j,m in theorem 2 are approximated using the midpoint rule (so that the total integral over ∂Ω is the trapezoidal rule).

We note that the color scale in fig. 7 is deliberately limited to exclude the large fields near the new source locations x j which are due to the singularity of u (M ) e at the x j . This may seem as an impediment to physically realize such cloaking devices. However, as noted in [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF], it is possible to use the Green exterior representation formula (valid for Im(k) ≥ 0, see [START_REF] Colton | Integral equation methods in scattering theory[END_REF]Theorem 3.3]) to replace the multipolar sources by a distribution of monopoles and dipoles on some boundary enclosing each of the x j . Since the cloak field u (M ) e is smaller, we expect it is easier to realize in practice. The drawback is that these "extended cloaking devices" leave only small gaps between the cloaked region and the exterior. By theorem 2 we expect The abscissa corresponds to the parameter θ for the segment of wavenumbers K 3 , as defined in fig. 3.

that |u (M ) e (x)| → ∞ for x ∈ R, so increasing M leads to smaller gaps. So there is a tradeoff between getting larger gaps and approximating the ideal cloak field u e .

Computation of error bounds and the maximum principle.

In order to use the error bounds from theorem 2 in our particular geometric setup, we define the radius r Ci (resp. r C O ) of the largest (resp. smallest) inscribed (resp. circumscribed) circle that is inside (resp. outside) the divergence region R (defined in theorem 2). The inscribing and circumscribing circles are represented in fig. 6. We recall from (3.1) that the cloak field truncation error can be bounded by the truncation error of a geometric series with ratio a that is determined by the relative positions of the ∂Ω j , the x j and the region of interest where we want to evaluate the fields, see [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF]. Since we expect the cloaking fields to diverge close to R, it does not make sense to evaluate the errors on the inscribing and circumscribing circles. We do it instead on slightly smaller or larger circles of radii r C0 + 0.1δ C and r Ci -0.1δ C . If we take the region X from theorem 2 to be the union of these two circles, then a simple geometric argument yields that there are eight points in X that attain the maximum over X in the definition of a [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF]. At each of these points the ratio of the geometric series ansatz is the same, so we can conclude the truncation error can be bounded by C(1a M +1 )/(1a), where a < 1, but the constants C depend on the point. We first estimate the constant C at a point x by using the "empirical method" we used in section 2.2. In other words we find the C for which |u

(3) e (x) -u e (x)| is equal to C(1 -a M +1 )/(1 -a) with M = 3.
Then we take the worst case scenario, i.e. the largest of such C for the eight points in X that we considered. We emphasize that this is a heuristic meant to simplify the exhaustive method, where would have to evaluate the largest C for all x ∈ X. We summarize in fig. 8 the application of this heuristic for wavenumbers k ∈ K 3 (as defined in fig. 3). In these experiments we used 128 equispaced discretization points for ∂Ω, δ C = 10/6 and δ D = 5 √ 2/3. Finally the incident field we used for this experiment was a point source located at x = (8, 5). As can be seen from fig. 8, the error bound we obtain for M = 22 overestimates the actual error and follows the same trend for varying wavenumber.

We illustrate in fig. 9 that when |Im(k)| > |Re(k)|, the maximum principle (section 3.2) can be used to predict where the maximum cloaking error occurs. In fact the wavenumbers we used in for fig. 8 also allow us to use the maximum principle to observe that a bound for the truncation error on the boundary of the circle with radius r Ci -0.1δ C automatically leads to a bound on the whole disk of same radius.

3.3.3.

Calculating scattered fields for Im(k) ≥ 0. To demonstrate cloaking, we recall how to calculate scattered fields from a sound-soft (or homogeneous Dirichlet) obstacle A. Here we follow the discussion in [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]. We assume for simplicity that the obstacle A is a bounded domain with C 2 boundary ∂A (for similar results in the more general case of Lipschitz boundary see [52, §9]). The scattering problem problem can be posed as the following exterior Dirichlet problem

∆u s + k 2 u s = 0, x ∈ R 2 \ A u s = -u i , x ∈ ∂A,
where u s also satisfies the Sommerfeld radiation condition

lim |x|→∞ |x| 1 2 ∂ ∂|x| -ik u s (x) = 0
where ∂/∂|x| denotes the radial derivative and the limit is uniform for all directions x/|x| (see [14, §3.4]). The exterior Dirichlet problem has a unique solution

u s ∈ H 1 loc (R 2 \ A) for Im(k) ≥ 0 and u i | ∂Ω ∈ H 1/2 (∂A)
, see e.g. [14, §3.2]. This is clearly the case under the assumptions in Remark 1, since u i ∈ C ∞ (∂Ω).

We seek the scattered field in the form of a mixed single and double-layer potential ψ ∈ H

(z -y; k)ϕ(y) = 1 2 ϕ(x) + [Kϕ](x) , (37)
where z → x + denotes the limit from the exterior of A. Taking the limit of ( 36) as we approach the boundary of A from the exterior and applying (37) yields ( 38)

ψ + Kψ -iηSψ = -2u i | ∂A ,
which has a unique solution ψ (see e.g. [14, §3.2]). We assume that ∂A admits a 2π-periodic parametrization of the form

q(τ ) = (x 1 (τ ), x 2 (τ )), 0 ≤ τ ≤ 2π,
that is q([0, 2π]) = ∂A and q is assumed smooth for our numerical experiments. Following [14, §3.5], one can reformulate [START_REF] Hao | High-order accurate methods for Nyström discretization of integral equations on smooth curves in the plane[END_REF] as an integral equation of the second kind with a weakly singular kernel. There are several methods to discretize such integral equations, see e.g. [START_REF] Hao | High-order accurate methods for Nyström discretization of integral equations on smooth curves in the plane[END_REF] for a review. Here we chose the Kapur-Rokhlin method [START_REF] Kapur | High-order corrected trapezoidal quadrature rules for singular functions[END_REF], which is based on the trapezoidal rule for periodic functions. In this method, the unknowns are the values of ψ at uniformly spaced points of [0, 2π]. To account for the singularity, the entries in a band of the system matrix are weighted so that the quadrature is exact for polynomials of a given order (6th order in our case). We do note that the Kress quadrature [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF] was used for in [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] for computing the scattered fields with k > 0 and is spectrally accurate. Unfortunately, accuracy of the Kress quadrature degrades for complex k. Indeed, the Kress quadrature is obtained by splitting the singular kernel into a singular and non-singular part. The latter requires the evaluation of J 0 (kr), which grows exponentially in Im(k) for fixed r > 0, see e.g. [21, §10.7]. The correction weights for the Kapur-Rohklin method only depend on the type of singularity and order of the method. Thus the Kapur-Rokhlin is better adapted for complex k. Convergence for k complex follows from convergence of the method for the real and imaginary parts, considered individually.

Active exterior cloaking for the heat equation

We now apply the single wavenumber exterior cloaking approach to the time domain heat equation. We recall in section 4.1 other cloaking approaches. We then use the Fourier-Laplace transform to obtain the Helmholtz equation from reasonable heat equation solutions section 4.2. The exterior cloaking approach is applied for different wavenumbers and then put together again in section 4.3 via the inverse Laplace transform. The details of the discrete Fourier transform based algorithm we used for this purpose are in section 4.4.

4.1.

Other cloaking approaches for the heat equation. Cloaking for the heat equation was originally introduced through a change of coordinate system [START_REF] Guenneau | Transformation thermodynamics: cloaking and concentrating heat flux[END_REF], inspired by transformation optics [START_REF] Greenleaf | Anisotropic conductivities that cannot be detected by eit[END_REF][START_REF] Pendry | Controlling electromagnetic fields[END_REF]. However, this approach leads to an extreme anisotropic thermal conductivity, and even a thermal cloak designed through a regularized geometric transform suffers from limited efficiency in the transient regime [START_REF] Guenneau | Transformation thermodynamics: cloaking and concentrating heat flux[END_REF]. A good thermal cloak efficiency requires as many as 10,000 isotropic concentric layers to finely approximate its spatially varying anisotropic conductivity [START_REF] Petiteau | Spectral effectiveness of engineered thermal cloaks in the frequency regime[END_REF]. Thus, fabricated metamaterial cloaks with a limited number of layers suffer from reduced efficiency in the transient regime [START_REF] Schittny | Experiments on transformation thermodynamics: molding the flow of heat[END_REF][START_REF] Narayana | Heat flux manipulation with engineered thermal materials[END_REF][START_REF] Xu | Ultrathin three-dimensional thermal cloak[END_REF][START_REF] Han | Experimental demonstration of a bilayer thermal cloak[END_REF][START_REF] Hu | Illusion thermotics[END_REF]. For other passive cloaking and mimicking approaches see e.g. [START_REF] Alwakil | Inverse heat mimicking of given objects[END_REF][START_REF] Diatta | Non-singular cloaks allow mimesis[END_REF]. Recent advances in thermal cloaking are thus underpinned by inverse homogenization problems that require heavy computational resources. On the other hand, thermoelectric devices have been proposed to pump the heat flow accurately from one side of a thermal cloak to the other side by adjusting the input current, so that the background temperature field can be restored in a stationary regime [START_REF] Nguyen | Active thermal cloak[END_REF][START_REF] Huang | A thermal cloak with thermoelectric devices to manipulate a temperature field within a wide range of conductivity ratios[END_REF]. In our former work [START_REF] Cassier | Active thermal cloaking and mimicking[END_REF], we envisioned using Peltier devices (surrounding the object to cloak) to control transient thermal fields generated by a source. Our approach can be viewed as a generalization of that in [START_REF] Xu | Dipole-assisted thermotics: Experimental demonstration of dipole-driven thermal invisibility[END_REF] that considered a single dipole source placed inside the object to cloak. There should be a trade-off between using a single dipole source and numerous monopole and dipole sources to achieve efficient thermal cloaking in the transient regime, which is what motivated the present work. Our analysis is performed in the frequency regime, where we can extend results of [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF] to the Helmholtz equation with complex wavenumbers. Results are then translated in the time domain through inverse Fourier-Laplace transform.

Remark 4. Since our approach is based on the Laplace transform of the time domain heat equation, it is more convenient to assume a zero initial condition. Indeed a non-zero initial condition would appear as a source term for the Helmholtz equation, which would prevent us from using the interior reproduction formula [START_REF] Friedlander | Introduction to the theory of distributions[END_REF]. However, as noted in [START_REF] Cassier | Active thermal cloaking and mimicking[END_REF], if the initial condition is a steady state solution to the heat equation (i.e. harmonic) we can use the linearity of the heat equation to apply our approach to u(x, t)u(x, 0). Remark 5. As we see next, we obtain solutions to the heat equation that achieve exterior cloaking but they can be large as we get close to the new source locations x j . However, as noted in [START_REF] Cassier | Active thermal cloaking and mimicking[END_REF], it is conceivable to use Peltier devices to physically implement the interior/exterior reproduction formula for the time domain heat equation [START_REF] Costabel | Boundary integral operators for the heat equation[END_REF]. This procedure allows to replace point-like sources by active surfaces that we call "extended cloaking devices", which would keep the temperatures at levels that would be practical to implement.

4.2.

From time domain to frequency domain. We now apply our frequency domain cloaking approach to the heat equation. The temperature u(x, t) (measured in Kelvin) in a homogeneous isotropic medium satisfies the heat equation, [START_REF] Hornung | Diffusion models for fractured media[END_REF] ρc ∂u ∂t = κ∆u + h, for t > 0, where t is the time (s), ρ is the mass density (kg.m -2 ), c is the specific heat (J.K -1 .kg -1 ) and κ is the thermal conductivity (W.K -1 ). Here we assume that ρ, c and κ are positive constants. The source term is h(x, t) (W.m -2 ) and assumed causal, i.e. h(x, t) = 0 for t < 0. For simplicity we assume a zero initial condition and consider (40)

∂u ∂t = σ∆u + h ρc , for t > 0,
where σ = κ/ρc is the thermal diffusivity (m 2 .s -1 ). Assuming further that the source term h(x, t) satisfies the growth condition (7) with L 2 (R 2 ) norm, α ≥ 0 and p = 0, we can see that e -ξt h(x, t)

∈ L 2 ([0, ∞), L 2 (R 2 )) for any ξ > α. Using [18, Corollary 2, p238] it is possible to conclude that (40) admits a unique solution u(x, t) satisfying e -ξt u(x, t) ∈ L 2 ([0, ∞), H 2 (R 2
)) for any ξ > α. This allows to define the Fourier-Laplace transform (6) of all terms in (40) on the half plane C + α , thus obtaining the Helmholtz equation [START_REF] Hu | Illusion thermotics[END_REF] ∆u

(x; ω) + iω σ u(x; ω) = - h(x, ω) ρc ,
where the wavenumber is k = i -iω/σ, using the principal value of the square root and h is the Fourier-Laplace transform of h. We note that Re(k 2 ) = Re(iω/σ) < 0, whenever Im(ω > 0), which guarantees that the Helmholtz equation satisfies a form of the maximum principle for any ω ∈ C + α (since α ≥ 0), for x outside of the support of the source h (see section 3.2 and (8) for the definition of C + α ).

Numerical experiments.

We show in fig. 10 a numerical simulation of active exterior cloaking of a Dirichlet object (a "kite" with constant zero temperature at its boundary) and compare it to the case where there is no cloaking devices. We used σ = 1.5 in [START_REF] Hsu | Numerical inversion of certain laplace transforms by the direct application of fast fourier transform (fft) algorithm[END_REF]. As can be seen from the time snapshot in fig. 10, the isotherm lines without the cloaking devices are significantly different from those of the point source that we used as the incident field u i , this is because of the field "scattered" by the object. For an observer far from the cloaking devices, the isotherms appear consistent with those of a point source, so it is hard for the observer to detect the object from thermal measurements. In our numerical experiments, the region Ω is a disk centered at [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF][START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF] and with radius δ C = 10/6 enclosing the "kite" object. We moved the distribution of monopoles and dipoles to 4 new source locations determined as in fig. 6 with δ D = 5 √ 2/3. We note that three new source locations would have been sufficient, as in [START_REF] Vasquez | Exterior cloaking with active sources in two dimensional acoustics[END_REF]. The fields are calculated on [0, 10] 2 using a 200 × 200 uniform grid. The incident field is generated by a point source at y = [START_REF] Cassier | Code to generate figures in "active exterior cloaking for the 2D Helmholtz equation with complex wavenumbers and appplication to thermal cloaking[END_REF][START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF]. The integral in theorem 2 is approximated with 256 uniformly placed points and the series in [START_REF] Gralak | Macroscopic maxwell's equations and negative index materials[END_REF] uses the truncation M = 22. The boundary of the scatterer is discretized using 512 equally spaced points on the parametric representation of ∂A and the scattered fields (in the frequency domain) are calculated according to the scheme in section 3.3. The frequency domain calculation is performed for 2050 frequencies and the Laplace transform is inverted using a Fast Fourier Transform based method (see section 4.4).

Because the multipolar sources are singular at the x j , the cloaking field diverges as we approach the x j . This could limit the physical implementation (e.g. because the material starts degrading with such high temperatures). Of course, we may use the Green exterior representation formula (see e.g. [START_REF] Costabel | Boundary integral operators for the heat equation[END_REF][START_REF] Cassier | Active thermal cloaking and mimicking[END_REF]) to replace each of the multipolar sources by a monopole and dipole distribution on surface containing the multipolar source. These active surfaces or "extended cloaking device" could be realized in practice using Peltier devices [START_REF] Nguyen | Active thermal cloak[END_REF] and the temperatures do not need be unreasonably large. To illustrate this, we display in fig. 10 black circles centered at the new source positions, outside of which we are guaranteed to have |u (M ) e (x, t)| ≤ u max . Our choice u max is [START_REF] Huang | A thermal cloak with thermoelectric devices to manipulate a temperature field within a wide range of conductivity ratios[END_REF] u max = 100 max

(x,t)∈Ω×[0,T ] |u i (x, t)|.
This choice is not a statement of what is feasible, but simply for illustration purposes. In fig. 10, we have T = 2 and u max ≈ 6.2.

To show that we are achieving exterior cloaking even when replacing the multipolar sources by extended cloaking devices (circular active surfaces), we changed the scale δ C (radius of Ω) of the cloaking configuration in fig. 6 with δ D = √ 2δ C , keeping Ω as a disk with a fixed center (10, 10) and a fixed point source positioned at (10, 1), which generates the incident field (see also section 3.3.1). For all values of δ C , we used the same diffusivity σ = 1.3 in [START_REF] Hsu | Numerical inversion of certain laplace transforms by the direct application of fast fourier transform (fft) algorithm[END_REF], truncation M = 22 and 128 points to discretize a parametric representation of ∂Ω. The temperature fields where evaluated using 130 wavenumbers. The computation was repeated for 40 equally spaced δ C in the interval [START_REF] Martínez | Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements[END_REF][START_REF] Cassier | Code to generate figures in "active exterior cloaking for the 2D Helmholtz equation with complex wavenumbers and appplication to thermal cloaking[END_REF], chosen so that the divergence region R from theorem 2 does not include the source location. For each δ C , the cloak field was evaluated on a 100 × 100 uniform grid of the square [0, 20] 2 and the circles outside of which |u(x, t)| ≤ u max were determined with T = 1 in [START_REF] Huang | A thermal cloak with thermoelectric devices to manipulate a temperature field within a wide range of conductivity ratios[END_REF]. We display in fig. 11 the radius of these circles relative to δ C and as a function of δ C . The dotted line in fig. 11 corresponds to the radius for which the circular active surfaces would touch and match the divergence region R. As we can see from fig. 11, the circular active 6). The dotted line corresponds to the value for which the circles touch. Since the blue data points are below the dotted line, there are gaps between the circular active surfaces, showing that even the "extended cloaking devices" do not completely surround the object (exterior cloaking). sources do not touch, and thus we have exterior cloaking even in this situation. Roughly speaking, according to fig. 11, the "urchins" have a radius that is about 70% of the radius of the gray circles in fig. 6.

4.4.

From the frequency domain to the time domain. For convenience we express the numerical algorithm we use to go from frequency domain to time domain in terms of the Laplace transform of u(x, t) rather than the Fourier-Laplace transform [START_REF] Bretin | On the green function in visco-elastic media obeying a frequency power-law[END_REF]. Under the same growth condition assumption (7) on u(x, t), its Laplace transform is [START_REF] Joseph | Heat waves[END_REF] u(x, s)

= ∞ 0 dt[e -st u(x, t)],
which is well defined for Re(s) > α, where α ≥ 0 is defined in section 4.2. Clearly we have u(x, s) = u(x, is), where the right hand side is the Fourier-Laplace transform of u given in [START_REF] Bretin | On the green function in visco-elastic media obeying a frequency power-law[END_REF]. The inverse Laplace transform is then given by ( 44)

u(x, t) = 1 2πi c+i∞ c-i∞
ds[e st u(x, s)],

for any c > α. 6 We follow the numerical method in [START_REF] Hsu | Numerical inversion of certain laplace transforms by the direct application of fast fourier transform (fft) algorithm[END_REF] for computing the inverse Laplace transform by approximating it with a discrete Fourier transform (or DFT, which can be evaluated efficiently with the Fast Fourier Transform or FFT see e.g. [START_REF] Frigo | The design and implementation of fftw3[END_REF] for definition). For a review of numerical inverse Laplace transform methods see e.g. [START_REF] Cohen | Numerical methods for Laplace transform inversion[END_REF]. The idea is that to a uniform grid of the time interval [0, T ] with N points, i.e. t p = p∆t ∈ [0, T ], p = 0, . . . , N -1, ∆t = T /N , we associate the discretization of a dual variable w given by w q = q∆w, q = 0, . . . , N -1, ∆w = 2π/T . The discretizations are chosen such that t p w q = 2πpq/N , which is the negative of the phase of the complex exponential in the DFT of length N . Using the change of variables s = c + iw in eq. ( 44) and approximating with a Riemann sum on a finite interval yields the following (we assume u(x, t) is real)

u(x, t) = 1 2π ∞ -∞ dw [exp(ct) exp(iwt)u(x, c + iw)] ≈ 1 2π N -1 q=-(N -1)
∆w [exp(ct) exp(iw q t)u(x, c

+ iw q )] = exp(ct) T N -1 q=-(N -1) [exp (iw q t) u (x, c + iw q )] = 2 exp(ct) T Re N -1 q=0 [exp (-iw q t) u (x, c -iw q )] - u(x, c) 2 ,
where the last equality follows by the the symmetry for Laplace transforms of real functions (u(x, s) = u(x, s)). By evaluating at the t p , our approximation can be written in terms of the real part of a DFT:

u(x, t p ) = 2 exp(ct p ) T Re N -1 q=0
[exp (-it p w q ) u (x, ciw q )] -u(x, c) 2 = 2 exp(ct p ) T Re fft(u (x, ciw q ) , q = 0, . . . , N -1) -u(x, c) 2 .

Here fft(v) represents the DFT of a vector v of length N , as defined in [START_REF] Frigo | The design and implementation of fftw3[END_REF]. Thus we end up evaluating the frequencies ω q = w q + ic, q = 0, . . . , N -1. As noted in [START_REF] Hsu | Numerical inversion of certain laplace transforms by the direct application of fast fourier transform (fft) algorithm[END_REF], the convergence only holds for [0, T /2), so in practice we use N = 2N + 2 and only use the first N + 1 time steps to give convergence on the interval [0, T ] .

Here the additional two frequencies mean the terminal time step is T as opposed to T -∆T. Since the heat equation solutions we consider decay, we take α = 0 and set c = α -∆w 2π ln(10 -6 ) > 0. We note that there are methods to speed up convergence of this class of numerical inverse Laplace transform, see e.g. [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF]. 6 We note that the assumptions of section 4.2 (i. e. the growth control of the source term, the zero initial condition and the causality of u) imply that for any ξ > α ≥ 0 : e -ξt u(x, t) ∈ L 2 ([0, ∞), H 2 (R 2 )) and (using the heat equation) that e -ξt ∂tu ∈ L 2 ([0, ∞), L 2 (R 2 )). Thus, one has that u(•, s) = s -1 L(∂tu)(•, s) is analytic with respect to s for Re(s) > α, where L stands for the Laplace transform. Since L(∂tu)(•, s) → 0 for Re(s) ≥ c > α and |s| → ∞, we get u(•, s) = o(|s| -1 ). We need a little more of decay to use the formula [START_REF] Kadic | Metamaterials beyond electromagnetism[END_REF] for the inverse Laplace transform. It is enough to assume that there exists ε > 0 such that u(•, s) ≤ C|s| -(1+ε) for Re(s) ≥ c > α. Then, [START_REF] Kadic | Metamaterials beyond electromagnetism[END_REF] is well defined as a Bochner integral valued in L 2 (R 2 ) (see e.g. the proof of [3, Theorem 2.5.1] or [START_REF] Sayas | Retarded potentials and time domain boundary integral equations[END_REF]) and thus also pointwise for a. e. x ∈ R 2 and all t ≥ 0.

Summary and perspectives

In our earlier work on active exterior cloaking for the parabolic heat equation [START_REF] Cassier | Active thermal cloaking and mimicking[END_REF], we noted in the concluding remarks that our approach could also be tied to the active exterior cloaking strategies for the Helmholtz equation by going to Fourier or Laplace domain in time. Here we have showed that it is possible to cloak objects from thermal measurements by using active heat sources, starting from a zero temperature condition. We believe that our work opens up a new path for active cloaking for a variety of physical situations, in addition to the class of differential equations of the form (5). Indeed, complex wavenumbers make it possible to model pseudodifferential operators in time such as fractional time derivatives [START_REF] Podlubnv | Fractional differential equations academic press[END_REF] and integro-differential equations. This opens new avenues in active exterior cloaking. To give an example of this flexibility, consider a (non-dimensionalized) heat equation with memory that arises when considering homogenized diffusion models in fractured media [

(x, ω) -(-iω + p(ω))u(x, ω) = -h(x, ω).
The present work allows us to also study active cloaking in the context of diffusive photon density waves governed by [START_REF] Kresin | Criteria for validity of the maximum modulus principle for solutions of linear strongly elliptic second order systems[END_REF] ∂u ∂t + µu = σ∆u + h, for t > 0, where µ > 0 is an absorption coefficient, σ is a conductivity and h represents the photon current density (photon flow per unit surface and per unit time). Making use of the Fourier-Laplace transform (6) and assuming zero initial conditions, [START_REF] Kresin | Criteria for validity of the maximum modulus principle for solutions of linear strongly elliptic second order systems[END_REF] takes the form of the Helmholtz equation in the context of diffusion wave scattering

(48) ∆u(x; ω) + iω -µ σ u(x; ω) = - h(x, ω) σ ,
where we note that (48) reduces to [START_REF] Hu | Illusion thermotics[END_REF], when ω µ. Scattering cancellation of such diffusive waves has been addressed in [START_REF] Schittny | Invisibility cloaking in a diffusive light scattering medium[END_REF][START_REF] Farhat | Cloaking through cancellation of diffusive wave scattering[END_REF].

Moreover, advection-diffusion problems play a prominent role notably in diffusion and mixing of fluid flow modelled by [START_REF] Constantin | Diffusion and mixing in fluid flow[END_REF] (49)

ρc ∂u ∂t = -v • ∇u + κ∆u + h, for t > 0,
where v is a constant velocity, ρ is a mass density, κ a conductivity and c the heat capacity. The same equation is known as the Fokker-Planck equation and is central to models for transport of salt, heat, buoys, and markers in geophysical flows [START_REF] Murphy | Spectral analysis and computation for homogenization of advection diffusion processes in steady flows[END_REF][START_REF] Zelati | Homogenization and hypocoercivity for fokker-planck equations driven by weakly compressible shear flows[END_REF]. It turns out that one can recast (49) using the exponential variable transform [START_REF] Li | Boundary element solution of heat convection-diffusion problems[END_REF] (50)

u(x, t) = exp[ v 2 κ • x]w(x, t) and h(x, t) = exp[ v 2 κ • x] g(x, t)
together with the Fourier-Laplace transform [START_REF] Bretin | On the green function in visco-elastic media obeying a frequency power-law[END_REF], into the Helmholtz equation (assuming zero initial conditions)

(51) ∆w(x; ω) -τ 2 w(x; ω) = - g(x, ω) κ ,
where τ 2 = (| v | /2 κ) 2iω/σ, and σ = κ/(ρc).

Finally, the exterior cloaking theory which we developed may allow us to also achieve exterior cloaking in the context of Maxwell-Cattaneo heat waves governed by [START_REF] Joseph | Heat waves[END_REF] (52)

τ ∂ 2 u ∂t 2 + ∂u ∂t = κ∆u + τ σ∆ ∂u ∂t + h, for t > 0,
where κ is the thermal conductivity, σ accounts for diffusive phenomena, τ is the thermal relaxation time (that corresponds to the time it takes for a medium to reduce its temperature to half). We assume that κ, τ , σ are positive constants. Making use of the Fourier-Laplace transform (6) with zero initial conditions, (52) takes the form of the Helmholtz equation in the context of diffusion wave scattering

(53) ∆u(x; ω) + ω(ωτ + i) κ -iωτ σ u(x; ω) = - h(x, ω) κ -iωτ σ ,
and thus one can define k 2 = ω(ωτ + i)/(κiωτ σ) that unlike for the Fourier heat equation ( 41) can lead to propagating features (when Re(k 2 ) > 0). Scattering cancellation of such diffusive waves has been addressed in [START_REF] Farhat | Scattering cancellationbased cloaking for the maxwell-cattaneo heat waves[END_REF].

Another situation where complex wavenumbers may play a prominent role is for in-plane pressure and shear elastodynamic waves propagating in passive, dissipative, active or even viscoelastic media (the case of anti-plane shear waves would be covered by the 2D scalar Helmholtz equation with complex wavenumber). The latter, viscoelastic media, would require using the Helmholtz decomposition proposed in [START_REF] Bretin | On the green function in visco-elastic media obeying a frequency power-law[END_REF] u = ∇Φ + ∇ × Ψ, with Ψ a divergence free vector field, in the vector Navier equation [START_REF] Miller | On perfect cloaking[END_REF] ( λ + 2μ)∇∇ • u + μ∇ 2 u + ρω 2 u = 0 where λ = λ + η p M, μ = µ + η s M, M being a convolution operator with certain power law (named after Szabo and Wu [START_REF] Szabo | A model for longitudinal and shear wave propagation in viscoelastic media[END_REF]), λ, µ are the usual Lamé parameters, ρ is the density and η s , η p 1. Our approach would then be applied to Helmholtz equations with the complex wavenumbers for shear (s) and pressure (p) waves ( 55)

k 2 m (ω) = ω 2 1 - ν m c 2 m M(ω) , m = s, p ,
with c p = (λ + 2µ)/ρ, c s = µ/ρ, the pressure and shear wave velocities, respectively, ν p = (η p + 2η s )/ρ, ν s = η s /ρ and where the multiplication operator M is the Fourier transform of the convolution operator M. The Hodge decomposition can also be applied to split the elastodynamic wave equation with isotropic viscoelasticity into two acoustic wave equations in time (P and S waves) that can be transformed into the Helmholtz equation, see e.g. [START_REF] Martínez | Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements[END_REF]. There may also be ways of adapting our approach to other differential operators in space. One example would be to consider constant anisotropic media (e.g. coming from a homogenization approach). We believe that our strategy could be also applied to [START_REF] Brančík | Matlab oriented matrix laplace transforms inversion for distributed systems simulation[END_REF] wherein the Laplacian is replaced by a bi-Laplacian in the right-hand side. This problem arises when modeling flexural waves in thin elastic plates, which are governed by the bi-harmonic (Kirchhoff-Love) equation. Active cloaking in this context has been considered in [START_REF] O'neill | Active cloaking of inclusions for flexural waves in thin elastic plates[END_REF]. We are also considering applying our theory to active cloaking for flexural gravity waves in floating thin elastic plates that would involve a tri-Laplacian in the right-hand side of (5) [START_REF] Farhat | Scattering theory and cancellation of gravity-flexural waves of floating plates[END_REF].

We conjecture that our approach can be adapted to the 3D Helmholtz equation for complex wavenumbers, which would allow us to address cloaking in general dispersive media (including media with losses or gain). Finally we note that open questions related to gain media (Im(k) < 0) remain. In particular what is a sensible functional space setting for the exterior Green representation formula and for the scattering problem in gain media. believe that the approach of exterior cloaking which we have developed in this article allows us to cover a broad range of problems for active cloaking of diffusion and wave phenomena.

Appendix A. Proof of Lemma 1 Proof.

Step 1: We prove first [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF]. The entire function J n is defined [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF][Chapter, formula 1 p. 15] via the following power series: (56) We point out that f (z) = 4(exp(z 2 /4) -1)/z 2 for z = 0 and f (0) = 1. Thus, ( 18) is an immediate consequence of ( 56) and (57).

J n (z) := ∞ k=0 (-1) k k! (n + k)! z 2 n+2k = 1 n! z 2 n + ∞ k=1 (-1) k k! (n + k)! z 2 n+2k , ∀z ∈ C.
Step 2: The asymptotic formula [START_REF] Diatta | Non-singular cloaks allow mimesis[END_REF] for J n follows immediately from the asymptotic expansion [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF] and the recurrence formula (see [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF] formula 2 page 45) J n (z) = 1 2 (J n-1 (z))-J n+1 (z)). An other way to obtain the formula ( 19) is to derive the power series that defines J n and applies the same method as for formula [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF].

Step 3: We now prove [START_REF] Disalvo | Thermoelectric cooling and power generation[END_REF]. By definition of the Hankel function H (see [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF], formula (1) page 73) , one has for all z ∈ C \ (-∞, 0]: H (1) n (z) -i (n -1)! π , where ln is the principal value of the logarithm function with a branch cut on (-∞, 0] and ψ := Γ /Γ with Γ the well-known Gamma function. We estimate now the two terms appearing in [START_REF] Norris | Source amplitudes for active exterior cloaking[END_REF]. For the first one, one obtains i π

2 z n n-1 k=1 (n -k -1)! k! z 2 2k ≤ (n -2)! π 2 |z| n-2 n-1 k=1 1 k! |z| 2 2k-2 . ≤ C K2 π (n -2)! 2 |z| n-2 , ( 60
)
with the constant C K2 defined by replacing the compact K 1 by K 2 in [START_REF] Nguyen | Active thermal cloak[END_REF]. As the function ψ evaluated on integers is given by (see [21, §5.4.14] Data Access. We provide the following supplementary materials. (i) A movie animating fig. 10. (ii) The Matlab code to reproduce the figures figs. 4, 5 and 7 to 11 is available in the repository [START_REF] Cassier | Code to generate figures in "active exterior cloaking for the 2D Helmholtz equation with complex wavenumbers and appplication to thermal cloaking[END_REF].

Author Contributions. SG initiated the project. MC proved Lemma 1. MC and TD proved the error estimates for the Graf addition formula. MC, TD and FGV proved the error estimates for the truncated cloaking field. TD and FGV performed the numerical experiments. SG found that the approach was valid even for media with gain. All authors contributed to the writing. All authors contributed other idea applications. All authors gave final approval for publication and agree to be held accountable for the work performed therein.

Funding. TD and FGV were supported by the National Science Foundation Grant DMS-2008610.

( 3 )

 3 D xj ,y := {x ∈ R 2 | |xx j | ≤ |yx j |}, where | • | denotes the Euclidean norm. By linearity, we should also be able to move a distribution of monopoles or dipoles on a compact portion ∂Ω j of the boundary ∂Ω to a new location x j , obtaining the same cloak field u c , provided we are outside of the closed disk (4) R j = {x ∈ R 2 | |xx j | ≤ max y∈∂Ωj |yx j |}.

  e. the annulus for which the ratio in the geometric series ansatz belongs to [1/2, 1].

Figure 3 .

 3 Figure 3. Wavenumber ranges used in the numerical experiments and their visualization in the complex plane.

Figure 4 .

 4 Figure 4. The monopole (a) and dipole (b) errors (logarithmic scale), R j,M and R j,M in (15) respectively, and the bounds from theorem 1 for a single point in space which was relatively close to the source. The color represents the set for k as given in fig.3. For the dipoles we cannot differentiate the different bounds as they all lie on top of each other. We observe that the dipole errors and bounds are larger. This is in line with theorem 1: the dipole error bound decays slower than the monopole error bound.

Figure 5 .

 5 Figure 5. The monopole (a) and dipole (b) errors (logarithmic scale), R j,M and R j,M in (15) respectively, and bounds for a range of points in space at a real k = 1, pure imaginary k = i, complex (dissipative) k = 1/ √ 3 + i √ 2/ √ 3, and complex (gain) k = (99i √ 199)/100 wavenumbers with the colors represented in fig. 3. The plots for different wavenumbers cannot be differentiated.

Remark 1 .

 1 The representation formula (28) is valid for example when Ω has Lipschitz boundary ∂Ω. To see this, we assume Ω ⊂ O where O is an open set and the incident field u i ∈ H 1 loc (O) solves (in the distributional sense) the homogeneous Helmholtz equation

Remark 3 . 3 . 1 .

 331 Controlling fields outside of a bounded open set Ω can be useful for the mimicking problem (making a scatterer inside Ω look like another one) or for cloaking a source inside Ω. This requires an exterior version of the Green representation formula[START_REF] Friedlander | Introduction to the theory of distributions[END_REF], which is valid for Im(k) ≥ 0, when the field to reproduce u i is a solution to the homogeneous Helmholtz equation outside of Ω and satisfies the Sommerfeld radiation condition (see e.g.[START_REF] Colton | Integral equation methods in scattering theory[END_REF] Theorem 3.3]). We are not aware of the validity of this result for gain media (Im(k) < 0). Therefore we anticipate that theorem 2 can be adapted to control fields outside of R for Im(k) ≥ 0 and k / ∈ (-∞, 0]. Proof of theorem 2.

3 . 2 .

 32 bound follows by letting C = N dev C max j |∂Ω j |. In addition, we have for x / ∈ R and k ∈ C \ (∞, 0] that |u c (x)u e (x)| = lim M →∞ |u c (x)u (M ) e (x)| = 0 since a < 1. Hence we have u c = u e for x / ∈ R. The fields do not agree for x ∈ R, because for at least one j ∈ {1, . . . , N dev }, the series in (29) diverges. Stability through the maximum principle. Homogeneous Helmholtz equation solutions satisfy a strong maximum principle if Re(k 2 ) < 0 or equivalently (35) |Im(k)| > |Re(k)|. Although we use this result for constant isotropic media, it has been proved in the very general context of the Helmholtz equation with anisotropic heterogeneous media [47, corollary 2.1]. Another proof in the case of isotropic heterogeneous media appears in [48, theorem 6]. We now state the strong maximum principle. By interior regularity (see Remark 1), a homogeneous solution u ∈ H 1 loc (O) to the Helmholtz equation in an open set O with a wavenumber k satisfying (35), is C ∞ (B) on any open bounded subset B with Lipschitz boundary ∂B satisfying B ⊂ O. Thus, one can apply the strong maximum principle on the set B to get on one hand that max x∈B |u(x)| = max x∈∂B |u(x)|,

  are C ∞ (B) solutions to the homogeneous Helmholtz equation on B, where B is a bounded open set such that B ⊂ R 2 \ R (see theorem 2). Therefore we can conclude from the maximum principle that the truncation error of exterior cloaking (|u e (x)u (M ) e (x)|) reaches its maximum over B only on the boundary ∂B.

3. 3 .

 3 Numerical experiments. We explain how we evaluate the truncated cloak field u (M ) e in section 3.3.1. Then the truncation errors are predicted in section 3.3.2 using the error bounds in theorem 2. Finally we explain in section 3.3.3 how we calculate scattered fields when Im(k) ≥ 0.

  Figure 6. (a) The configuration of the exterior cloak used for four exterior sources (x j ) to maximize the region where an object can be dissimulated (the blue striped region), we take δ D = √ 2δ C . (b) Cloaking region with smallest inscribed circle with radius r Ci and largest circumscribed circle with radius r C O .

Real part of field Real part of reproduction k = 10 k = i/ 2 k = 10 + i/ 2 k = 10 -i/ 2 Figure 7 .

 1021021027 Figure 7. Field reproductions (right) and the original field (left) at different wavenumbers on the square [0, 1] 2 with a point source located at (2, 5). The color scale was kept the same for each k and was chosen to highlight the different behaviors of the point sources for different k.

Figure 8 .

 8 Figure 8. Predicted and actual maximum error (logarithmic scale) on circles of radius slightly larger than r c O and smaller than r ci , which are outside of the divergence region of the cloak field. The abscissa corresponds to the parameter θ for the segment of wavenumbers K 3 , as defined in fig. 3.

Figure 9 .

 9 Figure 9. We display the cloaking field truncation error log 10 |u M e -u e | (outside of the "extended cloaking devices" in white) corresponding to (a) k = i/2 and (b) k = 10 + i/2. We note that (a) (resp. (b)) corresponds to the second (resp. third) row in fig. 7. By applying the maximum principle on the disk B (dashed curve), we see that the maximum error is attained on ∂B in (a) but not in (b). The difference is that the wavenumber in (a) satisfies |Im(k)| > |Re(k)| so a version of the maximum principle applies, see also fig. 2.

Figure 10 .

 10 Figure 10. (a) Temperature distribution at time t = 4 with σ = 1.5 (see (40)) resulting from a point source at location (8, 5) in the presence of a "kite" object with homogeneous Dirichlet boundary condition. (b) The same object and source as in (a), but with the cloaking devices activated. The temperatures outside of the black circles are bounded by u max ≈ 6.2 (see (42)) for t ∈ [0, 4]. The computational domain was [0, 10] 2 . The linear color scale spans temperatures in [-0.0133, 0.0133] and range from blue (negative) to red (positive), with zero been represented in white. [See also movie in supplementary material]

Figure 11 .

 11 Figure 11. The blue curve corresponds to circular active surface radii (relative to the scaling parameter δ C ) to achieve exterior cloaking of a point source located at (10, 1) for 40 different values of the scaling parameter δ C and for Ω centered at (10, 10) (see fig.6). The dotted line corresponds to the value for which the circles touch. Since the blue data points are below the dotted line, there are gaps between the circular active surfaces, showing that even the "extended cloaking devices" do not completely surround the object (exterior cloaking).

As 1 /|z| 2 n+2 2 ,

 122 (n + k)! ≤ 1/(n + 1)! for k ≥ 1, one gets[START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF] where the positive constant C K1 is defined by (58) C K1 = max z∈K1 f (|z|) with f entire defined by f (z) := ∀z ∈ C.

=n 2 )

 2 J n (z) + iY n (z) -i (n -1)! π 2 zUsing[START_REF] Narayana | Heat flux manipulation with engineered thermal materials[END_REF] and the series representation of Y n on C \ (-∞, 0] (see[START_REF] Watson | A treatise on the theory of Bessel functions[END_REF], formula (3) page 62 and (2) page 64 or formula (10.8.1) of[21]) on the previous expression leads to:(59)H (1) n (z) -i (n -1)! π ψ(k + 1)ψ(k + n + 1) + 1 (-1) k k!(n + k)! z 2 2k

z 2 ||I 2

 22 formula 5.4.14) ψ(m + 1) := m p=1 1 p + γ for m ≥ 1 with γ =: Γ (1) = ψ(1) the Euler constant, the second term of (59) can be bounded by |I 2 (n, z)| ≤ C = 2γ+1 (we point out that we use the inequality ψ(k+1) ≤ ψ(k+2) for k ≥ 0 to avoid the particularity of the case k = 0). Then, using the following inequality m p=1 1/p ≤ ln(m) + 1 for m ≥ 1 (obtained by comparison of the harmonic series with the integral) and the fact that ln(k + 1) ≤ ln(n + k), one gets that|I 2 (n, z)| ≤ + C + 2 ln(n + k) + 2 1 k!(n + k)! |z| 2 2k .Notice that ln(n + k)(n + k)! ≤ ln(n + k) (n + k)(n + k -(n, z)| ≤ CK2 (n-2)!2 |z| n-2 with CK2 = 1 π max K2 g(|z|)(2| ln z 2 |+C +4)h(|z|).

  a bounded open set with Lipschitz boundary ∂Ω. Assume u i is a H 1 loc (O) solution to the homogeneous Helmholtz equation, where O is an open set containing Ω. Define the

  = 0 with ηRe(k) ≥ 0 is a coupling parameter. This choice guarantees invertibility (see e.g.[13, §3.6]). Here it is necessary to define the corresponding boundary layer operators for x ∈ ∂A by for the single and double layer potential respectively. We note that the operators can be taken as bounded operators S, K : L 2 (∂A) → L 2 (∂A) (see e.g. [57, Theorem 4.4.1] for smooth ∂A or[52, chapter 6] for C 2 or even Lipschitz ∂A). We also need the following jump relations, letting z ∈ R 2 \ ∂A

	(36)	1/2 (∂A) satisfying (Kϕ)(x) := 2 lim z→x ∂A dS(y) [G(z -y; k)ϕ(y)] = ∂A dS(y) [G(x -y; k)ϕ(y)] , and ∂A dS(y) ∂G ∂ν(y) (x -y; k)ϕ(y) , [Sϕ](x) 2 , and u (Sϕ)(x) := 2 lim z→x + ∂A ∂G dS(y) ∂ν(y)

s (x; k) = ∂A dS(y) ∂G ∂ν(y) (xy; k) -iηG(xy; k) ψ(y)

where η

We recall that for p ≥ 1, L p loc (X) (resp. H 1 loc ) is the set of functions that are L p (resp. H 1 ) on any bounded open subset of X, with closure inside X. See e.g.[START_REF] Evans | Partial differential equations[END_REF][START_REF] Lions | Non-homogeneous boundary value problems and applications[END_REF].