Active thermal cloaking and mimicking

We present an active cloaking method for the parabolic heat (and mass or light diffusion) equation that can hide both objects and sources. By active, we mean that it relies on designing monopole and dipole heat source distributions on the boundary of the region to be cloaked. The same technique can be used to make a source or an object look like a different one to an observer outside the cloaked region, from the perspective of thermal measurements. Our results assume a homogeneous isotropic bulk medium and require knowledge of the source to cloak or mimic, but are in most cases independent of the object to cloak.


Introduction
Certain solutions to the heat equation can be reproduced inside or outside a closed surface by a source distribution on the surface determined by Green's identities.In particular, given a solution to the heat (or mass or light diffusion) equation in a homogeneous medium and with no sources inside of a domain, it is possible to reproduce it inside the domain with a distribution of sources on the surface of the domain, while also giving a zero solution outside.We call this the interior reproduction problem, see fig.1(a).Similarly, the exterior reproduction problem is to reproduce a solution to the heat equation in a homogeneous medium with no sources outside of a domain, while keeping a zero solution inside the domain (see fig. 1(b)).As we shall see, a growth condition for the heat equation solution is needed to guarantee that the exterior reproduction problem can be solved.
This growth condition plays the same role as a radiation boundary condition for the Helmholtz equation (section 2).By combining solutions to interior/exterior reproduction problems we can achieve cloaking or mimicking for the heat equation in the following scenarios.
Interior cloaking of a source: Given a localized heat source, find an active surface or cloak surrounding the source so that the source cannot be detected by thermal measurements outside the cloak (section 3.1).
Interior cloaking of an object: Given a passive object (e.g. an inclusion), find an active surface or cloak surrounding the object so that the temperature distribution outside the cloak is indistinguishable from having a region of homogeneous medium instead of the cloak and the object (section 3.2).
Source mimicking problem: Given a localized source, find an active surface or cloak surrounding it so that the source appears as a different source for an observer outside the cloak (section 4.1).
(a) (b) Figure 1.In (a) we illustrate the "interior reproduction problem" which consists of reproducing a solution to the heat equation in the interior of a bounded region Ω (in yellow), while enforcing a zero solution outside of Ω (in white) by placing heat sources on the boundary ∂Ω.In (b) we illustrate the "exterior reproduction problem" in a similar way.
Figure 2.An illustrative example of an arrangement of Peltier elements that could be used to cloak objects (e.g. a kite) inside a two-dimensional region Ω, illustrated here by a disk within a heat conducting plate.Each Peltier device is represented by two adjacent red and blue boxes, where the heat flux it can create is oriented in the direction normal to their interface.We have represented both Peltier elements that transport heat within the plate (across ∂Ω) and also between the exterior and the plate.
Object mimicking problem: Given a passive object inclusion, find an active surface or cloak surrounding it so that the object appears as a different object for an observer outside the cloak (section 4.2).
Active sources for the "interior cloaking of an object" problem have already been proposed and demonstrated experimentally for the steady state heat equation in [1].The idea is to use Peltier elements to dissimulate a hole in a conductive plate under a steady state temperature distribution.A Peltier element is a thermoelectric heat pump that moves heat by putting an electric current across a metal/metal junction (see e.g.[2]).Our approach deals with time varying solutions to the heat equation.It could be implemented with a setup similar to [1] by using Peltier devices that are either used to transport heat within the plate or act as heat sources/sinks on the plate by transporting heat between the environment and the plate, as illustrated in fig. 2.An alternative route using a single active dipole source placed inside the object to cloak in a constant gradient steady state regime is proposed in [3].
Although all the results are presented in the context of the heat equation, they also apply to the diffusion equation, which can be used to model e.g.diffusion of a species in a porous medium.In this case the active surface would consist of pumps that can transport the species either across the medium or between the medium and the environment.For the diffusion equation case, we can either hide or imitate a source or an inclusion with different diffusivity properties.
Controlling the heat flux may find applications in enhancing the efficiency of thermal devices in solar thermal collectors, protecting electronic circuits from heating, or the design of thermal analogs of electronic transistors, rectifiers, and diodes [4].Moreover, all our results could be easily adapted to control of mass diffusion with potential applications ranging from biology with the delay of the drug release for therapeutic applications [5,6] to civil engineering with the control of corrosion of steel in reinforced concrete structures [7].We further note that in many media, such as clouds, fog, milk, frosted glass, or media containing many randomly distributed scatterers, light is not described by the macroscopic Maxwell equations, but rather by the Fick's diffusion equation as photons of visible light perform a random walk.Cloaking for diffusive light was experimentally achieved in [8] using the transformed Fick's equation [5] and this suggests potential applications of our work in control of diffusive light as well.
Because we are using source distributions, we can achieve cloaking of inclusions and sources, regardless of how complicated they are and on arbitrarily large time intervals.One drawback of our method is that the source distribution completely surrounds the object or source that we want to cloak or mimic.Another drawback is that we assume perfect knowledge of the fields to reproduce on a surface, however we expect this can be relaxed, as we demonstrate numerically in section 2.3.Apart from the active cloaking strategies for the steady state heat equation in [1,3], there are passive cloaking methods for the heat equation that use carefully crafted materials to hide objects [9,10,11,12,5].Such materials, whose effective conductivity mimics that in the heat equation after a suitable change of variables has been made, are quite bulky.In [12], some proof of concept of passive thermal cloaking was achieved with a metamaterial cloak consisting of 10 concentric layers mixing copper and polydimethylsiloxane in a copper plate.However, it has been numerically shown using homogenization in [13] that one would require over 10,000 concentric layers with an isotropic homogeneous conductivity to accurately mimic the required anisotropic heterogeneous conductivity within a thermal cloak in order to achieve some markedly improved cloaking performance in comparison with [12].
The idea of using active sources based on the Green identities to cloak objects was first proposed for waves by Miller [14].The way of finding the sources for active cloaking is similar to that in active sound control for e.g.noise suppression [15,16].One problem with this approach is that the sources completely surround the cloak.However only a few sources are needed to cloak as was shown in [17,18,19,20,21,22] for the Laplace and Helmholtz equations.The approach can be extended to elastic waves [23] and flexural (plate) waves [24,25].The active cloaking approach can be applied to the steady state diffusion equation with the role of sources and sinks over a thick coating played by chemical reactions [26].Active cloaking has been demonstrated experimentally for the Laplace equation [27,28], electromagnetics [29] and the steady state heat equation [1,3].Finally the illusion/mimicking problem was first proposed using metamaterials via a transformation optics approach [30] and then with active exterior sources [31].
We start in section 2 by recalling results on representing solutions to the heat equation by surface integrals.This section includes a growth condition on heat equation solutions that is sufficient to ensure that the exterior reproduction problem is solvable and numerical experiments illustrating both the interior and exterior reproduction problems in two-dimensions.We also underline properties related to the maximum principle of the heat equation that we use on one hand to point out some stability of the two reproduction problems and on the other hand to interpret the numerical error in our simulations.Then in section 3 we explain how to cloak a source (section 3.1) or an object (section 3.2).The mimicking problem is presented in section 4 for both mimicking objects (section 4.2) and sources (section 4.1).The numerical method we use to illustrate our approach in two-dimensions is explained in section 5. Finally, our results are summarized in section 6.

