Analogies between optical propagation and heat diffusion: applications to microcavities, gratings and cloaks

A new analogy between optical propagation and heat diffusion in heterogeneous anisotropic media has been proposed recently by three of the present authors. A detailed derivation of this unconventional correspondence is presented and developed. In time harmonic regime, all thermal parameters are related to optical ones in artificial metallic media, thus making possible to use numerical codes developed for optics. Then, the optical admittance formalism is extended to heat conduction in multilayered structures. The concepts of planar microcavities, diffraction gratings and planar transformation optics for heat conduction are addressed. Results and limitations of the analogy are emphasized.

A novel association of optics and heat has been proposed recently [16][17][18][19][20][21][22][23][24][25][26][27][28]. It is based on an analogy between the governing equations of propagation for optics and diffusion for heat. The correspondence is not intuitive according to the markedly different nature of optical propagation and thermal diffusion phenomena. Nevertheless, an analogy can be unveiled in time harmonic regime when metallic media are considered in optics. Under these conditions, electromagnetic fields decay exponentially fast away from a metal-dielectric interface, so that their behaviour is somewhat reminiscent of a temperature field undergoing a fast decay away from a heat source. Such an analogy has been first used in the design of thermal cloaks, heat concentrators and rotators [16][17][18][19][20][21][22][23][24][25][26][27][28], using geometric transformation methods previously applied to electromagnetics [29,30]. A series of papers subsequently emerged and confirmed tools of transformation optics can be applied to control other diffusion processes [31,32].
In this work, we further explore the unconventional analogy between transformation optics and thermodynamics which has been first proposed in reference [16]. Correspondences are emphasized between optical and thermal fields, i.e. (vector) electric fields and (scalar) temperature, (vector) magnetic fields and (vector) heat flux. All thermal parameters are related to optical parameters in artificial metallic media, thus making it possible to reintroduce the optical concepts of effective index and complex admittance for heat conduction. More generally, all numerical codes developed for optics (multilayers, gratings, microcavities, scattering and diffraction, transformation optics . . .) can be used. The plan of this article is as follows: the whole optical admittance formalism [33][34][35][36][37][38] is first extended to thermal conduction, and deep analogies are drawn with multilayered electromagnetic structures in order to predict thermal properties in multilayers. Next, the concepts of planar microcavities, diffraction gratings and planar cloaks are addressed in the frame of heat conduction. A particular attention is paid to the mathematical correspondences and limitations between metal optics and heat conduction, including fields and energy balances, together with the extended optical admittance formalism for conduction.

