We adapt tools of transformation optics, governed by a (elliptic) wave equation, to thermodynamics, governed by the (parabolic) heat equation. We apply this new concept to an invibility cloak in order to thermally protect a region (a dead core) and to a concentrator to focus heat flux in a small region. We finally propose a multilayered cloak consisting of 20 homogeneous concentric layers with a piecewise constant isotropic diffusivity working over a finite time interval (homogenization approach).
© 2012 Optical Society of America
There is currently a renewed interest in gradient index materials (whose paradigm is the Luneburg lens ) in optics and acoustics owing to their links with transformational optics and acoustics. In 2006, two research groups (those of Pendry, Schurig and Smith  and Leonhardt ) independently proposed a systematic way to control wave trajectories in curvilinear coordinate systems. In order to attract attention of mass media, they designed a cloak that renders any object inside it invisible to electromagnetic radiation. The former team theoretized that a coating consisting of a meta-material whose physical properties are deduced from a coordinate transformation in the Maxwell system displays anisotropy and heterogeneity of permittivity and permeability working as a deformation of the optical space around the object. The physicists consider the blowup of a point, thereby tearing apart the metric space. Though this may seem haphazardous, this can be legitimated by making use of advanced mathematical treatments [4, 5]. The experimental validation  of their theoretical considerations was given in 2006 for a copper cylinder invisible to an incident plane wave at 8.5 GHz. In 2009, the group of de Lustrac proposed a non-magnetic cloak, with anisotropic and spatially varying permittivity, which was experimentally demonstrated at microwave frequencies . However, Leonhardt’s team used mathematical tools of complex analysis in order to conformally map the Helmholtz equation, thereby ensuring that the (spatially varying) refractive index be isotropic. These two approaches to cloaking (conformal and non-conformal) markedly enhance our capabilities to manipulate light, even in the intense near field limit . Other research groups looked at how governing equations behave under geometric transforms in other areas of physics, such as the conductivity equation , the Schrödinger equation for matter waves [10,11,12,13], the acoustic wave equation [14, 15, 16] and the Navier equation for elastodynamics [17, 18]. One should also note that it is possible to drastically reduce the scattering cross section of an object thanks to a homogeneous coating with close to zero , or negative  refractive index.
In the present paper, we extend the design of transformation based metamaterials to the area of thermodynamics. We first show that the diffusion equation retains its form under geometric transform and we derive the expression of its heterogeneous anisotropic diffusivity in a general coordinate system. We then apply specific transforms in order to design a cloak and a concentrator for heat flux. We use full wave computations to back up our claims of an unprecedented control of heat flux through a change of metric in the thermic space. We finally propose a multilayered device working as a thermic protection in the homogenization regime.
2. Transformed heat equation
We consider the two-dimensional diffusion equation in a bounded cylindrical domain Ω with a source pEq. (1)). In this paper, we consider a source with a time step (Heaviside) variation and a singular (Dirac) spatial variation, that is: p(x,t) = p0H(t)δ(x−x0), with H the Heaviside function and Delta the Dirac distribution. This means that the source term is constant throughout time t > 0, while it is spatially localized on the line x = x0.
Upon a change of variable x = (x,y) → x′ = (x′,y′) described by a Jacobian matrix J = ∂(x′,y′)/∂(x,y), this equation takes the form:Eq. (1) and Eq. (2) have the same structure, except that the transformed conductivity Eq. (2) is to multiply Eq. (1) by a smooth function (that is an infinitely differentiable function with a compact support on Ω) and to further integrate by parts, which leads to the following variational form: 21] and references therein.
We now apply to Eq. (4) the coordinate change x = (x,y) → x′ = (x′,y′) and noting that ∇ = J−1∇′, where ∇′ is the gradient in the new coordinates, we end up withEq. (2) which lays the foundation of transformation thermodynamics. Let us now apply the transformed Eq. (2) to the design of two metamaterials of potential pratical interest: an invisibility cloak and a concentrator for heat flux, as shown in Fig. 1.
3. Transformed based thermic metamaterials
Following the work by Pendry et al. , we consider the linear geometric transform:22] in order to model domains of a (smooth) arbitrary shape.
After some elementary algebra, we obtain the following transformation matrix:22], where a similar analysis was performed for control of transverse time-harmonic electromagnetic waves propagating in cylindrical invisibility cloaks.
One should note that when r′ = β in the transformed coordinates, i.e. r = 0 in the original coordinates, the transformation matrix becomes singular as its first coefficient vanishes and the other three tend to infinity.
3.1. An invisibility cloak for heat
Let us consider two embedded smooth closed curves defined by R1(θ) < R2(θ), for 0 ≤ θ < 2π, as shown in Fig. 1 and the transformation Eq. (6) with the parameters and β = R1(θ): It maps the field within the domain r ≤ R2(θ) onto the annular domain R1(θ) ≤ r′ ≤ R2(θ). The region r′ ≤ R1(θ) defines the invisibility zone, while the annulus in the cloak consists of a material with heterogeneous anisotropic conductivity given by Eq. (7) with
We start with the computation of the Jacobian matrix of the compound transformation (x,y) → (r,θ) → (r′θ′) → (x′,y′) which can be expressed as:
We then compute the transformation matrix using Eq. (3):
We now observe that from Eq. (6)r = (r′ − R1)/α, hence the transformed conductivity inside the circular coating of the cloak can be expressed asEq. (10)
We report in Fig. 2 some finite element computations with COMSOL MULTIPHYSICS for a circular thermal invisibility cloak when we apply a constant temperature normalized to 1 on the left side of a cell and a convective flux boundary condition on the right side i.e. n · κ∇u = h(uR − u), where uR is the free temperature on the right side, u = 0 that of the ambient medium, and h = 5W.m−2.K−1 a convective coefficient. We stress that the temperature is small but does not vanish on the right side of the cell. However, the invisibility region (inner disc) displays a specific protection for heat, as the field is constant there. This comes from the fact that from Eq. (2) taken in distributional sense, we are ensured that u and are continuous across the cloak’s inner circular boundary, where the unit outward normal n to the boundary is directed along the radial axis. We thus note that the radial flux κ′r′∂u/∂r′ is continuous across this boundary. However, its trace on the outer boundary is null as from Eq. (12) κ′r′∂u/∂r′ = (r′ − R1)/r′∂u/∂r′, which vanishes when r′ = R1. This means that by continuity, the flux which is simply κ∇u inside the invisibility region (the conductivity is a constant equal to the value outside the cloak) should also vanish there. It is therefore natural to have a constant temperature inside the invisibility region. Moreover, the constant value is allowed to vary with time, which is also in agreement with the four panels in Fig. 2.
For t > 0.005s, temperature is actually not vanishing in the inner disc, see panel (c) in Fig. 2 for t = 0.02s. This is an effect due to the non-harmonic regime. Interestingly, its value increases monotonically with time until t = 0.05s. We checked that when t > 0.05s, the normalized temperature no longer increases inside the invisibility region and does not exceed T/2, with T = 1 the temperature on the left side of the cell. Moreover, the isovalues of the field showing the diffusion of temperature within the metamaterial cloak clearly demonstrate that temperature inside the coating is smoothly detoured around the invisibility region, while going unperturbed elsewhere. Such a cloak could be used in thermal protection.
We note that on the inner boundary of the cloak, i.e. when r′ = R1, the conductivity becomes infinitely anisotropic: the radial coefficient κ′r′ vanishes while the azimuthal coefficient κ′θ′ tends to infinity: this means that the heat flux is detoured around the invisibility region. Moreover, the left-hand member of Eq. (2) vanishes i.e. the transformed heat equation is thermo-static on the inner boundary of the cloak. Indeed, using again r = (r′ − R1)/α, we note that the transformed density and heat capacity in Eq. (2) are such thatEq. (2), corresponding to the blow up of the point r = 0 onto the disc r′ ≤ R2. However, the boundary condition at r′ = R1 can only be satisfied approximately by the finite element package, and one can see that temperature increases steadily inside the invisibility region until it reaches a plateau for t > 0.05s. The temperature is then constant inside the invisibility region.
3.2. A concentrator for heat flux
Let us consider three embedded smooth closed curves defined by R1(θ) < R2(θ) < R3(θ), for 0 ≤ θ < 2π, as shown in Fig. 1 and the transformation Eq. (6) with and β = 0 for 0 ≤ r ≤ R2(θ) and with and for R2(θ) ≤ r ≤ R3(θ). This transformation maps the field in the region 0 ≤ r ≤ R2(θ) onto 0 ≤ r′ ≤ R1(θ′) (i.e. compression of thermal space) and the field in the region R2(θ) ≤ r ≤ R3(θ) onto R1(θ′) ≤ r′ ≤ R3(θ′) (i.e. extension of thermal space). Importantly, the compression and extension compensate each other for 0 < r′ ≤ R3(θ′) and the transformation should be the identity from R3(θ) < r to R3(θ′) < r′. The resulting heat concentrator (in transformed variables) then consists of two parts:
We now note that if the concentrator is circular, dR1(θ)/dθ = dR2(θ)/dθ = dR3(θ)/dθ = 0 so that c12 vanishes in Eq. (14) and Eq. (15). The derivation of the T matrix in that case follows mutatis mutandis that of the previous section with (α,β) = (R1/R2, 0) if 0 ≤ r ≤ R2 and (α,β) = ((R3 − R1)/(R3 − R2), R3(R1 − R2)/(R3 − R2)) if R2 ≤ r ≤ R3.
This leads us to T−1 = R(θ′)Diag(κ′r′, κ′θ′)R(θ′)T where the cloak is described by the following parameters
Moreover, the left-hand side of Eq. (2) has the following transformed density and heat capacity:Eq. (2).
We report in Fig. 3 some finite element computations for a circular thermal concentrator when we apply a constant heat source on the left side of a cell. We note that temperature isovalues are bent towards the center of the concentrator: heat diffuses faster in the left part of the device and slower in the right part. Such a concentrator could be used to focus heat in a tiny region, what might have potential applications in photovoltaics.
4. A broadband multilayered cloak via homogenization of the heat equation
We now wish to propose a practical design of a broadband cloak consisting of a large number of homogeneous layers with piecewise constant diffusivity (ratio of the conductivity by the product of density by specific heat capacity, which is in SI units of square meter per second). For this, we need to introduce the notion of anisotropic diffusivity.
4.1. Reduced set of parameters for anisotropic diffusivity
We first note that the coordinate transformation can also compress the region r < R2 into the shell R1 < r < R2, provided that the cloak is described by the following reduced anisotropic conductivity (written in its diagonal basis):Eq. (12) by the spatially varying determinant in Eq. (13). By doing so, the determinant is now equal to 1 in Eq. (2), and the rank-2 tensor is in SI units of diffusivity (m2.s−1).
However, such a simplified set of parameters leads to an approximation of Eq. (2) as
4.2. Homogenization model for the heat equation
We still face the obstacle of an anisotropic conductivity inside the cloak. This problem can be handled using techniques of homogenization. Let us consider that the temperature field is solution of the following governing equation (outside the source)Eq. (20) in terms of a macroscopic (or slow) variable x = (r,θ) and a microscopic (or fast) variable , where η is a small positive real parameter. The homogenization of this parabolic equation can be derived by considering the ansatz Eq. (20) is not rescaled.Eq. (21). Here, is the arithmetic mean over a unit cell along the radial axis and is a homogenized rank-2 diagonal tensor which has the physical dimensions of a homogenized anisotropic conductivity given by
To mimic the spatially varying anisotropic diffusivity where is given in Eq. (18), we proceed in two steps, following : In a first step, we normalise the homogenized conductivity by <ρc>, what we call anisotropic homogeneous diffusivity . This approach enables us to consider layers of diffusivities KA = κA/(ρAcA) and KB = κB/(ρBcB), what simplifies the numerical implementation. We then approximate the ideal cloak by a multi-layered cloak with M such anisotropic homogeneous concentric layers. In a second step, we approximate each layer i, i = 1,..,M by N thin isotropic layers through the homogenization process described above. This means the overall number NM of isotropic layers can be fairly large, in our case N = 2 and M = 10. A mathematical justification for this can be found in the context of acoustic and Schrödinger equations in  by applying the homogenization technique to and , so that the homogenized coefficients in Eq. (22) are still spatially varying on the slow scale r, but this leads to mathematical technicalities beyond the scope of the present paper.
4.3. Ilustrative numerical results
We show in Fig. 4 four snapshots of heat diffusing through a multilayered cloak working via homogenization. These plots have the same features as in Fig. 2. It should be noted that the diffusivity varies over a large range of values. However, this is within an easy reach with realistic materials, as one can rescale diffusivities in caption of Fig. 2 according to, say, Polyvinyl Chloride (PVC) for the bulk material (K = 8 10−8m2.s−1), in which case silver could be the inner most layer of the cloak (K = 1.65 10−4m2.s−1). Snapshots of heat distribution at t = 0.01s (a) and t = 0.05s (b) show that normalized temperature on the left side of the cell is non vanishing inside the inner disc (invisibility region), but it never exceeds a half of the maximum value of normalized temperature (even for t > 0.05s). The fact that the temperature inside the inner disc is still a constant for a given time as in Fig. 2, whereas in Fig. 4 the cloak no longer displays singular values of the diffusivity, requires a more sophisticated explanation (i.e. it is not legitimate to invoke a thermostatic problem as in section 3.1). In order to understand why this is the case, one needs to invoke a maximum principle for parabolic equations  that lies outside the scope of the present study.
In conclusion, we have studied analytically and numerically the extension of transformation optics to the domain of thermodynamics. We have proposed to design an invisibility cloak, which reduces the temperature inside an arbitrary region, and a concentrator which focuses heat flux in a given region. We stress that the efficiency of the thermal protection with the invisibility cloak depends upon the position of its center, because the maximum value of the constant field inside the invisibility region corresponds to the value of the temperature at this point before the geometric transform is applied. We have finally proposed a design of multi-layered cloak consisting of homogeneous isotropic layers of materials with realistic diffusivities. Applications may adress the field of microelectronics, for which reason the diameters of metamaterials are 400 to 600 micrometers and the time scale considered is 20 to 50 milliseconds in our numerical examples. However the key parameter is the ratio of time to length, so that the scale can be directly modified for most general thermodynamics applications. Notice here that we assumed the materials not to be temperature dependent, but this can be directly included within our approach. Moreover, analysis of cloaking for other types of heat sources, such as instantaneous ones, or harmonic ones, is also a possible extension of this work.
Finally, one should note that there is a strong analogy between the heat Eq. (1) and the time-dependent Schrödinger equation in the context of quantum mechanics. Some time-independent quantum cloaks for matter waves have been proposed by the groups of Zhang and Greenleaf using coordinate transforms [10, 11], and space foldings have also been investigated in this context . Such results would underpin cloaking in the context of time-harmonic heat equation. However, in the present case, the transformed heat Eq. (2) is rather analogous to the time dependent transformed Schrödinger equationEq. (2) is , which is a transformed (anisotropic heterogeneous) mass, whereas the coefficient in front of the source term in Eq. (2) is now replaced by det(J)V, which is a transformed (heterogeneous) potential. Moreover, det(J) in the left-hand side is a time-scaling parameter, h̄ is the Planck constant and Ψ is the wave function. In order to facilitate the physical interpretation of Eq. (25), one can e.g. divide throughout by det(J), and consider a reduced set of parameters for the transformed mass.
S.G. is thankful for an ERC funding through ANAMORPHISM project. The authors are also grateful for stimulating discussions with M. Zerrad and for insightful comments from the anonymous reviewer on the general form of the heat equation and on useful analogies with the Schrödinger equation, which helped improving the manuscript.
References and links
1. R. K. Luneburg, Mathematical Theory of Optics (University of California Press, 1964).
4. A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, “Full-wave invisibility of active devices at all frequencies,” Commun. Math. Phys. 275(3) 749–789 (2007). [CrossRef]
5. R. V. Kohn, H. Shen, M. S. Vogelius, and M. I. Weinstein, “Cloaking via change of variables in electric impedance tomography,” Inverse Probl. 24015016 (2008). [CrossRef]
6. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314, 977–980 (2006). [CrossRef] [PubMed]
7. B. Kanté, D. Germain, and A. de Lustrac, “Experimental demonstration of a nonmagnetic metamaterial cloak at microwave frequencies,” Phys. Rev. B 80, 201104(R) (2009). [CrossRef]
9. A. Greenleaf, M. Lassas, and G. Uhlmann, “On nonuniqueness for Calderon’s inverse problem,” Math. Res. Lett. 10, 685–693 (2003).
10. A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, “Isotropic transformation optics: approximate acoustic and quantum cloaking,” New J. Phys. 10, 115024 (2008). [CrossRef]
12. A. Diatta and S. Guenneau, “Non singular cloaks allow mimesis,” J. Opt. 13, 024012 (2011). [CrossRef]
13. A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, “Schrödinger’s Hat: Electromagnetic, acoustic and quantum amplifiers via transformation optics,” (preprint:arXiv:1107.4685v1).
14. S. A. Cummer and D. Schurig, “One path to acoustic cloaking,” New J. Phys. 9, 45–45 (2007). [CrossRef]
15. H. Chen and C. T. Chan, “Acoustic cloaking in three dimensions using acoustic metamaterials,” Appl. Phys. Lett. 91, 183518 (2007).
16. A. Norris, “Acoustic cloaking theory,” Proc. R. Soc. London 464, 2411–2434 (2008). [CrossRef]
17. G. W. Milton, M. Briane, and J. R. Willis, “On cloaking for elasticity and physical equations with a transformation invariant form,” New J. Phys. 8, 248–268 (2006). [CrossRef]
18. M. Brun, S. Guenneau, and A.B. Movchan, “Achieving control of in-plane elastic waves,” Appl. Phys. Lett. 94, 061903 (2009). [CrossRef]
19. A. Alu and N. Engheta, “Achieving transparency with plasmonic and metamaterial coatings,” Phys. Rev. E 95, 016623 (2005). [CrossRef]
20. G. W. Milton and N. A. Nicorovici, “On the cloaking effects associated with anomalous localised resonance,” Proc. R. Soc. London, Ser. A462,3027–3059 (2006). [CrossRef]
21. V. V. Jikov, S. M. Kozlov, and O. A. Oleinik, Homogenization of Differential Operators and Integral Functionals (Springer-Verlag, New-York, 1994). [CrossRef]
22. A. Nicolet, F. Zolla, and S. Guenneau, “Electromagnetic analysis of cylindrical cloaks of arbitrary cross-section,” Opt. Letters 33, 1584–1586 (2008). [CrossRef]
24. M. Rahm, D. Schurig, D. A. Roberts, S. A. Cummer, D. R. Smith, and J. B. Pendry, “Design of electromagnetic cloaks and concentrators using form-invariant coordinate transformations of Maxwell’s equations,” Photonics Nanostruct. Fundam. Appl. 6, 87–95 (2008). [CrossRef]
25. L. Nirenberg, “A strong maximum principle for parabolic equations,” Commun. Pure Appl. Math. 6, 167–177 (1953). [CrossRef]