Cloaking via mapping for the heat equation

This paper explores the concept of near-cloaking in the context of time-dependent heat propagation. We show that after the lapse of a certain threshold time instance, the boundary measurements for the homogeneous heat equation are close to the cloaked heat problem in a certain Sobolev space norm irrespective of the density-conductivity pair in the cloaked region. A regularised transformation media theory is employed to arrive at our results. Our proof relies on the study of the long time behaviour of solutions to the parabolic problems with high contrast in density and conductivity coefficients. It further relies on the study of boundary measurement estimates in the presence of small defects in the context of steady conduction problem. We then present some numerical examples to illustrate our theoretical results.

1. Introduction 1.1. Physical motivation. This work addresses the concept of near-cloaking in the context of the time-dependent heat propagation. Our study is motivated by experiments in electrostatics [LMMC12] and thermodynamics [SKGW13], which have demonstrated the markedly different behaviour of structured cloaks in static and dynamic regimes. It has been observed, in the physics literature, that the thermal field in [SKGW13] reaches an equilibrium state after a certain time interval, when it looks nearly identical to electrostatic field measured in [LMMC12]. There have already been some rigorous results for the electrostatic case [GLU03]. However, during the transient regime the thermal field looks much different from the static result, and a natural question that arises is whether one can give a mathematically rigorous definition of cloaking for diffusion processes in the time domain. This has important practical applications as cloaking for diffusion processes have been applied to mass transport problems in life sciences [GP13] and chemistry [ZS13] as well as to multiple light scattering whereby light is governed by ballistic laws [SKBW14]. Interestingly, the control of wave trajectories first proposed in the context of electromagnetism [PSS06] can also be extended to matter waves that are solutions of a Schrödinger equation [ZGSZ08] that is akin to the heat equation and so the near-cloaking we investigate has broad implications and application.
1.2. Analytical motivation. Change-of-variables based cloaking schemes have been inspired by the work of Greenleaf and co-authors [GLU03] in the context of electric impedance tomography and by the work of Pendry and co-authors [PSS06] in the context of time-harmonic Maxwell's equations. The transformation employed in both those works is singular and any mathematical analysis involving them becomes quite involved. The transformation in [GLU03,PSS06] essentially blows up a point to a region in space which needs to be cloaked. These works yield perfect-cloaks, i.e., they render the target region completely invisible to boundary measurements. Regularised versions of this singular approach have been proposed in the literature. In [KSVW08], Kohn and co-authors proposed a regularised approximation of this map by blowing up a small ball to the cloaked region and studied the asymptotic behaviour as the radius of the small ball vanishes, thus recovering the singular transform of [GLU03,PSS06]. An alternate approach involving the truncation of singularities was employed by Greenleaf and co-authors in [GKLU08] to provide an approximation scheme for the singular transform in [GLU03,PSS06]. It is to be noted that the constructions in [KSVW08] and [GKLU08] are shown to be equivalent in [KLS13]. We refer the interested reader to the review papers [GKLU09a,GKLU09b] for further details on cloaking via change of variables approach with emphasis on the aforementioned singular transform.
Rather than employing the singular transformation, we follow the lead of Kohn and co-authors [KSVW08]. This is in contrast with some works in the literature on the time-dependent thermal cloaking strategies where singular schemes are used e.g., [GAV12,PGB + 14,Pet15]. Note that the evolution equation which we consider is a good model for [SKGW13], which designs and fabricates a microstructured thermal cloak that molds the heat flow around an object in a metal plate. We refer the interested reader to the review paper [RLLM16] for further details on transformation thermodynamics. The work of Kohn and co-authors [KSVW08] estimates that the near-cloaking they propose for the steady conduction problem is ε d -close to the perfect cloak, where d is the space dimension. Our construction of the cloaking structure is exactly similar to the construction in [KSVW08]. In the present time-dependent setting, we allow the solution to evolve in time until it gets close to the associated thermal equilibrium state which solves a steady conduction problem. This closeness is studied in the Sobolev space H 1 (Ω). We then employ the ε d -closeness result of [KSVW08] to deduce our near-cloaking theorem in the timedependent setting. To the best of our knowledge, this is the first work to consider near-cloaking strategies to address time-dependent heat conduction problem.
In the literature, there have been numerous works on the approximate cloaking strategies for the Helmoltz equation [KOVW10,Ngu10,NV11,Ngu12] (see also [NV12] for the treatment of the full wave equation). The strategy in [NV12] to treat the time-dependent wave problem is to take the Fourier transform in the time variable. This yields a family of Helmholtz problems (family indexed by the frequency). The essential idea there is to obtain appropriate degree of invisibility estimates for the Helmholtz equation -the estimates being frequency-dependent. More importantly, these estimates blow-up in the frequency regime. But, they do so in an integrable fashion. Equipped with these approximate cloaking results for the Helmholtz equations, the authors in [NV12] simply invert the Fourier transform to read off the near-cloaking result for the time-dependent wave equation. Inspired by [NV12], one could apply the Laplace transform in the time variable for the heat equation and try to mimic the analysis in [NV12] for the family of elliptic problems thus obtained. Note that this approach does not require, unlike ours, the solution to the heat conduction problem to reach equilibrium state to obtain approximate cloaking. We have not explored this approach in detail and we leave it for future investigations.
Gralak and co-authors [GAA + 16] recently developed invisible layered structures in transformation optics in the temporal regime. Inspired by that work, we have also developed a transformation media theory for thermal layered cloaks, which are of practical importance in thin-film solar cells for energy harvesting in photovoltaic industry. These layered cloaks might be of importance in thermal imaging. We direct the interested readers to [AIKK05] and references therein.
In the applied mathematics community working on meta-materials, enhancement of near-cloaking is another topic which has been addressed in the literature. Loosely speaking, these enhancement techniques involve covering a small ball of radius ε by multiple coatings and then applying the push-forward maps of [KSVW08]. These multiple coatings which result in the vanishing of certain polarization tensors help us improve the ε d -closeness of [KSVW08] to ε dN -closeness where N denotes the number of coatings in the above construction. For further details, we direct the readers to [AKLL12b,AKLL12a] in the mathematics literature and [AE05,AE07] in the physics literature (the works [AE05, AE07] employ negative index materials). One could employ the constructions of [AKLL12b,AKLL12a] in the time-independent setting to our temporal setting to obtain enhanced near-cloaking structures. This again is left for future investigations.
In this present work, we are able to treat time-independent sources for the heat equation. As our approach involves the study of thermal equilibration, we could extend our result to time-dependent sources which result in equilibration. This, however, leaves open the question of near-cloaking for the heat equation with genuinely time-dependent sources which do not result in thermal equilibration. For example, sources which are time harmonic cannot be treated by the approach of this paper. The approach involving Laplace transform mentioned in a previous paragraph might be of help here. There are plenty of numerical works published by physicists on these aspects, but with no mathematical foundation thus far. The authors plan to return to these questions in the near future.
1.3. Paper structure. The paper is organized as follows. In section 2, we briefly recall the change-of-variable-principle for the heat equation. This section also makes precise the notion of near-cloaking and its connection to perfect cloaking followed by the construction of cloaking density and conductivity coefficients. Our main result (Theorem 2) is stated in that section. Section 3 deals with the long time behaviour of solutions to parabolic problems; the effect of high contrast in density and conduction on the long time behaviour of solutions is also treated in that section. The proof of Theorem 2 is given in section 4. In this section, we also develop upon an idea of layered cloak inspired by the construction in [GAA + 16]. Finally, in section 5, we present some numerical examples to illustrate our theoretical results.