Optical propagation versus heat conduction (a) Spatio-temporal regime
Free space optical propagation in linear, isotropic and homogeneous materials can be classically modelled using the governing equation [38][39][40][41][42][43]: with E the vector electric field, ε and μ the temporal permittivity and permeability, J and q the density of currents and charges and S opt a source term resulting from these densities. In equation (2.1), the convolution product (* t ) involves time (t) and takes into account the material inertia, that is, the time delay between the excitation of material and its optical response. Note the second-order derivation versus time of the electric field (elliptic equation).
On the other hand, heat conduction [43][44][45][46] involves a first-order derivation versus time (parabolic equation). We will here consider that conduction satisfies the following equation in linear, isotropic and homogeneous media as with T the temperature and a, b the thermal parameters, that is, diffusivity and conductivity, respectively. The source term S is the bulk density of heat power. Note that in equation ( In addition to the discrepancy in time derivation order, other straightforward differences can be emphasized between optical propagation and heat diffusion, which are -heat equation is scalar, whereas optics equation is vectorial; however, one can consider optics equations on each basis vector, and -optical parameters (ε, μ) vary with time and hence will exhibit a frequency dispersion, whereas this is not the case for the thermal parameters in equation (2.2).
Another difference lies in the fact that the temperature (or a flux condition) is most often fixed at the domain frontier, a supplementary condition not present in free space optics. However, in the absence of sources, heat flux and temperature are continuous at interfaces, as it is for the tangential electromagnetic fields. In equations (2.1) and (2.2), all fields vary with time (t) and space location ρ = (r, z) = (x, y, z). It is interesting to recall the Green's functions G of equations (2.1) and (2.2) that give the particular solutions in the form G * t S opt or G * t S th , with H(t) in conduction, (2.4) with ρ = |ρ|, H the Heaviside function and δ the Dirac distribution. Equation (2.3) is given under the assumption that the optical material is not dispersive, and shows that the optical wavefront propagates at velocity v = 1/ √ (εμ) and is located at time t at the sphere surface of radius ρ = vt where the whole energy is located. On the other hand, equation (2.4) is for a diffusion process, so that the energy is spread within the whole heated volume of radius ρ = 2 √ (at) which increases with time. Note that (2.4) does not involve any Dirac distribution: strictly speaking, heat would have already diffused everywhere at arbitrary time (t = 0 ⇒ G th = 0 ∀ρ), owing to the Gaussian nature of G th .
In what follows, we now consider the propagation or conduction equations in isotropic homogeneous regions free of sources, what we write as the homogeneous equations Moreover, we will be interested in a temperature elevation resulting from a transient heat source, similar to what happens when the heat source results from the absorption of a modulated laser beam. Following the superposition properties of linear systems, the total temperature satisfies T = T a + T with T a owing to the environment, and T resulting from the heat source given as with |h(t)| < 1. This last temperature T can be written as T = T 1 + T 2 , where T 1 is a steady-state temperature and results from S 0 , whereas T 2 is a transient temperature and results from S 0 h(t). In this work, we are actually interested in this last transient temperature T 2 .

(b) Harmonic regime
For the sake of simplicity, we use Fourier transform rather than Laplace transform, despite the differences in the two transforms. Another reason is that a modulated laser beam will create optical absorption ( with k the wavenumber (or conduction number by analogy with optics), u the scalar electric field or temperature and u its Fourier transform defined aš shows that both fields (u = E or T) satisfy a similar harmonic equation in the Fourier space, taking into account the specific dispersion law of the k function, that is and In other words, the dispersion law k(ω) is the memory of the pth order of time derivation (p = 1 for conduction, p = 2 for optics). Note in (2.9) thatε(ω) andμ(ω) are the Fourier transforms of the optical functions (ε, μ), whereas in (2.10), the diffusivity parameter a is not dispersive under the thermal approximation. Hence, the k dispersion behaves as √ ω for conduction, whereas in optics, it depends on the permeability and permittivity frequency dispersion law; however, as a first approximation, the k dispersion in optics can be considered to follow ω in a narrow bandwidth.
Above all, a major point to be emphasized in the harmonic regime concerns the k value at a given frequency. Indeed, in optics, such a value can be real or complex, depending upon whether dielectrics or metals are under study. On the other hand, the k value is necessarily complex in conduction for reasons following (2.7)-(2.10), so that heat diffusion can be considered to be analogous to optical propagation (more exactly, electric field attenuation) in artificial metallic media. Furthermore, following (2.10), the equivalent metal should be specific with identical real and imaginary parts in the refractive index n (defined as k = k 0 n in optics).

Optical admittance formalism for multilayers (a) Multilayer optics: basics
In multilayer optics (figure 1), the admittance formalism is commonly used to predict optical properties and to design antireflective coatings, beam splitters, polarizing devices, mirrors and narrow-band filters [33][34][35][36][37][38]. It is based on the complex admittance, a key function that often simplifies the calculation and allows for analytical design. This admittance function Y is given for a particular polarization (transverse electric S or transverse magnetic P) aŝ withĤ andÊ the tangential components of the electric and magnetic polarized field, and z the direction normal to the stack. Note here that all fields are harmonic (monochromatic) and result from a single illumination incidence at the top surface of the stack; in other words, the admittance is defined after double Fourier transform (versus time t and versus transverse coordinate r) of the fields; we will come back to this point in the next section, where we re-derive the admittance formalism. In the case of a stationary wave within the stack, the admittance varies with the z altitude; on the other hand, when the wave is only progressive or retrograde, what may happen in free space in the surrounding media, it is constant and reduced to the effective index ±m, that iŝ Equations (3.1) and (3.2) are given for each polarization of the optical field, which means that both effective indices and admittances are polarization (and incidence, wavelength) dependent. Moreover, we have assumed to be in the far field, so that the electric field E is driven by two transverse polarization modes. The effective indices are known from the illumination conditions, and enable to initiate a sequence for the admittance extraction throughout the whole stack (see next section). They are given as with η 0 = √ (μ 0 /ε 0 ) the vacuum impedance, n = √ (ε rμr ) the refractive index, (ε r ,μ r ) the complex relative permittivity and permeability, α = √ (k 2 − σ 2 ) and σ the spatial pulsation which can be connected in transparent media to a propagation direction θ (illumination angle), that is with ν the spatial frequency; in the case of plane waves (σ < k), relation (3.3) is turned into: m = (1/η 0 μ r ) n cos θ for TE polarization and m = (1/η 0 μ r ) n/ cos θ for TM polarization. All relations (3.1)-(3.4) are valid for plane, evanescent or dissociated waves. the next sections), we start with a problem invariant with respect to the transverse coordinates (x, y), so that relation (2.7) only involves z and ω variables ∂ 2ǔ ∂z 2 + k 2 (ω)ǔ = 0. (3.5) In optics, the unique z-dependence corresponds to the case of a plane wave illuminating the stack at normal incidence (σ = 0). The solution of (3.5) is given in each layer (i = 1, N) by

(b) Conduction in multilayers
where the field is E or T in the layer. The next step in optics to complete the admittance formalism is based on the fact that the magnetic fieldv follows the same analytical expression (3.6), and that the electric and magnetic components are linked through (3.1) and The last step consists of combining (3.8)-(3.9) and writing the continuity of the tangential electric and magnetic fields in the absence of sources, which results in a matrix formalism to pass from one interface to the next as and with e i the thickness of layer i. Now, we can check that the similarity with the conduction process is immediate. First of all, we need two quantities that are continuous throughout the whole stack; in the absence of singular sources, these quantities are the harmonic temperatureŤ and the scalar (z-projection) heat flux densityȟ given byȟ Hence, the temperature plays the role of the electric field in (3.6), whereas the heat flux replaces the magnetic field in (3.7). However, this last condition must be further clarified to determine which quantity replaces the effective index m in (3.8). By analogy with optics (relation (3.2)a,b)), this last quantity is defined for conduction aš The result is immediate from (3.7), that is 14) and the admittance follows within the stack aš At this stage, the correspondence is complete and the matrix formalism (3.10) and (3.11) can be used either in optics or conduction in the harmonic regime. Note that two slight modifications were necessary to reach this one-to-one correspondence, which are is replaced by the conduction parameter (c) Wavelength, diffusion length and phase velocity From equations (2.8) and (3.6), the field (T or E) can be reconstructed in the layers as the sum of elementary components. In the substrate, these components take the form (3. 16) withǔ + = |ǔ + | exp(jψ + ) and k = k + j k . The real exponential is for the field envelope, whereas the cosine term allows one to retrieve a spatial period λ = 2π/k (the wavelength in optics) and a phase velocity v = ω/k (light velocity in optics). Equation (3.16) shows that in conduction the (pseudo)spatial period λ th is proportional to the diffusion length is analogous to the phase velocity of optics. However, we do not extend the analogy to the refractive index, because this index involves the vacuum property in optics (where propagation occurs), whereas conduction fails in vacuum (another reference material should be used).
(d) Two-dimensional formalism for spatial wave packet Now, we take into account the spatial distribution of the field, which means that we drop the unique z-dependence to recover a ρ = (r, z) variable, with r = (x, y). Such an extension to a threedimensional geometry is immediate if we consider, in addition to the frequency wave packet (Fourier transform versus time), a spatial wave packet (Fourier transform versus space coordinate r). This second Fourier transform is given aŝ with ν the spatial frequency, the Fourier variable conjugated of r. Such double Fourier transform u is governed by the same equation (3.5) as the single Fourier transformǔ, except that k(ω) must be replaced by α(ω) = √ (k 2 − σ 2 ), with σ = 2π ν the spatial pulsation. We obtain In each layer, the field in this second Fourier plane takes the form so that the elementary components that enable to rebuild the original (spatio-temporal) field can once again be written in the substrate aŝ withû + = |û + | exp (jφ + ) and α = α + jα . In other words, this works like optics at oblique incidence and the k parameter only has to be replaced by α = √ (k 2 − σ 2 ). The salient consequence for the matrix formalism (3.10) and (3.11) following the previous section is that the effective index for heat conduction (relation (3.14) must also be modified as m = −jαb. (3.23) Regarding the spatial and temporal period, k must be replaced by real(α). At last, the original (real) field is reconstructed with a double inverse Fourier transform over temporal (f ) and spatial (ν) frequency.

Energy balance
To be complete, our correspondences should also involve an energy balance for both optical propagation and heat conduction. This condition is written below within a domain Ω limited by a closed surface Σ with unit outward normal n (figure 1).

(a) Temporal regime
In optics, we obtain from Maxwell's equations in the spatio-temporal regime with P = E ∧ H the real Poynting vector, F the optical power provided by sources within Ω, Φ the Poynting flux through the surface Σ, and A the variation of internal energy (or absorption) within Ω, with dA/dΩ = ∂w/∂t = H.∂/∂t(μ * H) + E.∂/∂t(ε * E). As usual, the density of electromagnetic energy reduces to w = (1/2)[E 2 + H 2 ] when the material inertia is neglected. For purely heat conduction, the balance is directly obtained from the governing equation integrated over Ω and with F the heat power provided by the source, Φ the heat flux through Σ and A the variation of internal energy, with: In other words, the density of electromagnetic energy is here replaced by a quantity proportional to the temperature, with b/a = γ C the product of mass bulk density (γ ) and heat capacity (C).

(b) Fourier analysis
Such balances (4.1) and (4.2) must now be rewritten in the Fourier planes. In optics, quadratic (intensity) detectors are known to realize a time average process of the square fields, owing to their low temporal response with respect to optical periods. Such specificity allows one to transform the spatio-temporal balance (4.1) into a single harmonic balance at each temporal pulsation ω (first Fourier plane), that is with Π the complex harmonic Poynting vector: Then, to reach a balance at a single spatial frequency ν in the second Fourier plane, one may use Parseval's theorem and other specific properties to find − ν,zĴ In a last step, considering the hermiticity of the fields, we recover the harmonic balance per unit of spatial frequency dν x dν y with respect to a surface Σ surrounding the whole planar multilayer with sides at infinity (figure 1) withε = Im(ε),μ = Im (μ) and the notation [u] z = u(z 2 ) − u(z 1 ), where z i are the upper and lower altitudes of Σ and z 2 > z 1 . Note here that we used equations (3.1) and (4.6) to reach the property: 2Π · n = Y|Ê| 2 , a flux quantity which is continuous throughout the stack in the absence of sources.
Addressing the harmonic balance in conduction is a priori simpler, because temporal frequencies are much lower than in optics, thus disregarding temporal integration by the detector is legitimate. Moreover, here the thermal parameters are assumed to be non-dispersive, which simplifies the calculation. One can directly use the Fourier transform versus time of the governing equation (4.2) to reach a complex harmonic balance at temporal pulsation ω Per unit of surface area dxdy parallel to the interfaces of the stack, we obtain Owing to the Hermitian property ofŤ such a balance can also be limited to positive frequencies ω as follows Now, to write the conduction balance at a single spatial frequency (second Fourier plane), a z integration after a double Fourier transform of (4.2) directly leads to with α 2 = k 2 − σ 2 . Note again in (4.11) that YT is a flux quantity which is continuous within the stack.

Optical admittance calculation for heat conduction
This section is devoted to numerical calculation. For that, we use an optical thin film admittance software where the k dispersion and the effective index were rewritten following (2.10) and (3.23).
As usual [33][34][35][36][37], the effective index value in the substrate is the starting value of the admittance sequence within the stack, given at each interface by A similar sequence is given for the field at interfaceŝ A complete z-variation of admittance and field can also be obtained in the layers if we replace δ i = α i e i by δ i = α i z, with 0 < z < e i . Finally, reflection (r) and transmission (t) factors can be defined in a way similar to optics, that is at the bottom interface In optics, the flux balance related to a far field surface surrounding the whole multilayer leads to with Φ the Poynting fluxes which are incident (Φ + 0 ), reflected (Φ − 0 ) and transmitted (Φ + S ), and A the normalized absorption. For an elementary plane wave, these fluxes are proportional to Φ = real(m)|ũ ± | 2 , so that the flux balance can be rewritten for the surface Σ as The heat balance is slightly different in the sense that, owing to heat diffusion, we cannot consider a source located at infinity (z = −∞), so that the far field position of the source is not valid and may underpin the separation between incident and reflected fluxes. However, the top interface temperature T 0 is assumed to be forced (and known), so that the whole field distribution can be obtained within the stack, from the sequence (5.2). Then, the spectral quantity which we will analyse here is the ratio t th =T + S /T 0 = t/(1 + r), withT + S the field at the substrate interface. Hence, in the figures of this section, all fields are normalized toT 0. Now, some elementary or preliminary remarks are useful before numerical calculations. Indeed, numerous tools were developed in optics to design optical coatings and address specific challenges, which can be briefly addressed here in conjunction with heat conduction.

(a) Preliminary remarks (i) Single boundary
Consider first the case of a single boundary. As in optics, reflection is given by For a unique z-dependence (analogous to normal incidence), we obtain with β i the effusivity This last equation shows that, when the calculation of reflection is considered, the effusivities in heat conduction play the role of the refractive indices in optics. This is valid at normal incidence, that is, for a zero spatial frequency (σ = 0). It is interesting to note that the range of effusivities is much larger in conduction (more than two decades) than in optics (half a decade). Besides from that, temperature reflection from a single boundary is a real quantity. We also obtain 1 + r = t ⇒ t th = 1.

(ii) Specific multilayers
Absentee (half-wave) layers and quarter-wave stacks play a key role in multilayer optics [33][34][35][36][37][38], owing to their stationary properties. For these stacks, the phase term δ in the matrix (3.11) must be real and a multiple of π /2, which strongly simplifies the matrix and allows for an analytical design. However, this is only valid for quasi-transparent materials (δ is a real number), whereas the analogy between optics and conduction is for metallic layers. Therefore, the concept of quarter-wave stacks does not hold for conduction. We could also wonder whether other concepts in metal optics could still be used for conduction. For instance, induced transmission filters [33][34][35][36][37] combine dielectric and metallic films to produce narrow-band filters with large rejection bands; with such devices, the optical field can be confined within the spacer (cavity) layer in a way similar to alldielectric Fabry-Perot filters. Unfortunately, this will not be allowed for conduction multilayers: indeed, induced transmission filters require mixing metallic and dielectric layers, whereas for the optics/conduction analogy, all materials should be metallic with identical real and imaginary indices.

(iii) Field extrema
Another question concerns the ideas of field confinement and enhancement [47][48][49], which are major challenges in optics. Indeed, in optics, the harmonic field can spatially oscillate in the stack and specific coatings can be designed to reach huge enhancement of excitation in particular layers. In conduction, this would involve a local maximum of temperature within the stack, which is not possible in the absence of sources according to thermodynamics.
This result can be directly recovered from the governing equation (3.20) in the second Fourier plane, which can be turned into bû * ∂ 2û ∂z 2 + bα 2 (σ , ω)|û| 2 = 0, which, after z integration, leads in each medium to where the first term within the brackets is continuous. Then, the presence of a maximum would allow one to choose these brackets between the z 1 location of the extremum and z 2 rejected at infinity where the field is zero. The result would be z b ∂û ∂z 2 = z bα 2 (σ , w)|û| 2 (5.11) which cannot be satisfied either for real or imaginary parts, because α 2 = jω/a − σ 2 . Such a result (no temperature extrema in the absence of sources) does not contradict the fact that thin metallic layers create interferences in optics and allow the field to oscillate with depth location; indeed, one can check that optical oscillations vanish when the real and imaginary indices of the metal are identical, which underpins the analogy with conduction. For similar reasons, one can show that optical phenomena such as total internal reflection and total absorption [47][48][49] are not achieved in conduction. Then, using the optical admittance software, we investigated the stack thermal properties versus z location and temporal frequency (ω), at a given spatial frequency σ = 0. Figure 2 is given for the spectral transmittance t th (ω) and emphasizes a classical behaviour versus increasing frequencies. The frequency range is [10 −2 , 10 Hz].   The temperature modulus |û i | is defined in the second Fourier plane from relation (3.21). Because a single temporal (ω) and spatial frequency (ν) are here involved, it gives in degree Kelvin the amplitude of the time variations of the one-dimensional harmonic temperature u i (z, t) as followŝ From a practical way, the static temperature should be superimposed to u i (z, t). Furthermore, as mentioned in the introduction, the temperatureû i is normalized to its value at the top interface where it is forced, so that knowledge of the source is not required. Note that forcing the temperature can be reached with surface or bulk sources, as discussed in §6. It was also necessary to check the energy balance of conduction given in (4.11). To do that, we chose a region free of sources, that is, the stack delimitated by the top interface (in contact with the superstrate) and the bottom interface (in contact with the substrate). Within this domain, relation (4.11) becomesŜ (5.12) which means that because no heat power is provided within the domain, the bulk temperature elevation is responsible for the surface flux variation between the top and the bottom surfaces of the stack. In figure 6, both real and imaginary parts of these two terms (flux and temperature integral) are plotted and the results show a perfect agreement. We also checked that the agreement holds regardless of the spatial frequency σ . Note again that the temperature was normalized to its value at the top interface, and that the flux modulus is given in W m −2 .

(c) From metallic to transparent optical multilayers
To complete this section, figure 7a − d illustrates how the admittance diagrams and the field distribution of metal optics (which we used for conduction) approach those of dielectric multilayers (where the analogy with conduction breaks down) when the imaginary index is decreased. To achieve this correspondence between admittance diagrams in heat conduction and light propagation, we considered a thermal multilayer with conduction number k = k + jk , and plotted the admittance diagrams at Ω These panels confirm that specific admittance circles are progressively recovered when the media become transparent, which is a well-known fact in thin film optics [33]. Also, the field oscillations within the stack are recovered when the imaginary index vanishes, in agreement with (5.11).

Microcavities
Optical formalisms developed for luminescent microcavities also remain valid for heat conduction. Until now, multilayers were excited with a source in free space in the superstrate. By opposition with this far field geometry, the same devices are known as microcavities when they support the sources in their bulks. The optical (linear) solution [48,49] remains strongly similar whatever the position of the source (far or near-field), provided that the field discontinuities are corrected, that is, z ∧ δE = z ∧ δH = 0 for optical coatings (far field), and z ∧ δE = 0, z ∧ δH = J for optical microcavities (near-field) with J a surface electric current. Extension of the admittance formalism [48,49] then gives the solution in the second Fourier plane aŝ whereû i andû i are the fields on each side of interface i in media i − 1 and i, respectively, andĴ i the electric current density at surface i which creates these fields. The admittances characterize each half part of the stack; they are still calculated from (5.1), but their initial values of effective index are taken in the substrate for Y and in the superstrate for Y [48,49]. The field sequence (5.2) can then be used to obtain the field valuesû ± 0 in the extreme media. After summation over all currents, the result isû where C ± i are optical factors which are design dependent, and that allow to reach the superstrate (C − i ) or the substrate (C + i ). Consider now the case of thermal microcavities, in the form of planar multilayers supporting thermal sources at theirs interfaces. The conduction process of such devices will again behave like optical cavities in metals, because the discontinuities for temperature and flux are again given in a way similar to optics, that is, for a surface source S i at interface i This last relationship is identical to (6.1) and allows one to address thermal microcavities with the same optical softwares, provided that the scalar thermal sources (S i ) replace the algebraic electric currents J i . Contrary to the far field geometry, the thermal field now can oscillate within the stack, owing to the presence of thermal sources. In the case of bulk sources, the field would follow the source at high frequencies. Note, however, that optics involves the intensity (the field square), which introduces additional interaction terms (Ĵ iĴ * j ) between the currents [48,49]. To conclude this section, we note that the power (optical or thermal) provided by the confined source within the stack depends on the multilayer geometry [48,49].

Diffraction gratings
Optical diffraction by a periodic surface is well known to spread the incident energy into a discrete series of spatial frequencies given by with d the grating period, σ 0 the incident spatial pulsation and q a relative integer. In terms of scattering angles θ p in air, the result is with λ the wavelength and i 0 the illumination angle. Several techniques exist to predict the angular efficiency at each diffraction order m, like modal and differential methods, finite-elements and boundary integrals. The same computer codes can again be applied for diffraction of heat conduction, a process that makes temperature only significant at specific frequencies or directions given by (7.1) and (7.2). However, unlike for optics, one must take account of the attenuation resulting from thermal diffusion, which confines heat diffraction in the vicinity of the grating. Moreover, because diffraction does not occur for subwavelength gratings, the grating period at normal illumination (excitation) would be greater than one wavelength (d > λ); such a condition is seemingly trivial in optics but not in conduction, where wavelength and thermal length are similar (λ = 2π L). Indeed, the consequence is that the conduction diffraction process will be significant at a distance z of the average surface lower than its period, that is, z < L < d/2π . Finally, depending on the thermal length value L, far field effects can be masked or not, which forces diffraction to be emphasized in the near-field.
Numerical calculation is given in figures 8 and 9, at a thermal pulsation of 10 Hz. The grating is a surface with a sinusoidal shape graved at the top of a silver (Ag) substrate, with an amplitude of 1 cm and a period of 10 cm. The superstrate is air. The heat source results from the absorption of an optical beam of 37.5 cm in size, at normal incidence and 633 nm wavelength. The complex refractive index of Ag at this wavelength is chosen as n Ag = 0.134 + 3.99i. Owing to the low depth penetration of optics within the metal, such heat source can be assumed to be located at the grating surface. Moreover, the thermal parameters of Ag are 1.71 × 10 −4 m 2 s −1 (diffusivity) and 418W (mK) −1 (conductivity).
The boundary integral formalism for the electromagnetic wave scattering from onedimensional rough surfaces [50,51] has been adapted to this heat conduction problem. The temperature, which is set continuous across the grating surface, and the normal component of the heat flux, which discontinuity identifies with the heat source, are the unknown functions of a set of coupled boundary integral equations. These equations are numerically solved with the method of moments [50,51] using piecewise-constant basis functions and point-matching testing functions. The discretized grating surface is 3 m long, with a sampling step of 1.5 mm. Figure 8 gives the two-dimensional modulus of the temperature distribution in the substrate and superstrate. As expected, the field shows periodic variations in both media near the centre of the beam footprint. However, the harmonics which are specific of the grating law cannot be  figure 8, due to the fact that calculation is for the near-field, which requires to superimpose all harmonicsǔ (x, z, ω) = qǔ q (x, z, ω). (7.3) To go further, we have plotted in figure 9 the modulus of the different grating harmonicsǔ q in the superstrate and in the substrate. The results emphasize the validity of the grating law (7.2) in conduction.
In figures 9 and 10, the incident optical power is chosen to 1 W. Because the silver absorption is around 3%, the absolute absorption responsible for the heat source is around 30 mW. From a practical point of view, the harmonic temperature could be increased by a factor 30 when using an absorbing black coating (99% absorption), provided no damage threshold occurs owing to the total temperature. Note that beam focusing would not help, because several periods of the grating must be illuminated. Furthermore, near-field probing would make easier the detection of temperature variations in the grating vicinity, where the thermal attenuation is weak. The grating period can also be increased to allow classical measurements.
Other diffraction processes can be performed in a similar way, including scattering from inhomogeneities (roughness or bulk) and Mie theory for spherical particles. The diffraction problem for heat conduction could also be addressed within the formalism of microcavities. Indeed, a surface thermal source can be created by the deposition of a very thin (less than λ/50) metallic layer on a dielectric material, followed by an etching process with a period d. Under optical illumination, optical absorption in this device will only be concerned with the metallic patterns. When the beam is modulated at frequency ω, the resulting heat source is with δ the Dirac distribution, h the shape of isolated pattern and A an absorption term. Then, Fourier transformation allows one to retrieve the grating law in the form and heat diffraction can be directly calculated following (6.3).

Towards a multilayer planar cloak
The last example in the analogy concerns transformation optics, a tool currently used for cloaking in different fields including heat [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]52,53]. Invisibility cloaks usually have circular or spherical geometries, while we keep here the planar geometry of previous sections and consider a single layer that should be 'hidden' with the help of additional symmetrical stacks on its left and right sides. In order to design such a planar cloak, we apply a geometric transform in the spirit of Kohn et al. [53] [29]. While the two transforms make sense in two-dimensional and three-dimensional cases, the former one (Kohn's transform) seems to be more natural in our one-dimensional configuration. With a linear transformation (x = αx + β, y = y), we obtain an ideal anisotropic planar cloak surrounding a central (spacer) anisotropic layer. All media (cloak and spacer layer) are anisotropic following the transformation. The conductivity matrix (b ij ) of the cloak is diagonal with b 11 = (e 1 − e 2 )/(e 1 − η) and b 22 = (e 1 − η)/(e 1 − e 2 ), whereas that (b ij ) of the spacer layer is also diagonal with b 11 = e 2 /η and b 22 = η/e 2 . In order to achieve these required anisotropic media, we consider an alternation of layers with isotropic conductivities which behave like an effective medium with anisotropic conductivity, see references [16,18] for construction of circular counterparts of the planar cloak with homogenization techniques. The result for the cloak is an alternation of isotropic However, although this homogenization process works for the cloak, this is not the case for the spacer layer. Because the spacer conductivity matrix involves two real diagonal terms that are real infinite and zero, respectively, the transparency effect with an isotropic spacer layer can only be reached for very large conductivities inside the spacer layer. This prediction is confirmed in figure 10, because the temperature modulus (similar to that of figures 3-5) outside the cloak is nearly superimposed to that of the initial homogeneous medium. In figure 10, the spacer layer conductivity is far greater (several decades) than the other conductivities (i.e. 10 5 higher than that of the initial medium). Complete geometry of the multilayer stack is given in figure 11. These results confirm that transformation optics still remain valid in the case of heat diffusion. Clearly, here, cloaking is limited to some very specific limiting case, owing to the high conductivity required within the spacer layer, but the mathematical technique remains the same.

Static versus harmonic regimes
The matrix formalism (3.10) developed above works whether the regime is static (Ω = 0) or dynamic (Ω = 0). In the static regime, the matrix coefficients take a single form, and the resulting product matrix for the multilayer becomes M = 1 e t /b eq 0 1 (9.1) with e t the total thickness and b eq the equivalent conductivity e t b eq = k e k b k , (9.2) which is a well-known result. However, there is a key difference with the harmonic regime which concerns the admittances. Indeed, at non-zero frequencies, the admittances only depend on the multilayer design (geometry, optical or thermal parameters), whereas in the static regime, they also depend on the extreme fields in the form This point recalls why the frontier temperatures (or fluxes) must be forced in the static regime, which is not required in the harmonic regime. This is related to the fact that the one-dimensional static solution in the extreme media free of sources involves two parameters (T = αz + β), whereas the harmonic solution requires one parameter (T = T s exp (jαz)); in other words, with the harmonic regime the derivative of T is proportional to T, which is not the case for the static solution.
The analogy with optics still holds in the static regime, because electrostatics also involves a potential difference (rather than a potential), which is classical. However, in the harmonic regime (electromagnetism), there is no supplementary condition at the domain frontiers, so that the multilayer problem is entirely solved by the sequence relationship (5.1) of the admittances, which allows one to retrieve the top admittance Y 0 (m s ) from the substrate one. However, in some situations like modal optics (guided waves), another condition must be fulfilled to force the field to only merge both in the superstrate and substrate, that is Y 0 (m s ) = −m 0 . The result is a discretization of (5.1) which gives the guided modes as the poles of the reflection factor. A similar discretization would be obtained in the harmonic regime when temperature (or flux) is forced in both extreme media.

Conclusion
We have explored an analogy between optical propagation and heat conduction in isotropic, homogeneous and linear media. Both fields (scalar electric field and temperature) were shown to follow the same harmonic equation in the second Fourier plane, provided that the optical wavenumber k opt = ω √ [ε(ω)μ(ω)] is replaced by the heat conduction number k th = (1 + j) √ [ω/2a]. The dispersion law of this wavenumber has the memory of the time derivation order of the governing spatio-temporal equation. Moreover, the analytic expression of the fields, conduction and wavenumbers show that thermal diffusion is strongly analogous to optical propagation in metallic media; however, the optical refractive index of these artificial metallic media should have identical real and imaginary parts, a specific property where applications of metal optics can be reduced.
The case of planar multilayers was addressed, and the optical admittance formalism was shown to be directly applicable to thermal conduction. Within this framework, the temperature and algebraic heat flux played the role of electric and magnetic field tangential components, respectively. The only notable difference in the admittance formalism is the required modification of the optical effective index in order to take into account the heat conduction number (m th = −jαb replaces m opt = nα/k). Numerical calculations were presented to validate all results, including the optical and thermal energy balances.
The analogy between metal optics and conduction was then extended to microcavities and diffraction gratings, and to transformation optics. All results show that most softwares developed for metal optics can be directly used for conduction, whether they address interference filters or microcavities, diffraction gratings or photonic crystals, metamaterials or transformation optics. We finally note that Green's functions for reflected and thermal fields at planar interfaces were computed in reference [52] and they are reminiscent of what researchers compute in optics. However, one must keep in mind that these analogies only hold for specific metals (real index = imaginary index). We hope our work will foster theoretical and experimental efforts in thermal metamaterials.
Data accessibility. All data and calculation methods are given in the manuscript, so that all results are fully reproducible.
Authors' contributions. All authors gave substantial contribution to theory, conception and design, numerical calculation, analysis and interpretation of data, drafting and revising the article and final approval of the version to be published.