Integral representation of heat equation solutions
We start by recalling results on boundary integral representation of solutions to the heat equation.Concretely we show in sections section 2.1 and section 2.2 how to use a distribution of monopole and dipole heat sources on a closed surface, an "active surface", to reproduce a large class of solutions to the heat equation inside and outside the surface.We include a two-dimensional numerical study (section 2.3) of how the reproduction error is affected by the time step, the boundary discretization and errors in the boundary density we use to represent the fields.
The temperature u(x, t) of a homogeneous isotropic body satisfies the heat equation (1) where u is in Kelvin, x is the position in meters, t is the time in seconds, κ is the thermal conductivity (W m −1 K −1 ), c is the specific heat (JK −1 kg −1 ), ρ is the mass density (kg m −3 ) and h(x, t) is a source term (W m −3 ).To simplify the exposition we consider instead (2) where k = κ/ρc is the thermal diffusivity (m 2 s −1 ) and h = h/ρc is the source term (K s −1 ).In dimension d the Green function or heat kernel for ( 2) is where We point out that outside the origin (x, t) = (0, 0), K(x, t) is a smooth function even on the line t = 0.This can be shown directly or by using the hypoellipticity property of the heat operator.Namely, as the heat kernel solves the homogeneous heat equation in the distributional sense on any open set that does not contain the origin, by hypoellipticity (see [32] theorem 1.1 page 192) K(x, t) can be extended as a ) is to reproduce a solution of the homogeneous heat equation in the space-time cylinder Ω × (0, ∞) by placing appropriate source densities on its boundary ∂Ω × [0, ∞) and possibly on Ω × {0} (initial condition).The sources are chosen such that the fields vanish in (R d − Ω) × (0, ∞) (where Ω = Ω ∪ ∂Ω denotes the closure of Ω).If the initial condition is harmonic, it is sufficient to have sources on ∂Ω × [0, ∞) only.For the exterior reproduction problem (section 2.2) we seek to reproduce a solution of the homogeneous heat equation in (R d − Ω) × (0, ∞) using sources on ∂Ω × [0, ∞) and possibly on (R d − Ω) × {0}.The fields are required to vanish in Ω × (0, ∞).

2.1.
Reproducing fields in the interior of a bounded region.The goal here is to reproduce solutions u to (2) in Ω by controlling sources on ∂Ω while leaving the exterior unperturbed.More precisely, for some temperature field u(x, t), we wish to generate u Ω (x, t) such that for t > 0: Notice that u Ω (x, t) is not defined for x ∈ ∂Ω, as usual in boundary integral equations.For initial condition u(x, 0) = f (x), x ∈ Ω and a source term h(x, t) supported in R 2 − Ω, this can be achieved via the Green identities (see e.g.[33,34,35]) here n(y) is the outward pointing unit length normal vector at y ∈ ∂Ω.The Green identities guarantee that u Ω (x, t) = 0 for x / ∈ Ω, as desired.The first term in the integral ( 5) is the single layer potential (a collection of monopole heat sources) and the second is the double layer potential (a collection of dipole heat sources).For zero initial conditions, the solution u is completely represented within Ω by the single and double layer potentials, with densities given by the field to be reproduced and its normal derivative on ∂Ω.We point out that the representation formula (5) holds for instance if u ∈ C 2 (Ω × [0, +∞)).A less restrictive condition, for zero initial condition, is to assume some Sobolev regularity for u, see e.g.[34, theorem 2.20].In this less regular setting, the first integral in (5) becomes a duality pairing between boundary Sobolev spaces.
Remark 2.1 (Causality and instantaneous control).The boundary integral representation (5) is causal in the sense that to reproduce u Ω (x, t), t > 0, we only need information about u in the past, i.e. for times before the present time t.Moreover, the time convolution (5) suggests that at the present time t, we only require control of heat sources localized in time to the present time t.Indeed, the integral over ∂Ω in (5) is a collection of monopole and dipole sources localized in time to s and that depends only on knowing u and ∂u/∂n at time s.Moreover the contribution of past s, i.e. with s < t, amounts to the memory effect of the bulk.Thus for experimental purposes, the boundary integral representation could be approximated by e.g.Peltier devices.
A numerical example is given in fig. 3.Here the field u is generated by a point source at x = (0.25, 0.25) and t = 0.For the heat equation we took k = 0.3 and the domain is Ω = B(x 0 , r 0 ), the open ball with center x 0 = (0.5, 0.5) and radius r 0 = 0.25.It should be noted that in all of the numerical examples the thermal diffusivity k is chosen for the convenience of computations, and may be The field is reconstructed inside a disk of radius 1/4 centered at (1/2, 1/2) by using only heat sources on the corresponding circle.
In (c) we show the log 10 of the reconstruction error (the absolute value of the difference between the exact u Ω and its numerical approximation).We generated a plot similar to (c) by taking the maximum over the time interval [0.2, 0.3] at each grid point (the plot being very similar to (c), we include the code to generate it as supplementary material).This indicates that the maximum error is attained near ∂Ω in space (and in time at t = 0.2 on the time interval [0.2, 0.3]), conforming to the maximum principle (applied with initial time t = 0.2) for the interior problem (remark 2.2) and the exterior one (remark 2.8).We can use the maximum principle on the computed fields because the numerical method we use generates solutions to the heat equation (see section 5 for more details).different depending on the numerical experiment.We computed the fields on the unit square [0, 1] 2 with a 200 × 200 uniform grid at time t = 0.2 s.The integral (5) is approximated using the midpoint rule in time with 200 equal length subintervals of [0, t] and the trapezoidal rule on ∂Ω with 128 uniformly spaced points.A more detailed explanation of the numerical method appears in section 5.The accuracy of our numerical method with respect to discretization changes and noise is evaluated in section 2.3.Figure 3(c) shows a plot of the log 10 error between the computed field and the desired one.It can be observed that the accuracy of the numerical method improves as we move away from ∂Ω.
We point out that the term in (5) involving the initial condition is very different from the other terms because it is an integral over Ω rather than just the boundary ∂Ω.If the initial condition f (x) is non-zero but harmonic, this integral over Ω can also be expressed as an integral over ∂Ω by using Green identities as follows (see [33]) where φ(x, t) is given in two dimension by [33]: where Ein(z) = E(z) + ln z + γ, γ is the Euler constant and is the exponential integral (see e.g.[36, eq.6.2.2 and 6.2.4]).The expression for φ in three dimensions is given in [33].From the experimental perspective, it is not clear whether the kernels in the boundary integrals ( 6) can be achieved with heat monopole and dipole sources.So we assume from now on, that the initial condition is harmonic but does not need to be reproduced using sources on ∂Ω.
Under this assumption, we can simply subtract the initial condition to obtain the heat equation with zero initial condition, which is the case that we focus on.
Remark 2.2.As formula (5) uses the boundary data and the initial condition to express u Ω , this reconstruction of u satisfies some stability due the maximum principle applied to the heat equation (a property that does not hold e.g. for the wave equation).Indeed, for any T > 0, the maximum principle [37,38,39] states that a continuous function on Ω × [0, T ] that solves the homogeneous heat equation (2) on Ω × (0, T ) (in the distributional sense) reaches its minimum and maximum either at t = 0 or at any time t ∈ [0, T ] on the boundary ∂Ω.For instance, the solution of the initial-Dirichlet boundary value problem (with no sources) described in the chapter 7 page 171-172 of [39] satisfies the above conditions.In particular, this solution is continuous up to the boundary, i.e. on Ω × [0, T ].To obtain such continuity, the continuous initial data and the continuous Dirichlet data have to match on ∂Ω at t = 0, see e.g.[39].
To understand the stability, we take two solutions u j ∈ C 0 (Ω × [0, T ]), j = 1, 2, that satisfy the homogeneous heat equation (2) in Ω × (0, T ) for T > 0 in the distributional sense.In this setting solutions in the distributional sense are also smooth solutions of the homogeneous heat equation (2) (since by hypoellipticity, u j ∈ C ∞ (Ω × (0, T )) for j = 1, 2, see e.g.[32] theorem 1.1 page 192).Thus, one can apply the maximum principle (see [37], theorem 10.6 page 334) to obtain that (7) max Moreover if the initial conditions are harmonic, using the maximum principle for the Laplace equation gives Thus, in the space of solutions of the homogeneous heat equation (with the regularity described above), an error committed on the initial condition u(x, 0) or the boundary Dirichlet data of a solution u (i.e. the dipole distribution) on ∂Ω × [0, T ] controls the error (in the supremum norm) in the reconstruction of u Ω in Ω × (0, T ].Finally, we point out an important property that constrains the behavior of u 2 −u 1 if the maximum in ( 7) is attained at a point not located at the boundary or at the initial time.Under the additional assumption that Ω is connected, one shows based on the mean value property of the heat operator and a connexity argument (as in [38], section 2.2.3, theorem 4 page 54-55), that if there exists a point (x 0 , t 0 ) ∈ (0, T ] × Ω such that |u 2 − u 1 | reaches its maximum at (x 0 , t 0 ) in (7), then u 2 − u 1 has to be constant in Ω × [0, t 0 ].Indeed, for formula (8), this property holds also for t 0 = 0 since the initial condition is harmonic and the Laplace operator satisfies a mean value property.