Mathematical setting
Let Ω ⊂ R d (d = 2, 3) be a smooth bounded domain such that B 2 ⊂ Ω. Throughout, we use the notation B r to denote an euclidean ball of radius r centred at the origin.
2.1. Change-of-variable principle. The following result recalls the principle behind the change-of-variables based cloaking strategies. This is the typical and essential ingredient of any cloaking strategies in transformation media theory.

Proposition 1. Let the coefficients
where the coefficients are given as ; with the understanding that the right hand sides in (1) are computed at x = F −1 (y). Moreover we have for all t > 0, The proof of the above proposition has appeared in the literature in most of the papers on "cloaking via mapping" techniques. It essentially involves performing a change of variables in the weak formulation associated with the differential equation. We will skip the proof and refer the reader to [ Following Kohn and co-authors [KSVW08], we fix a regularising parameter ε > 0 and consider a Lipschitz map F ε : Ω → Ω defined below Note that F ε maps B ε to B 1 and the annulus B 2 \ B ε to B 2 \ B 1 . The cloaking strategy with the above map corresponds to having B 1 as the cloaked region and the annulus B 2 \ B 1 as the cloaking annulus. The Lipschitz map given above is borrowed from [KSVW08, page 5]. Remark that taking ε = 0 in (3) yields the map which is the singular transform of [GLU03,PSS06]. The map F 0 is smooth except at the point 0. It maps 0 to B 1 and B 2 \ {0} to B 2 \ B 1 .

Essential idea of the paper.
Let us make precise the notion of near-cloaking we will use throughout this paper. Let f ∈ L 2 (Ω) denote a source term such that supp f ⊂ Ω \ B 2 . Let g ∈ L 2 (∂Ω) denote a Neumann boundary data. Suppose the initial datum u in ∈ H 1 (Ω) be such that supp u in ⊂ Ω \ B 2 . Consider the homogeneous (conductivity being unity) heat equation for the unknown u hom (t, x) with the aforementioned data. (5) in Ω.
Here n(x) is the unit exterior normal to Ω at x ∈ ∂Ω.
Our objective is to construct coefficients ρ cl (x) and A cl (x) such that for some arbitrary bounded positive density η and for some arbitrary bounded positive definite conductivity β. This construction should further imply that the evolution for the unknown u cl (t, x) given by in Ω, is such that there exists a time instant T < ∞ so that for all t ≥ T , we have The above closeness will be measured in some appropriate function space norm.
Most importantly, this approximation should be independent of the density-conductivity pair η, β in B 1 . Note, in particular, that the source terms f (x), g(x), u in (x) in (5) and (6) are the same.
2.3. Cloaking coefficients & the defect problem. The following construct using the push-forward maps is now classical in transformation optics and we use it here for thermodynamics.
The density coefficient η(x) in (7) is any arbitrary real coefficient such that The conductivity coefficient β(x) in (8) is any arbitrary bounded positive definite matrix, i.e., there exist positive constants κ 1 and κ 2 such that The following observation is crucial for the analysis to follow. Consider the densityconductivity pair Next let us compute their push-forwards using the Lipschitz map F ε -as given by the formulae (1) -yielding Then, the assertion of Proposition 1 (see in particular the equality (2)) implies that the solution u cl (t, x) to (6) satisfies for all t > 0, with u ε (t, x) being the solution to with the coefficients in the above evolution being given by (9)-(10). The coefficients ρ ε and A ε are uniform except for their values in B ε . Hence we treat B ε as a defect where the coefficients show high contrast with respect to their values elsewhere in the domain. Due to the nature of these coefficients, we call the evolution problem (11) the defect problem with high contrast coefficients or defect problem for short. The change-of-variable principle (see Proposition 1) essentially says that, to study cloaking for the transient heat transfer problem, we need to compare the solution u ε (t, x) to the defect problem (11) with the solution u hom (t, x) to the homogeneous problem (5) for x ∈ Ω \ B 2 .

Main result.
We are now ready to state the main result of this work.
x) be the solution to the defect problem (11) and let u hom (t, x) be the solution to the homogeneous conductivity problem (5). Suppose the data in (5) and (11) are such that Let us further suppose that the source terms satisfy Then, there exists a time T < ∞ such that for all t ≥ T we have where the positive constant C depends on the domain Ω and the L ∞ bounds on the density-conductivity pair (η, β) in B 1 .
Thanks to the change-of-variable principle (Proposition 1), we deduce the following corollary.
Corollary 3. Let u cl (t, x) be the solution to the thermal cloak problem (6) with the cloaking coefficients ρ cl (x), A cl (x) given by (7)-(8). Let u hom (t, x) be the solution to the homogeneous conductivity problem (5). Suppose the data in (5) and (6) are such that Let us further suppose that the source terms satisfy (12) and (13). Then, there exists a time T < ∞ such that for all t ≥ T we have where the positive constant C depends on the domain Ω and the L ∞ bounds on the density-conductivity pair (η, β) in B 1 .
Whenever the source terms f (x), g(x) and the initial datum u in (x) satisfy the compatibility conditions (12)-(13), we shall call them admissible. Our proof of the Theorem 2 relies upon two ideas (i) long time behaviour of solutions to parabolic problems.
(ii) boundary measurement estimates in the presence of small inhomogeneities. (5) and (11). Those elliptic boundary value problems demand that the source terms be compatible to guarantee existence and uniqueness of solution. These compatibility conditions on the bulk and Neumann boundary sources f, g translates to the zero mean assumption (12). Note that the compatibility condition (12) is not necessary for the well-posedness of the transient problems.

Remark 4. Our strategy of proof goes via the study of the steady state problems associated with the evolutions
A remark about the assumption (13) on the initial datum u in in Theorem 2 will be made later after the proof of Theorem 2 as it will be clearer to the reader.

Study of the long time behaviour
In this section, we deal with the long time asymptotic analysis for the parabolic problems. Consider the initial-boundary value problem for the unknown v(t, x). (16) We give an asymptotic result for the solution to (16) in the t → ∞ limit.

Proposition 5. Let v(t, x) be the solution to the initial-boundary value problem
where v in denotes the initial average, i.e., Proof of the above proposition is standard and is based in the spectral study of the corresponding elliptic problem, see e.g. [Paz83]. Alternatively, we could also use energy methods based on a priori estimates on the solution of the parabolic partial differential equations, see e.g. [Chi00, Chapter 13]. As the proof of Proposition 5 is quite standard, we omit the proof.
Let us recall certain notions necessary for our proof. Let {ϕ k } ∞ k=1 ⊂ H 1 (Ω) denote the collection of Neumann Laplacian eigenfunctions on Ω, i.e., for each k ∈ N, with {µ k } denoting the eigenvalues. Recall that the spectrum is discrete, nonnegative and with no finite accumulation point, i.e., The next proposition is a result similar in flavour to Proposition 5, but in a more general parabolic setting with high contrast in density and conductivity coefficients. More precisely, let us consider the heat equation with uniform conductivity in the presence of a defect with high contrast coefficients. We particularly choose the conductivity matrix to be (10) and the density coefficient to be (9). For an unknown v ε (t, x), consider the initial-boundary value problem We give an asymptotic result for the solution v ε (t, x) to (19) in the t → ∞ limit.
As before, we are interested in the · H 1 -norm rather than the · L 2 -norm.