2.2.
Reproducing fields exterior to a bounded region.For the exterior reproduction problem we seek to reproduce solutions to (2) outside of Ω by controlling sources on ∂Ω, while leaving the interior unperturbed.That is, for some temperature field v(x, t) solving the heat equation we wish to generate v Ω (x, t), such that ∈ Ω, for t > 0. Notice that we have used the subscript Ω differently in (9) than in (4).We adhere to the convention that u Ω always refers to the interior reproduction problem of heat equation solution u and v Ω refers to the exterior reproduction of a field v.
The problem of reproducing v Ω can be solved when the source term associated with v(x, t) is supported in Ω for all time.We only consider the case where v(x, 0) = 0. Non-zero initial conditions are left for future studies.Without further assumptions the exterior reproduction problem may not have a unique solution as can be illustrated by the one-dimensional non-uniqueness example by Tychonoff [40].Uniqueness for the Dirichlet problem can be guaranteed via the maximum principle [38, chapter 2, section 3, theorem 7] (see also remark 2.8) or via Sobolev regularity estimates in space and time [41,42].For the transmission problem, growth restrictions in the Laplace domain are used to prove uniqueness in [43].Here we want to establish a boundary representation formula for exterior solutions to the heat equation, which uses both Dirichlet and Neumann data.Such exterior representation formula has already been mentioned in [42,43,44], but without giving an explicit growth condition on the heat equation solution that guarantees its validity.We give a growth condition for the heat equation, analogous to the Sommerfeld radiation condition for the Helhmholtz equation [45] (a comparison between these two conditions is in remark 2.4).
To prove a boundary potential formula for the exterior reproduction problem, we follow the same steps as in [45] for the Helmholtz equation.Namely we use the interior reproduction problem (section 2.1) on the complement of Ω truncated to a ball of sufficiently large radius r (see fig. 4), and then give a sufficient condition guaranteeing that the contribution from the sources at |x| = r vanishes as r → ∞, allowing heat equation solutions outside of Ω to be reproduced by only controlling sources at ∂Ω.The sufficient condition that we impose on the growth of v(x, t) is close to the growth condition for the exterior Dirichlet problem uniqueness, see remark 2.4.
is said to satisfy the "growth condition" if there exists an r 0 > 0, such that if r > r 0 , where m is an integer, C > 0, a ≥ 0 are constants, the exponent b satisfies 0 ≤ b < 2 and S d (0, 1) is the sphere of radius 1 centered at the origin in d dimensions.
Remark 2.4.The bound in condition 2.3, is not as restrictive as the Sommerfeld radiation condition for the Helmholtz equation because of the Gaussian spatial decay of the heat kernel at fixed positive time.Indeed the Green function for the Helmholtz equation decays in space faster than r −m as r → ∞, where m > 0 depends on the dimension and r = |x|.On the other hand, the heat kernel decays faster than r m e −a r b as r → ∞, with a > 0, 0 ≤ b < 2 and m ∈ Z.This is how we motivate the bound in condition 2.3.However, we conjecture that it may be possible to improve the bound to allow b = 2 (Gaussian decay over a polynomial), which would bring it on the same par as the growth condition guaranteeing uniqueness for the Dirichlet problem outside of a bounded domain (remark 2.8).We point that condition 2.3 is satisfied by a large class of solutions of the heat equations (1) in (R d − Ω) × (0, +∞) that includes in particular the heat kernel and any of its spatial derivatives (see lemma 2.7).
, with zero initial condition and a source term spatially supported in the compact set Ω. Furthermore, assume that v satisfies the growth condition 2.3, then v can be reproduced for t > 0 in the exterior of Ω by the boundary representation formula where v Ω is as in (9).
Proof.Let (x, t) ∈ (R d − Ω) × (0, +∞) be fixed.Without loss of generality we assume that 0 ∈ Ω and take r to be large enough such that Ω ⊂ B(0, r) and ) satisfies the heat equation with a zero source term outside of Ω, we can use (5) with zero initial condition to reproduce v(x, t) on the bounded open set B(0, r) − Ω, giving where ( 13) The minus sign in I 2 (x, t) is because we defined n as the outward pointing unit normal to Ω.The goal is now to show that I 1 (x, t; r) → 0 as r → ∞, leaving us with only I 2 (x, t) which gives the desired result (11).We rewrite I 1 (x, t; r) using and switching the convolutions in time to get ( 14) • n v(y, t − s) K(x−y, s).
We define ξ = x/|x| (x = 0 since x / ∈ Ω).Thus, for y ∈ S(0, r), we can bound the heat kernel by ( 15) Noticing that we also have y • n = r for |y| = r, we can use condition 2.3 to bound I 1 (x, t; r) for sufficiently large r.Thus, using ( 14), the bound (15) and applying condition 2.3 leads to ( 16) where A d (r) is the surface of a sphere of radius r in d dimensions, which is given in terms of the Gamma function (see e.g.[36, eq.5.2.1]) by Now using the change of variables u = |x − rξ| 2 /4ks on the integral appearing in the right hand side of (16) yields: where the upper incomplete Gamma function Γ(d/2 − 1, •) is defined for all y > 0 by: In our case y = |x − rξ| d−2 /(4kt) → +∞ as r → +∞.Thus, we need an equivalent of Γ(d/2 − 1, y) as y → +∞.To this aim, we do an integration by parts on (18) to get : Then, as u ≥ y, one observes that and thus concludes from (19) that: Combining ( 16), (18) and the equivalent of the incomplete Gamma function (20) for y = |x − rξ| d−2 /(4kt) gives that for r large enough: where C x,t is positive constant that depends only x and t that are here fixed.To conclude, observe that the upper bound in (21) goes to 0 as r → ∞ since b < 2. This statement holds for any x outside of Ω and any t > 0, yielding the representation (11).
Remark 2.6.We point out that the regularity assumption v ∈ C 2 (R d − Ω) × [0, +∞)) of the solution in theorem 2.5 can be relaxed.Indeed, our proof still holds with weaker regularity assumptions but for a smooth bounded open set Ω (i.e. with a C ∞ boundary ∂Ω).For instance, our proof works under a Sobolev local regularity, namely if the solution v (in the sense of distributions) belongs to H 2,1 (O × (0, T )) for any T > 0 and any open bounded set O ⊂ R d − Ω (we refer to [34,46] for the definition of H 2,1 ).This assumption and the zero initial condition of v allows us to apply the representation formula (12) in the proof for O = B(0, r) − Ω by applying the theorem 2.20 of [34].In this setting the integrals in (12) have to be interpreted as duality pairings between Sobolev spaces of the boundary ∂Ω (see [34] for more details).By the trace theorem 2.1 page 9 in [46], the assumed local Sobolev regularity ensures that v and ∂v/∂n belong to L 2 (∂Ω × (0, t)) and L 2 (∂B r × (0, t)) for t > 0 and r sufficiently large.Thus, the integrals I 1 (x, t; r) and I 2 (x, t; r) can be interpreted not only as a duality paring but as integrals.Furthermore, by interior regularity (and even hypoellipticity) of the differential operator in the heat equation (see [32] theorem 1.1 page 192), one has v ∈ C ∞ (R d − Ω) × (0, +∞).Thus, as v is smooth on this set, the growth condition 2.3 is still well-defined.The proof of theorem 2.5 follows similarly and yields the representation formula (11) in this new setting.
The heat kernel and its spatial derivatives (which all solve the heat equation) satisfy the growth condition 2.3 as we see next.
Since the heat kernel K is in C ∞ (R d × (0, +∞)), an induction on the degree of differentiation reveals that for (x, t) ∈ Ω × (0, +∞), K and any of its spatial derivatives have the form, (22) ∂ α x K(x, t) = P (x 1 /t, x 2 /t, . . ., x d /t, 1/t) K(x, t), where P is a multivariate polynomial.In the particular case α = (0, 0, . . ., 0), we have ∂ α K = K and thus P = 1.Using the triangle inequality on the expression (22) of ∂ α x K leads to: , where Q is a polynomial which has the same monomial terms as P , but whose coefficients are given by the modulus of the coefficients of P .
Let r 0 ≥ 1 and |x| = r > r 0 .We set u = |x| 2 /t in the right hand of side of (23).
) and the coefficients of Q are positive, it follows from (23) and the expression (3) of 4k) for |x| > r 0 and t¿0.As Q is a polynomial, due to the exponential term, the right hand side of ( 24) is clearly bounded for u > 0, thus there is a constant C 1,α > 0 (depending only on α, d and r 0 ) such that: for |x| > r 0 and t¿0.Now, as the bound (24) holds for any α, one immediately deduces that there is a C 2,α such that ( 26) , r > r 0 and t > 0. Thus, by (26), (24) and the bound r/t ≤ r 2 /t = u (as r > r 0 ≥ 1), there is a C 3,α > 0 such that: ( 27) for any r > r 0 , ξ ∈ S d (0, 1) and t > 0.
Remark 2.8.As in remark 2.2 for bounded Ω, the representation formula (11) satisfies some stability due to the maximum principle.However, as R d − Ω is unbounded, it requires a bound that controls the growth of the functions when |x| → +∞, namely, one assumes that there exist A, a > 0 such that: for some finite T > 0. This last condition allows Gaussian growth and is similar to condition 2.3.More precisely, let v j ∈ C 0 (R d − Ω) × [0, T ] for j = 1, 2 be two solutions of the homogeneous heat equation in (R d − Ω) × (0, T ) in the distributional sense that satisfy the growth condition (28).Again by hypoellipticity (see [32] theorem 1.1 page 192), u j is indeed a smooth solution of the homogeneous heat equation on (R d − Ω) × (0, T ) for j = 1, 2 since it is C ∞ on this set.Then by the maximum principle (see e.g.[38, chapter 3, section 3, theorem 6]) one has: (29) sup Note that the proof of [38, chapter 3, section 3, theorem 6] is done on all of R d , but can be adapted to R d − Ω with Ω a bounded Lipschitz domain.Furthermore if the initial condition is harmonic and decays to 0 as |x| → +∞, using the maximum principle for the Laplace equation in unbounded domains, one can simplify (29) to include only surface terms in the right hand side Thus, in the space of solutions of the homogeneous heat equation (that satisfy (28) and the regularity described above), both (29) and (30) tell us that an error committed on the initial condition or on the Dirichlet boundary data controls the reconstruction error of v Ω on (R d − Ω) × (0, T ], in the supremum norm.Furthermore, uniqueness on (R d − Ω) × [0, T ] for the heat equation exterior Dirichlet problem follows from the maximum principle equality (29), provided the growth condition (28) and regularity assumptions hold.Moreover if (28) is satisfied for any t > 0, this uniqueness result extends to (R d − Ω) × [0, +∞), assuming the same regularity assumptions but with an infinite time.Finally, as in remark 2.2, under the additional assumption that the open set Furthermore, one shows that if one considers the formula (30), this last property holds also for t 0 = 0 and the constant C has to be zero (since in formula (30), one assumes that the initial conditions decay to 0 when |x| → +∞ which imposes that C = 0).
A numerical example of the exterior reproduction of a field can be seen in fig. 5.The details of the example are the same as those in fig.3, except the point source has been moved to (0.5, 0.55).