Proposition 6.
Let v ε (t, x) be the solution to the initial-boundary value problem where m ε denotes the following weighted initial average Remark 7. From estimate (20) we conclude that the estimate on the right hand side is a product of an exponential decay term and a term of O(ε −d ). So, if γ ε 1, uniformly in ε, then the solution converges to the weighted initial average in the long time regime. We will demonstrate that the decay rate γ ε is bounded away from zero uniformly (with respect to ε) in subsection 3.1. If the density coefficient ρ ε (x) is given by (9), then where Γ(·) denotes the gamma function. Substituting for ρ ε in the expression for the weighted initial average yields Using the assumption on the initial datum v in that it belongs to H 1 (Ω) ∩ L ∞ (Ω), we get the following uniform bound (uniform with respect to ε) on the weighted initial average Proof of Proposition 6. Let µ ε k and ϕ ε k be the Neumann eigenvalues and eigenfunctions defined as where the conductivity-density pair A ε (x), ρ ε (x) is given by (10)-(9). Here again the spectrum is discrete, non-negative and with no finite accumulation point, i.e., The general representation formula for the solution to (19) becomes The first term in the above representation is nothing but the weighted initial average where we have used the following notation for the weighted Lebesgue space norm: The density coefficient ρ ε (x) defined in (9) appears as the weight function in the above Lebesgue space. It follows that thanks to the definition of ρ ε (x) in (9). Computing the weighted norm of the Next, it follows from the representation formula (22) that Repeating the above computations for the difference v ε (t, x) − m ε , we obtain Gathering the inequalities (23)-(24) together proves the proposition with the constant γ ε = µ ε 2 , i.e., the first non-zero Neumann eigenvalue in the spectral problem (21) on Ω with the high contrast density-conductivity pair (ρ ε (x), A ε (x)).
3.1. Reduction to a Schrödinger operator. The constants γ and γ ε in Propositions 5 and 6 respectively give the rate of convergence to equilibrium. As the proofs in the previous subsection suggest, these rates are nothing but the first non-zero eigenvalues of the associated Neumann eigenvalue problems.
The constant decay rate γ ε in Proposition 6 may depend on the regularisation parameter ε. In our setting, ε is nothing but the radius of the inclusion. Hence, to understand the behaviour of γ ε in terms of ε, we need to study the perturbations in the eigenvalues caused by the presence of inhomogeneities with conductivities and densities different from the background conductivity and density. More importantly, we need to understand the spectrum of transmission problems with high contrast conductivities and densities. The spectral analysis of elliptic operators with such conductivity matrices is done extensively in the literature in the context of electric impedance tomography -see [AM03,AKL09] and references therein for further details. In [AM03], for example, the authors give an asymptotic expansion for the eigenvalues µ ε k in terms of the regularising parameter ε (even in the case of multiplicities). We refer the readers to [AM03, Eq.(23), page 74] for the precise expansion. For high contrast conductivities such as A ε , refer to the concluding remarks in [AM03, pages 74-75] and references therein. Another important point to be noted is that the above mentioned works of Ammari and co-authors do not treat high contrast densities while addressing the spectral problems. For our settingmore specifically for the spectral problem (21) -the readers are directed to consult the review paper of Chechkin [Che06] which goes in detail about the spectral problem in a related setting. Furthermore, the review paper [Che06] gives exhaustive reference to literature where similar spectral problems are addressed. Rather than deducing the behaviour of the first non-zero eigenvalue of the spectral problem (21) from [Che06], we propose an alternate approach. Studying (21) is the same as the study of the spectral problem for the operator The idea is to show that the study of the spectral problem for L is analogous to the study of the spectral problem for a Schrödinger-type operator. The essential calculations to follow are inspired by the calculations in [Pav14, section 4.9, page 125]. Note that the operator L defined by (25) with zero Neumann boundary condition is a symmetric operator in L 2 (Ω; ρ ε ), i.e., Let us now define the operator with the coefficient

An algebraic manipulation yields
In our setting, with the high contrast coefficients A ε and ρ ε from (10)-(9), the coefficients Σ ε and W ε become By definition (26), the operators L and H are unitarily equivalent. Hence, they have the same eigenvalues. Note that the operator H is a Schrödinger-type operator where the coefficients are of high contrast. By the Rayleigh-Ritz criterion, we have the characterisation where ϕ ≡ 0 and is orthogonal to first k−1 Neumann eigenfunctions ϕ ε 1 , . . . , ϕ ε k−1 . Note, in particular, for k = 2 (i.e. the first non-zero eigenvalue) we have Note that we have used the fact that W ε (x) is supported on B ε . Note further that if we assume β ∈ W 1,∞ (B 1 ; R d×d ) and η ∈ W 2,∞ (B 1 ). Note further that the spectral problem (18) comes with the following characterisation of the first non-zero eigenvalue (again by the Rayleigh-Ritz criterion) Let us now take the test function ϕ to be the normalised eigenfunction ϕ 2 (x) associated with the first non-zero eigenvalue µ 2 for the Neumann Laplacian: Using the observation (28) that the potential W ε is of O(1) and that the Neumann eigenfunctions are bounded in W p,∞ for any p < ∞ [Gri02], we have proved that In this subsection, we have essentially proved the following result.
Proposition 8. Suppose that the high contrast conductivity-density pair A ε , ρ ε is given by (10)-(9) with β ∈ W 1,∞ (B 1 ; R d×d ) and η ∈ W 2,∞ (B 1 ). Let µ ε 2 and µ 2 be the first non-zero eigenvalues associated with the Neumann spectral problems (21) and (18) respectively. Then we have Hence, as a corollary to the above result, we can deduce that the decay rate for the homogeneous transient problem (16) and that for the high contrast transient problem (19) are close to each other in the ε 1 regime.

Near-cloaking result
In this section, we will prove the main result of this paper, i.e., Theorem 2. To that end, we first consider the steady-state problem associated with the homogeneous heat equation (5). More precisely, for the unknown u eq hom (x), consider Note that the last line of (31) is to ensure that we solve for u eq hom (x) uniquely. Here we assume that the source terms are admissible in the sense of (12). This guarantees the solvability of the above elliptic boundary value problem. Next we record a corollary to Proposition 5 which says how quickly the solution u hom (t, x) to the homogeneous problem (5) tends to its equilibrium state.
Corollary 9. Let u hom (t, x) be the solution to (5) and let u eq hom (x) be the solution to the steady-state problem (31). Suppose the source terms f (x), g(x) and the initial datum u in (x) are admissible in the sense of (12)-(13). Then for some positive constant γ.
Proof. Define a function w(t, x) := u hom (t, x) − u eq hom (x). We have that the function w(t, x) satisfies the evolution equation in Ω, which is the same as (16). The estimate (32) is simply deduced from Proposition 5 (see in particular (17)).
Now we record a result, as a corollary to Proposition 6, demonstrating how quickly the solution u ε (t, x) to the defect problem with high contrast coefficients (11) tends to its equilibrium state. Consider the steady-state problem associated with the defect problem (11). More precisely, for the unknown u ε eq (x), consider the elliptic boundary value problem Note that the normalisation condition in (33) makes use of the density coefficient ρ ε (x) defined by (9).
Corollary 10. Let u ε (t, x) be the solution to (11) and let u ε eq (x) be the solution to the steady-state problem (33). Suppose the source terms f (x), g(x) and the initial datum u in (x) are admissible in the sense of (12)-(13). Suppose further that Proof. Define a function w ε (t, x) := u ε (t, x) − u ε eq (x). We have that the function w ε (t, x) satisfies the evolution equation in Ω which is the same as (19). By assumption, the initial datum u in is supported away from B 2 . This along with the admissibility assumption (13) implies that the weighted average m ε defined in Proposition 6 vanishes. Then the estimate (34) is simply deduced from Proposition 6 (see in particular (20)).