Numerical sensitivity study of field reproduction.
We study the sensitivity of the numerical approximation of the boundary representation formulas ( 5) and (11) to the following factors: (a) the spatial discretization (number of points on ∂Ω), (b) the temporal discretization (number of time steps) and (c) errors in the densities appearing in the boundary representation formulas.As can be seen in fig.3(c) and fig.5(c), the reproduction error peaks close to the boundary, so we decided to exclude a neighborhood of ∂Ω from the error measures we present.The numerical approximation of the boundary reproduction formulas is explained in detail in section 5.Here we keep the same domain Ω as in the examples of sections 2.1 and 2.2.The boundary integral representations were used to approximate the field on a 100 × 100 uniform grid of [0, 1] 2 and the thermal diffusivity was taken to be k = 0.2.
The first case we consider is that of the interior reproduction problem, i.e., when the source distribution is supported in R 2 − Ω.We expect using (5) that u Ω = u inside Ω and u Ω = 0 outside Ω.To evaluate the quality of the numerical we include code to generate it as supplementary material).This indicates that the maximum error is attained near ∂Ω in space (and in time at t = 0.2), conforming to the maximum principle (remarks 2.2 and 2.8).We can use here the maximum principle on the computed error because the numerical methods we use generates solutions to the heat equation (see section 5).approximation we make, we calculate the relative reproduction error on a slightly smaller domain Ω −s = B(x 0 , (1 − s)r), where we used the L 2 norm of a function over some set R, namely We also calculate the absolute error outside of a slightly larger domain Ω s = B(x 0 , (1 + s)r), i.e. err + (u; The L 2 norms appearing in the error quantities that we consider are approximated using Riemann sums on the 100 × 100 grid of [0, 1] 2 .The domains of interest, [0, 1] 2 −Ω s and Ω −s , are illustrated by the blue and red regions of fig.6 respectively. In our numerical experiments we chose s = .05to get a buffer annulus at ±5% of r.The field u is generated by a delta source δ(x, t).
Similarly for the exterior reproduction problem, where we want to reproduce a field v satisfying the heat equation with a source term supported in Ω and satisfying the radiation type boundary condition (10), we calculate the absolute interior error and the relative exterior error For the exterior reproduction studies, the field u is produced by a delta source located at x = (0.5, 0.55) and t = 0.

2.3.1.
Sensitivity to spatial discretization.In fig.7 we illustrate the changes in reproduction error for both the interior (fig.7 first row) and exterior (fig.7 second row) reproduction problems.For both studies a uniform discretization of ∂Ω is used and 1000 uniform time steps.While increasing the number of points on ∂Ω decreases the error in all cases, the decrease from 50 to 100 points is modest.We think this is due to the temporal discretization error being dominant.

2.3.2.
Sensitivity to temporal discretization.We report in fig.8 the change in reproduction error as we increase the number of time steps while keeping the number of uniformly spaced points used to discretize ∂Ω fixed and equal to 100.This is done for both the interior (fig.8 first row) and exterior (fig.8 second row) reproduction problems.For a fixed time the errors decrease with the number of time steps, as expected.
2.3.3.Sensitivity to errors in the densities.In practice it cannot be assumed that the field to be reproduced is perfectly known so we report in fig. 9 how the reproduction error is affected by errors in the monopole and dipole densities appearing in ( 5) and ( 11) when discretized with 1000 time steps and 100 points on ∂Ω.Say φ (n) ∈ R 100 is a vector representing the values of either the monopole or dipole density in one of the boundary representation formulas at time n∆t, where ∆t is the time step.We perturb φ (n) with a vector δφ (n) ∈ R 100 with independent identically distributed zero mean Gaussian entries with standard deviation being a fraction (3%) of φ (n)  2 .For clarity we only show the error for a single realization of the perturbation δφ (n) .The errors we observe in fig. 9 oscillate rapidly because we introduced a random perturbation at every single time step.As expected, the error we introduced in the densities increases the overall error at every single time step.We include code to generate the spatial distribution of the reconstruction errors as supplementary material (see README file), we observe that the maximum error over a time window is attained near ∂Ω as predicted by the maximum principle (remarks 2.2 and 2.8).Since the additive noise comes from perturbing the monopole interior error exterior error

Interior reproduction problem
Exterior reproduction problem and dipole densities, this perturbation is also a smooth solution of the heat equation (see section 5), which makes the maximum principle applicable.
Remark 2.9.In the numerical results we present, the exterior errors, relative and absolute, are smaller than the interior errors.We also see a transient effect in the absolute exterior error.We believe this is similar to the temperature distribution near a point source, which also increases and then decreases with time.Since the errors that we report in figs.7 to 9 are based on the L 2 norm, the maximum principle considerations of remarks 2.2 and 2.8 do not apply directly to these numerical experiments.

Cloaking
The goal here is to use the results from section 2 to cloak sources or objects inside a cloaked region, by placing sources on the surface of the region.By cloaking we mean that it is hard to detect the object or source from only thermal measurements made outside the cloaked region.The boundary representation formulas of section 2 give us the appropriate surface source distribution.We start in section 3.1 with the interior cloaking of a source, directly applying the boundary representation formula in section 2.2.The interior cloaking of an object is illustrated in section 3.2 by using the boundary representation formula in section 2.1.The boundary representation formulae impose restrictions on what can be cloaked and how.In either case the field must be known for all time and with no sources in the region where it is interior error exterior error

Interior reproduction problem
Exterior reproduction problem Figure 9. Influence of random perturbations added to the monopole and dipole densities on the reproduction error for the interior reproduction problem of a point source located at x 0 = (0, 0), t = 0 (first row) and for the exterior reproduction problem of a point source located at x 0 = (0.5, 0.55), t = 0 (second row), both with thermal diffusivity k = 0.2.Here ∂Ω is the circle of radius 0.25 centered at (0.5, 0.5).Since the errors for small times are large, we only show the errors for t ≥ 0.1.In orange: error with random perturbations.In blue: error obtained with the unperturbed densities.
reproduced.For the interior cloaking of a source, the temperature field generated by this source must also satisfy condition 2.3.

3.1.
Cloaking a source in an unbounded domain.Given certain kinds of localized heat source distributions, we can find an active surface surrounding the source so that the source cannot be detected by an observer outside the surface.Let v i (x, t) be a free space solution to the heat equation ( 2) with zero initial condition and compactly supported source distribution h(x, t).Let Ω be an open bounded set (with Lipschitz boundary) that contains the support of the source h(x, t) for t > 0.
In an analogy with wave problems, we call v i the "incident field" and we further assume that it satisfies the growth condition 2.3.By theorem 2.5, we can find monopole and dipole densities on ∂Ω so that the boundary representation formula (11) gives −v i outside of Ω and 0 inside Ω.We call this the cloaking field v c and it is given for t > 0 by In this manner the total field v tot = v i + v c is zero outside of Ω and equal to v i inside Ω.Because the active surface ∂Ω perfectly cancels the effect of the source h(x, t) for x / ∈ Ω, the source cannot be detected by an observer.Figure 5 shows a numerical example of v c .

3.2.
Cloaking passive objects in an unbounded domain.One way to detect an object in free space using only thermal measurements would be to generate an incident or probing field u i (x, t) with a source distribution h(x, t), i.e. a solution to the heat equation ( 2) in free space with zero initial condition and h as its source term.In the presence of an object, the total field is given by u tot = u i +u s , where u s is the field "scattered" by the object, borrowing terminology from the wave equation.The scattered field is produced by the interaction between the incident field and the object and depends on the properties of the object (boundary condition, heat conductivity, . ..).We point out that u s (x, 0) = 0 because u tot (x, 0) = u i (x, 0).Having u s = 0 reveals the presence of an object.In the following we assume that the object is "passive", meaning that the scattered field is linear in the incident field.In particular this means that u s = 0 when u i = 0. Examples of passive objects include objects with homogeneous linear boundary conditions (e.g.Dirichlet, Neumann or Robin) or objects with a heat conductivity that is different from that of the surrounding medium (see e.g.[47,48] for transmission problems for the heat equation).We point out that the object is assumed to be open with Lipschitz boundary.
The results in section 2 can be used to cloak a passive object R by placing it inside a cloaking region Ω (i.e. a bounded open set Ω with smooth boundary such that R ⊂ Ω) and makes this whole region invisible from probing incident fields u i generated by a source h spatially supported in R d − Ω.Indeed, by controlling monopoles and dipoles on ∂Ω, the region Ω and the object within can be made indistinguishable from a patch of homogeneous medium, from the perspective of thermal measurements outside of Ω.In a similar manner to section 3.1, the idea is to use (5) to cancel the incident field in Ω, while leaving the outside of Ω unperturbed.The cloaking field, u c , produced by this active surface ∂Ω is then for t > 0: In principle, this cloaking field can be used to perfectly cancel the incident field in Ω for all t > 0. Since the temperature of the "modified incident field": u i + u c is zero in Ω, the temperature field surrounding the object vanishes and no scattered field is produced.In practice, the field u i + u c in the vicinity of the object does not perfectly vanish, but we expect it to be sufficiently close to zero so that the scattered field u s is very small (because of linearity).
Our technique is illustrated with an object with homogeneous Dirichlet boundary conditions in fig.10.Here the field u i is generated by a point source at x = (0.9, 0.3) and t = 0.For the heat equation we took k = 0.2 and the cloaked region is Ω = B(x 0 , r) with x 0 = (0.5, 0.5) and r = 1/3.We computed the fields on the unit square [0, 1] 2 with a 200 × 200 uniform grid.The field u c is found by approximating the integral (5) using the midpoint rule in time with 600 equal length subintervals of [0, 0.5] and the trapezoidal rule on ∂Ω with 128 uniformly spaced points.A more detailed explanation, including how the scattered fields are calculated, is in section 5. We represent in figure fig. 10 the total fields respectively generated by the incident field u i (left column) and by the "modified incident field": u i + u c (right column).As can be seen in the right column, the temperature fields, outside of the cloaked region are indistinguishable from the incident field u i .We do not include a detailed error plot for this configuration, as the error is similar to the one we encountered when studying the interior reproduction problem.Remark 3.1.Here are two ways of dealing with active objects, i.e. that are not passive.First if the object produces a non-zero scattered field u s when u i = 0, the object acts as a source and u s can be cancelled using the technique in section 3.1.This presumes perfect knowledge of u s and u i .Second, if the object does not produce a scattered field when immersed in some harmonic field u 0 , we can use as a cloaking field u c (x, t) = −u i (x, t) + u 0 (x) for x ∈ Ω and u c (x, t) = 0 for x / ∈ Ω, instead of (32).An example of such an object would be one with a constant c = 0 Dirichlet boundary condition.By our assumption, the field u 0 (x) = c does not create any scattering, regardless of the shape of the object.

Mimicking
Another possible application of the boundary representation formulas ( 5) and ( 11) is to mimic sources or passive objects.This is done in two steps.First we cancel out the original source or suppress the scattering of the original object using a source distribution on a surface ∂Ω surrounding the object.Second we adjust the source distribution so that the object or source appears to the observer as another object or source.We illustrate this idea with two cases: making sources look like other sources (section 4.1) and making a passive object look like a different passive object (section 4.2).Other combinations are possible but are not presented here.4.1.Source mimicking.For source mimicking, we consider the problem where there is a compactly supported source distribution, f (x, t) which we seek to make appear as a different compactly supported source distribution, g(x, t) from thermal measurements outside of a region Ω (where Ω is an open bounded set with Lipschitz boundary).The support of both distributions is assumed to be contained in Ω for all t > 0. The two corresponding solutions to the heat equation (2) are v(x, t; f ) and v(x, t; g), and we further assume they satisfy condition 2.3.
Mimicking can be achieved by simultaneously canceling the field v(x, t; f ) outside Ω and adding v(x, t; g) also outside Ω.Both can be done using the results in section 3.1.That is, using (31) we can find a monopole and dipole density on ∂Ω that generates a field (33) v c (x, t) = 0 x ∈ Ω, v(x, t; g) − v(x, t; f ) x / ∈ Ω.In this way the field v(x, t; f ) + v c (x, t) is equal to v(x, t; g) outside of Ω, as desired.
A numerical example to illustrate the method is given in fig.11.Fields are calculated in [0, 1] 2 using a uniform grid of 200 by 200 points at t = 0.2 s.Here a point source at y (1) = (0.6, 0.4) and t = 0 is made to appear as a point source at y (2) = (0.39, 0.6) and t = 0, from thermal measurements outside of Ω. Figure 11 (a) and (b) represent the fields v(x, t; f ) and v(x, t; g), where f (x, t) = δ(x − y (1) , t) and g(x, t) = δ(x − y (2) , t). has been constructed by applying (33).An error plot is shown in fig.11(d) where the error is largest near the boundary as we expect based on the sensitivity analysis of section 2. We believe the diagonal line with small error in fig.11(d) is an artifact of our choice of sources.

4.2.
Mimicking a passive object.Consider a passive object R that is completely contained in an open set Ω and a source exterior to Ω which produces the incident field u i .The goal is to make R look like a different passive object, S, from thermal measurements outside Ω.To achieve this we can use linearity and both the interior and exterior boundary representation formulas to find monopole and dipole densities on ∂Ω producing a field (34) where v s (x, t) is the scattered field corresponding to the object S, included also in Ω, resulting from the incident field u i (x, t).In this fashion the total field is u i (x, t) + v s (x, t) outside of Ω and 0 inside of Ω.Its associated scattered field is v s (x, t) outside of Ω and 0 inside of Ω − R as desired.
To illustrate the method we consider a "kite" object with homogeneous Dirichlet boundary conditions and make it appear as a "flower" object with identical boundary conditions.Here the field u i is generated by a point source at x = (0.25, 0.5) and t = 0.For the heat equation we took k = 0.2 and the domain is Ω = B(x 0 , r) with x 0 = (0.5, 0.5) and r = 0.25.We computed the fields on the unit square [0, 1] 2 with a 200 × 200 uniform grid.The field u c is found by approximating the integral (5) using the midpoint rule in time with 180 equal length subintervals of [0, 0.05] and the trapezoidal rule on ∂Ω with 128 uniformly spaced points.A more detailed explanation, including how the scattered fields are calculated, is in section 5. Figure 12(a) and (b) show the scattered fields from two different objects and fig.12(c) shows the mimicked scattered field.Because we are approximating the fields numerically, the field u i + u c is very close to zero in Ω, but not exactly zero.The errors (which are not reported here) are larger near the boundary and decay as we move outwards, as we observed in section 2.
Remark 4.1.Although we have only shown mimicking of passive objects, the same techniques outlined in remark 3.1, could be used to suppress the field generated by an active object, which is part of what needs to be done to mimic an active object.

A simple numerical approach to potential theory for the heat equation
The boundary representation formulas ( 5) and ( 11) can expressed very efficiently in terms of the single and double layer potentials for the heat equation defined for t > 0 by as well as the corresponding boundary layer operators given for t > 0 by where φ(x, t) and ψ(x, t) are time dependent densities defined on ∂Ω.A full derivation of these and other potential theory operators for the heat equation and their properties can be found in [34,39,49].For a review of classic potential theory see e.g.[50,45,39,51].
With the potential theory notation, the interior boundary representation formula (5), with zero initial condition, becomes (35) u Whereas the exterior boundary representation formula (11) becomes Galerkin methods are commonly used to approximate the spatial integrals in eqs.(35) and (36), with a number of different approaches to deal with the integration in time.For instance time marching [33], time-space Galerkin methods [42,44] convolution quadrature [43] and collocation [52].For simplicity we opted for an approach based on the trapezoidal rule for the integration on ∂Ω and the midpoint rule for the time convolution.To be more precise, there are two convolutions that need to be evaluated in order to calculate the boundary representations in ( 35) and (36); a convolution in space and a convolution in time.For the spatial integration, trapezoidal rule is used with uniformly placed points on a parametric representation of ∂Ω.Since this amounts to integrating a periodic function, we can expect that the convergence rate of the trapezoidal rule depends explicitly on the rate of decay of the Fourier coefficients of the function to be integrated [53].Due to the smoothness of the heat kernel, exponential convergence is expected.For the integration in time, the convolutions are of the form t 0 ds g(s)f (t − s).
These convolutions are approximated with the midpoint rule as follows (37) which avoids evaluating f at t = 0.This is handy in our case because of the singularity of the heat kernel at (x, t) = (0, 0), which only occurs for boundary layer operators.After this space-time discretization, the resulting approximations are by construction an interpolation on the space-time grid of a finite distribution of monopoles and dipoles located on ∂Ω.These distributions are smooth solutions of the homogeneous heat equation on any open set of (R d − ∂Ω) × R. Moreover lemma 2.7 (see eq. ( 25)) ensures that they are bounded for t > 0 and x / ∈ B(0, r) for r large enough.Thus, as they satisfy eq. ( 28), we can apply the maximum principle over any finite time window [t 1 , t 2 ] (with 0 ≤ t 1 < t 2 ) on any closed set of R d − Ω that does not intersect ∂Ω.We leave an accuracy study of the numerical approximation we use to future work.In particular there are more accurate ways of dealing with the approximation for s ∈ [0, ∆t] than the midpoint rule we use, see e.g.[33,42,52,54].
5.1.Scattered field computation.In the case of inclusions, finding the scattered field requires the use of the boundary layer operators.For simplicity we only consider a homogeneous Dirichlet inclusion, R, i.e.where the temperature on ∂R is held constant at 0. Neumann inclusions, and inclusions with varying thermal diffusivity, require the introduction of other boundary integral operators (adjoint double layer and hypersingular [34]), but a similar numerical approach can be applied to this case.
The idea is to look for a scattered field u s of the form (36) outside of R. Since the temperature at ∂R is constant and equal to 0, the scattered field satisfies u s | ∂R = −u i | ∂R .Hence we know the Dirichlet data on ∂R in the representation formula (36), but not the Neumann data.We treat this as an unknown boundary density ψ in (38) u The latter formula is only valid for x / ∈ R. To obtain an integral equation on ∂R we take the limit of (38) as x approaches ∂R.The limit could be different if we approach ∂R from the inside or from the outside.The limits are given by the socalled jump relations.For x ∈ ∂R, the jump relation needed for the single layer potential is (39) lim z→x K 0 ψ(z, t) = Vψ(x, t), which holds for z tending towards ∂R from both the interior and exterior of R. The jump relations needed for the double layer potential are (40) where z → x + denotes approaching ∂R from the exterior of R and z → x − from the interior.Using these jump relations and (38) we have (41) lim z→x + u s (z, t) = lim z→x + K 1 (−u i (z, t)) − lim z→x + K 0 ψ(z, t) Rearranging terms yields a boundary integral equation for ψ We discretize (42) as a linear system where the unknown is ψ evaluated on a uniform grid of ∂R and of [0, T ].The boundary integral operators V and K in (42) are discretized using the trapezoidal rule in space and the midpoint rule in time.
For instance, V is approximated by a M N × M N matrix which is lower triangular by blocks, with each block of size N × N .Here M is the number of time steps and ∂R is approximated by a polygon with N sides: .
The V i are N × N matrices with entries (V i ) jk = k K(x j − x k , i∆t), where x k is the center of the k−th segment of length k .Clearly the matrix in ( 43) is guaranteed to be invertible if V 1 2 is invertible.Though we have not studied the invertibility of V 1 2 , we observe that it may become singular for a given spatial discretization if the temporal discretization is not fine enough.

Summary and perspectives
We proposed a strategy for active cloaking for the time-dependent parabolic heat (or mass, or diffusive light) equation.Similar to previous work for active cloaking for e.g. the Helmholtz or Laplace equation (e.g.modelling time-harmonic waves and thermostatic problems), our results rely on active sources coming from Green identities to reproduce solutions inside or outside a bounded domain.The idea is to use a source distribution on a closed surface to reproduce certain solutions to the heat equation inside the surface and the zero solution outside or vice versa.We give a growth condition which is sufficient to guarantee that a solution to the heat equation can be reproduced outside of a closed surface.We apply these theoretical results in four ways: interior cloaking of a source, interior cloaking of an object, source mimicking, and object mimicking.For the cloaking problems the idea is to find an active surface that surrounds the object or source to make the object or source undetectable by thermal measurements outside the surface.In the mimicking problems, instead of making the object undetectable, we make the source or object appear as a different source or object from the perspective of thermal measurements outside the cloak.Our solution to these problems inherits the limitations of the reproduction method, namely that the fields must be known for all time on the active surface surrounding the object or source we want to cloak or mimic.However the maximum principle guarantees some stability of our approach, as it is based on boundary representation formulas.This is a special feature of the heat equation.We illustrate our method with simple potential theory based simulations which are consistent with the heat equation and thus allow us to interpret the numerical errors using the maximum principle.Our study was limited to the case where the initial condition is zero or harmonic.It may be possible to use (6) to cancel out the initial condition with a boundary integral, but (a) it is not clear whether this mathematical construct has a physical interpretation and (b) this representation formula is only valid inside a bounded domain.These are questions we plan on exploring.Although we focused on the isotropic heat equation, it may be possible to carry out a similar cloaking strategy for anisotropic media and when an advection term is added to the heat equation.Our approach could also be tied to the active cloaking strategies for the Helmholtz equation in [18,17] by going to Fourier or Laplace domain in time, so that the active sources do not completely surround the object to achieve partial cloaking.For the Fourier domain, the physical interpretation is to study time-harmonic sources.Another open question is whether the growth condition we provided is optimal and how it is related to the uniqueness question for the exterior problem associated with different boundary condition types for the heat equation.
We note that this work could be also adapted to cloaking [55,56] and mimicking [57] of quantum matter waves.We believe this approach could be further generalized to the Fokker-Plank equation arising in statistical mechanics, which could applied to gravitational systems for which some passive cloaking theory has been proposed [58].Finally, as some scattering cancellation based cloaking approach has been proposed for Maxwell-Cataneo heat waves [59], we believe that our method for active cloaking could be also applied to such pseudo-waves.
Data Access.The Matlab code to reproduce the figures figs.Author Contributions.SG suggested the problem and applications to other physical phenomena.MC, TD, and FGV proved the exterior representation formula and the other theoretical results.All authors designed the numerical experiments.TD and FGV performed the numerical experiments.All authors contributed to the writing.
Funding.TD and FGV were supported by National Science Foundation Grant DMS-2008610.

Figure 3 .
Figure 3. Numerical example of the interior field reproduction problem for a point source located at (1/4, 1/4) and t = 0 s.A snapshot of the original field at time t = 0.2 s appears in (a).The field is reconstructed inside a disk of radius 1/4 centered at (1/2, 1/2) by using only heat sources on the corresponding circle.In (c) we show the log 10 of the reconstruction error (the absolute value of the difference between the exact u Ω and its numerical approximation).We generated a plot similar to (c) by taking the maximum over the time interval [0.2, 0.3] at each grid point (the plot being very similar to (c), we include the code to generate it as supplementary material).This indicates that the maximum error is attained near ∂Ω in space (and in time at t = 0.2 on the time interval [0.2, 0.3]), conforming to the maximum principle (applied with initial time t = 0.2) for the interior problem (remark 2.2) and the exterior one (remark 2.8).We can use the maximum principle on the computed fields because the numerical method we use generates solutions to the heat equation (see section 5 for more details).

Figure 4 .
Figure 4. To show that the exterior reproduction problem has a solution outside of a bounded region Ω we use the interior reproduction problem in the yellow region.Theorem 2.5 shows that for fields satisfying condition 2.3, the contribution of the sources on the sphere S d (0, r) of radius r vanishes as r → ∞ for d ≥ 2.

Figure 5 .
Figure 5. Numerical example of the exterior field reproduction problem.The details of this example are the same as fig.3, but with the point source moved to (0.5, 0.55).We generated a plot similar to (c) by taking the maximum over the time interval [0.2, 0.3] for each grid point (the plot being very similar to (c), we include code to generate it as supplementary material).This indicates that the maximum error is attained near ∂Ω in space (and in time at t = 0.2), conforming to the maximum principle (remarks 2.2 and 2.8).We can use here the maximum principle on the computed error because the numerical methods we use generates solutions to the heat equation (see section 5).

Figure 7 .
Figure 7. Influence of the number of points used to discretize ∂Ω on the reproduction error for the interior reproduction problem of a point source located at x 0 = (0, 0), t = 0 (top row) and for the exterior reproduction problem of a point source located at x 0 = (0.5, 0.55), t = 0 (bottom row), both with thermal diffusivity k = 0.2.Here ∂Ω is the circle of radius 0.25 centered at (0.5, 0.5).Since the errors for small times are large, we only show the errors for t ≥ 0.1.

Figure 8 .
Figure 8. Influence of the number of time steps on the reproduction error for the interior reproduction problem of a point source located at x 0 = (0, 0), t = 0 (first row) and for the exterior reproduction problem of a point source located at x 0 = (0.5, 0.55), t = 0 (second row), both with thermal diffusivity k = 0.2.Here ∂Ω is the circle of radius 0.25 centered at (0.5, 0.5).Since the errors for small times are large, we only show the errors for t ≥ 0.1.

Figure 10 .
Figure 10.Numerical example of the object cloaking problem, with object to hide having homogeneous Dirichlet boundary conditions.The incident field is produced by a point source at x = (0.9, 0.3) and t = 0s.The left images (a), (c) and (e) show time snapshots of the object without the cloak and the right images (b), (d) and (f) the corresponding snapshots when the cloak is active.Placing the cloak close to the object is a challenging simulation because we expect that errors will be highest near the boundary as in fig. 3. [See also movie in supplementary material]

Figure 11 .
Figure 11.Time snapshots from a numerical example of the source mimicking problem with point sources at t = 0.2s.In (c) the original point source in (a) is made to appear as the point source in (b) from the perspective outside of the cloaking region Ω.The point source in (a) is at x = (0.6, 0.4) and t = 0 and in (b) is at x = (0.39, 0.6) and t = 0.A plot of the log 10 errors in the exterior appears in (d).The diagonal line is an artifact due to the symmetry of the problem.A similar problem with point sources at x = (0.5, 0.4) and x = (0.5, 0.6) produces a horizontal line.

Figure 12 .
Figure 12.Numerical example of the object mimicking problem with a "kite" and "flower" object, both with homogenous Dirichlet boundary conditions.A snapshot the scattered field from the original object at time t = 0.05s appears in (a).In (c), the scattered field from the original object is made to appear as the scattered field from the object in (b) from the perspective of thermal measurements outside of the cloaking region Ω.

3 , 5
and 7 to 12 is available in repository [TBA].For example, fig. 3 can be reproduced by running the Matlab script fig3.m.The scripts for figs.3 and 5 generate additional plots that are mentioned in the respective figure captions.Figure 10 is animated in movie1.mp4.The spatial distribution for the monopole/dipole added noise (section 2.3) can be generated with the scripts fig9supint.m and fig9supext.m.
t) is a smooth solution of the heat equation on any open subset of this set.We consider a bounded non empty open set Ω with Lipschitz boundary ∂Ω and possibly multiple connected components in R d , d ≥ 2. The interior reproduction problem (section 2.1