Proof of Theorem 2. From the triangle inequality, we have
.
The boundary trace inequality gives the existence of a constant c Ω -depending only on the domain Ω -such that Using the result of Corollary 9, Corollary 10 and the DtN estimate from [FV89, Lemma 2.2, page 305] (see also [KSVW08, Proposition 1]), we get that the right hand side of the above inequality is bounded from above by Hence the existence of a time instant T follows such that for all t ≥ T , we indeed have the estimate (14).
Remark 11. We make some observations on why the admissibility assumption (13) and the assumption on the support of the initial datum in Corollary 10 were essential to our proof. In the absence of these assumptions, in the proof of Theorem 2, we will have to show that Let us compute the difference We have not managed to characterise all initial data that guarantee the asymptote (35). The difficulty of this task becomes apparent if you take initial data to be supported away from B 1 , but not satisfying the zero mean assumption (13). Note that our assumption on the initial data in Corollary 10 guarantees the above difference is always zero, irrespective of the value of the parameter ε.
4.1. Layered cloaks. Advancing an idea from [GAA + 16], we develop a transformation media theory for thermal layered cloaks, which are of practical importance in thin-film solar cells for energy harvesting in the photovoltaic industry. The basic principle behind this construction is the following observation.
Proposition 13. Let the spatial domain Ω := (−3, 3) 2 . Let the density conductivity pair ρ ∈ L ∞ (Ω; R), A ∈ L ∞ (Ω; R 2×2 ) be such that they are (−3, 3)-periodic in the where the coefficients are given below with the understanding that the right hand sides in (36) are computed at (x 1 , x 2 ) = (y 1 , f −1 (y 2 )). Furthermore, we have Next, we prove a near-cloaking result in this present setting of layered cloaks. It concerns the following evolution problems: The homogeneous problem (37) in Ω and the layered cloak problem where the coefficients ρ cl and A cl in (38) are defined using the Lipschitz mapping The precise construction of the layered cloaks is as follows where the push-forward maps are defined in (36). The density coefficient η(x) in (40) is any arbitrary real positive function. The conductivity coefficient β(x) in (41) is any arbitrary bounded positive definite matrix. Let us make the observation that the cloaking coefficients ρ cl and A cl given by (40) and (41) respectively can be treated as push-forward outcomes (via the pushforward maps (36)) of the following defect coefficients It then follows from Proposition 13 that comparing u hom and u cl is equivalent to comparing u hom and u ε where u ε (t, x) solves the following defect problem with the aforementioned ρ ε and A ε as coefficients: in Ω Theorem 14. Let u ε (t, x) be the solution to the defect problem (44) with high contrast coefficients (42)-(43) and let u hom (t, x) be the solution to the homogeneous conductivity problem (37). Then, there exists a time instant T < ∞ such that for all t ≥ T we have The proof is similar to the proof of Theorem 2. More specifically, we show first that the solutions to the transient problems (44) and (37) converge exponentially fast to their corresponding equilibrium states. We can then adapt the energy approach in the proof of [FV89] to show that the equilibrium states are ε 2 close in the H 1 2 (Γ)-norm. We note that our analysis of near-cloaking for thermal layered cloaks can be easily adapted to electrostatic and electromagnetic cases. It might also find interesting applications in Earth science for seismic tomography, in which case one could utilise results in [AKKL13] to prove near cloaking results related to imaging the subsurface of the Earth with seismic waves produced by earthquakes.

Numerical results
This section deals with the numerical tests done in support of the theoretical results in the paper. The tests designed in the subsections to follow make some observations with regards to the near-cloaking scheme designed in the previous sections of this paper. The numerical simulations are done in one, two and three spatial dimensions. It seems natural to start with the one dimensional case, but as it turns out, see subsection 4.1, the physical problem of interest for a one-dimensional (so-called layered) cloak requires a two-dimensonal computational domain, and thus we start with the 2D case. We refer the reader to [GAA + 16] for the precise physical setup and importance of the layered cloak described in the context of Maxwell's equations (this is easily translated into the language of conductivity equation). These numerical simulations were performed with the finite element software COMSOL MULTIPHYSICS.
We choose the spatial domain to be a square Ω := (−3, 3) 2 . We take the bulk source f (x), the Neumann datum g(x) and the initial datum u in (x) to be smooth and such that supp f ⊂ Ω \ B 2 , supp u in ⊂ Ω \ B 2 . We further assume that This guarantees that the data is admissible in the sense of (12)-(13).

Near-cloaking.
The numerical experiments in this subsection analyse the sharpness of the near-cloaking result (Theorem 2) in this paper. We solve the initial-boundary value problem for the unknown u ε (t, x): with the density-conductivity coefficients: Note that in two-dimensions, there is no high contrast in the conductivity coefficient A ε (x). There is, however, contrast in the density coefficient ρ ε (x). The evolution (46) is supplemented by the Neumann datum and the initial datum The bulk force in (46) is taken to be We next solve the initial-boundary value problem (Neumann) for u hom (t, x) with the above data.
Let us define We compute the function G ε (t) defined by (47) as a function of time for various values of ε, see Figure 1, and numerically observe that after 110s, G ε (t) reaches an asymptote that tends towards a numerical value close to 8.6 when ε gets smaller, in agreement with theoretical predictions of Theorem 2 for space dimension d = 2. The simulations were performed using an adaptive mesh of the domain Ω consisting of 2×10 5 nodal elements (with at least 10 2 elements in the small defect when ε = 10 −4 ), and a numerical solver based on the backward differential formula (BDF solver) with initial time step 10 −5 s, maximum time step 10 −1 s, minimum and maximum BDF orders of 1 and 5, respectively. Many test cases were ran for various values of β (including anisotropic conductivity) and η in the small inclusion, and we report some representative curves in Figure 1.
Insert shows a zoom-in.

Contour plots.
Let us start with the contour plots for the layered cloak developed in subsection 4.1, see Figure 2. More precisely, we compare solutions u hom and u cl to the evolution problems (37) and (38) respectively, where periodic boundary conditions are imposed at the boundary x 1 = ±3 and homogeneous Neumann datum at the boundary x 2 = ±3. Here, we illustrate the layered cloak by numerically solving an equation for the unknown u(t, x 1 , x 2 ), but the conductivity and density only depend upon the x 2 variable. We take the source f (x 1 , x 2 ) = x 2 sin(x 2 ) for x 2 ∈ (−3, −2) ∪ (2, 3), 0 for x 2 ∈ (−2, 2) and the initial datum is taken to be Following the layered cloak construction in (41), we take the cloaking conductivity to be the push-forward of identity in the cloaking strip: We report in Figure 2 some numerical results that exemplify the high level of control of the heat flux with a layered cloak: in the upper panels (a), (b), (c), one can see snapshots at representative time steps (t=0, 1 and 4s) of a typical one-dimensional diffusion process in a homogeneous medium for a given source with a support outside x 2 ∈ (−2, 2). When we compare the temperature field at initial time step in (a) with that when we replace homogeneous medium by a layered cloak in 1 < |x 2 | < 2 in (d), we note no difference. However, some noticeable differences are noted for the temperature field between the homogeneous medium and the cloak when comparing (b) with (d) and (c) with (e). The gradient of temperature field is dramatically decreased in the invisibility region x 2 ∈ (−1, 1), leading to an almost uniform (but non zero) temperature field therein, and this is compensated by an enhanced gradient of temperature within the cloak in 1 < |x 2 | < 2. One notes that the increased flux in 1 < |x 2 | < 2 might be useful to improve efficiency of solar cells in photovoltaics. Our next experiment is to compare the homogeneous solution and the cloaked solution where the cloaking coefficients are constructed using the push-forward maps as in (7)-(8). The data (f , g, u in ) are chosen as in subsection 5.1. We report in Figure 3 some numerical computations performed in COMSOL that illustrate the strong similarity between the temperature fields in homogeneous (a-d), small defect (e-h) and cloaked (i-l) problems. Obviously, the fields are identical outside B 2 at initial time step, then they differ most outside B 2 at small time step t=1s, see (b),(f),(j), and become more and more similar with increasing time steps, see (c), (g), (k) and (d), (h), (l). These qualitative observations are consistent with the near-cloaking result, Theorem 2, of this paper.
We have also performed a similar experiment in three dimensions, see Figure  4. Refer to the caption in Figure 4 for the parameters considered in the three dimensional problem. Note that for these 3D computations, we mesh the cubical domain with 40, 000 nodal elements, we take time steps of 0.1s and use BDF solver with initial time step 0.01s, maximum time step 0.1s, minimum and maximum BDF orders of 1 and 3, respectively. We use a desktop with 32 Gb of RAM memory and a computation run takes around 1 hour for a time interval between 0.1 and 10s. We are not able to study long time behaviours, nor solve the high contrast small defect problem in 3D, as this would require more computational resources. Nevertheless, our 3D computations suggest that there is a strong similarity between the temperature fields in homogeneous and cloaked problems, outside B 2 . Note also that we consider a source not vanishing inside B 2 , which motivates further theoretical analysis for sources with a support in the overall domain Ω. 5.3. Cloaking coefficients. The cloaking coefficients ρ cl and A cl in the annulus B 2 \ B 1 play all the essential roles in thermal cloaking phenomena. So, it is important in practice -to gain some physical intuition and start engineering and manufacturing processes of a meta-material cloak -to analyse the coefficients defined by (7)-(8). For readers' convenience we recall them below (only for the part Ω \ B 1 as coefficients inside B 1 can be arbitrary) for y ∈ B 2 \ B 1 Figure 2. Contour plots of a layered cloak (with a source outside 3)} and f (x) = 0 for x 2 ∈ (−2, 2). Upper panel: Plots of u for data u in (x) = x 2 and g(x) = 0 for x 2 = ±3 and such that u(−3, x 2 ) = u(3, x 2 ) for a homogeneous medium with diffusivity A = 1 at time steps t = 0s (a), 1s (b) and 4s (c). Lower panel: Same for the medium with diffusivity A = 1 outside x 2 ∈ (−2, 2) and a layered cloak inside x 2 ∈ (−2, 2) with diffusivity A(x 2 ) = diag A 11 (x 2 ), 1 Figure 3. Contour plots of a 2D cloak (with a source outside B 2 ): Upper panel: Plots of u for data u in (x 1 , x 2 ) = x 1 x 2 and g(x 1 , x 2 ) = −3 for x 1 = ±3 and 0 for x 2 = ±3 and diffusivity A = 1 at time steps t = 0s, 1s, 4s 10s and 20s. Middle panel: same for a small defect of radius ε = 10 −1 , density ρ ε = 2ε −2 and diffusivity A ε = 2 Id. Lower panel: Same for a cloak.
Bear in mind that only the radial coordinate r gets transformed by F ε and the angular coordinate θ remains unchanged. Reformulating the push-forward maps in terms of the polar coordinates yield the following expressions for the cloaking coefficients in the annulus: where A 11 (r ) is given by for r ∈ [1, 2] (50) and the rotation matrix We plot in Figure 5 the matrix entry A 11 (r ) given above and the push-forward density ρ 2D cl given in (49). We observe that when ε gets smaller, the radial conductivity and density take values very close to zero near the inner boundary of the cloak, which is unachievable in practice (bear in mind that the azimuthal conductivy being the inverse of the radial conductivity, the conductivity matrix becomes extremely anisotropic). Therefore, manufacturing a meta-material cloak would require a small enough value of epsilon so that homogenization techniques can be applied to approximate the anisotropic conductivity with concentric layers of isotropic phases, e.g. as was done for the 2D meta-material cloak manufactured at the Karlsruher Institut für Technologie [SKGW13].
In three spatial dimensions, we can recast the cloaking coefficients in the spherical coordinates (r, θ, ϕ). As in the cylindrical coordinate setting, only the radial variable gets modified by the Kohn's transformation F ε and the variables θ, ϕ remain unchanged. The transformed radial coordinate r is given by (48). The push-forward maps of interest for cloaking are (51) for r ∈ (1, 2) Note that the matrix entry B(r ) in the push-forward conductivity coincides with the push-forward density ρ 3D cl . We plot in Figure 6 the radial conductivity and density in the cloaking annulus given in (51). One should recall that the 3D spherical cloak only has a varying radial conductivity, the other two polar and azimuthal diagonal entries of the conductivity matrix being constant (i.e. independent of the radial position). Besides, the radial conductivity has the same value as the density within the cloak. All this makes the 3D spherical cloak easier to approximate with homogenization techniques. However, no thermal cloak has been manufactured and experimentally characterised thus far, perhaps due to current limitations in 3D manufacturing techniques; a possible route towards construction of a spherical cloak is 3D printing. Similar computations in spherical coordinates, but for Pendry's singular transformation are found in [NZAG08], where the matrix M was introduced to facilitate implementation of perfectly matched layers, cloaks and other transformation based electromagnetic media with spherical symmetry in finite element packages, see also [NRM + 94] that predates the field of transformational optics. −div A ε (x)∇φ ε = µ ρ ε (x) φ ε in Ω := (−3, 3) d , ∇φ ε · n = 0 on ∂Ω.
The coefficients are of high contrast and take the form The first non-zero eigenvalue for the Neumann Laplacian is µ 2 = (π/6) d . We illustrate the result of Proposition 8 by showing that the first non-zero eigenvalue µ ε 2 for the defect spectral problem is ε d -close to µ 2 for various value of ε and in dimensions one, two and three. The results are tabulated in Table 5.4. The associated eigenfields are also plotted in Figure 7. Note that the spectral problems in 1D, 2D and 3D are solved with direct UMFPACK and PARDISO solvers using adaptive meshes with 250, 150 and 50 thousand elements, respectively. We made sure that there are at least 10 2 elements in the small defect for every spectral problem solved.  Table 1. Numerical estimate of the difference | µ 2 − µ ε 2 |, versus the parameter ε = 10 −m with m = 1, · · · , 7 for dimension d = 1, with m = 1, · · · , 3 for d = 2 and with m = 1 for d = 3. Same source, Neumann data, diffusivity and density parameters as in Figure 7. Numerical results for contour plots of eigenfields φ and φ ε associated with µ 2 (a) and µ ε 2 (b-d), in dimension 1 (a,b), 2 (c) and 3 (d).

Concluding remarks
This work addressed the question of near-cloaking in the time-dependent heat conduction problem. The main inspiration is derived from the work of Kohn and co-authors [KSVW08] which quantified near-cloaking in terms of certain boundary measurements. Hence the main result of this paper (see Theorem 2) asserts that the difference between the solution to the cloak problem (6) and that to the homogeneous problem (5) when measured in the H 1 2 -norm on the boundary can be made as small as one wishes by fine-tuning certain regularisation parameter. To the best of our knowledge, this is the first work to consider near-cloaking strategies to address time-dependent heat conduction problem. This work supports the idea of thermal cloaking albeit with the price of having to wait for certain time to see the effect of cloaking. We also illustrate our theoretical results by some numerical simulations. We leave the study of fine properties of the thermal cloak problem for future investigations: • Behaviour of the temperature field inside the cloaked region.
• Designing certain multi-scale structures (à la reiterated homogenization) to achieve effective properties close to the characteristics of ρ cl and A cl . • Study thermal cloaking for time-harmonic sources.