Layered thermal metamaterials for the directing and harvesting of conductive heat

The utility of a metamaterial, assembled from two layers of nominally isotropic materials, for thermal energy re-orientation and harvesting is examined. A study of the underlying phenomena related to heat flux manipulation, exploiting the anisotropy of the thermal conductivity tensor, is a focus. The notion of the assembled metamaterial as an e ff ective thermal medium forms the basis for many of these investigations and will be probed. An overarching aim is to implement in such thermal metamaterials, functionalities well known from light optics, such as reflection and refraction, which in turn may yield insights on e ffi cient thermal lensing. Consequently, the harness and dissipation of heat, which are for example, of much importance in energy conservation and improving electrical device performance, may be accomplished. The possibilities of energy harvesting, through exploiting anisotropic thermopower in the metamaterials is also examined. The review concludes with a brief survey of the outstanding issues and insights needed for further progress. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http: // dx.doi.org / 10.1063 / 1.4916220]


I. INTRODUCTION
A significant contributor to energy wastage is the inevitable generation as well as the uneven dissipation of heat. Practical methods to adeptly channel heat flux would have many applications enabling improved energy utilization and thermal energy management. The first step would be to understand the notion of heat flux. From a thermodynamics viewpoint, the transfer of thermal energy due to temperature differences has been well codified. Based on traditional engineering antecedents, the subdivisions into particular mechanisms of heat transfer, as related to conductive, convective, and radiative pathways have also been made. While such partitioning has never been completely clear, scientific investigations have long had the objective of unifying such divisions. For instance, thermal conduction and convection has almost always been treated through a continuum basis, such as through the deployment of the diffusion or the energy conservation equations, while radiative phenomena have used quantum mechanical formulations incorporating both particle-and wave-like characteristics. In this context, much progress has been made at the nanoscale (in both length and time), where ballistic transport of heat could be used to describe phonon-mediated conduction. However, there seem to be very little practical consequences for the macroscale in terms of practical spatial and temporal scales, at which heat management is typically an issue.
The aim of this paper is to propose the use of relatively large-scale thermal metamaterials for thermal energy harness. A scientific objective is to describe thermal energy transfer on a unified basis, which may be provided through a metamaterials type description. Considering the correspondence of phonons to photons, both from a fundamental particle point of view and their role as energy carriers in and out of materials respectively, and given the wide scale usage of optical elements such as lenses, waveguides, etc., it is the thesis of this paper that the manipulation of thermal energy would also be possible through such metamaterials. a Electronic mail: pbandaru@ucsd.edu The metamaterials terminology itself seems to have originated from the electromagnetism (EM) community, for the manipulation of EM fields through arrangement of specific materials. For example, the judicious placement of materials with a negative dielectric constant and negative magnetic permeability yielded unique phenomena related to negative refraction, with the possibility of bending light the "wrong way". 1-4 Consequently, any diverging characteristics could be minimized yielding perfect lenses. The possibility of thus controlling electromagnetic fields 5 has yielded in the past decade, even more fascinating possibilities ranging from higher efficiency antennae to invisibility cloaks. 6 There is then an obvious imperative to mimic such aspects to thermal applications as well, e.g., in thermal cloaks, [7][8][9] which would forbid the entry of thermal flux, or nanoscale heat sources.
In this article, we take a more fundamental perspective in investigating the rationale for thermal metamaterials from the bottom-up. The notion of arranging materials with nominal isotropic attributes to yield unique characteristics would be a focus. The necessity of considering intrinsic materials properties was also felt in the gradual progress of technology to investigating nanoscale materials. For example, while the wavelength of light was considered to be a limiting characteristic for the optical lithography patterning of electrical transistors per the Rayleigh resolution limit, it was gradually realized that such limits might not necessarily hold. 10 For example, in the electrostatic limit, when electric (or magnetic field) interactions at the dipole level are considered 11 the limits could be considerably extended. The agency of surface plasmon excitations 12 could now yield smaller length scales. Related quantum mechanical considerations on tunneling could also allow for classically forbidden size limits. Such new size scales, as related to the new excitations being investigated, in concert with intrinsic materials properties, have then formed the basis for the arrangement into metamaterials.
While quantum mechanical or wave-like characteristics of thermal conduction have not yet been considered, the field of thermal metamaterials has evolved over the past few years more on the basis of the judicious arrangement of materials, which serve to direct heat depending on the thermal conductivity characteristics of the constituents. The thermal metamaterials area has generally exploded in the past decade, almost wherever wave motion has been considered, with a range extending from low frequency elastic 13,14 and acoustic waves 15,16 to high frequency EM waves. 2, 17 While the extension to diffusive transport typical to heat motion was at first tenuous, nevertheless, engineered materials arrangements were shown in the year 2012 to accomplish functions such as heat flux focusing and re-direction. 18 The motivation in the related seminal paper by Narayana and Sato 18 yet seemed to rely on the concepts of coordinate transformation as first enunciated for EM metamaterials. 6,19,20 While relevant and interesting analogies do exist, e.g., in two dimensions, there could be an equivalence 21 between acoustics (with very long wavelength phonons/heat transport) and electromagnetics in isotropic media, the relationship is not very clear. 16 A major issue is that the wave-like character and the symmetry implicit in the Maxwell equations is not apparent in the Fourier law of heat conduction, 16,17 which is governed by diffusive phenomena. Consequently, much of thermal transport seems to be governed by non-resonant elements.
It is then the aim of this paper (which reviews many recent results and incorporates an overview of recent publications [22][23][24][25][26][27][28][29] from our group) to indicate that significant insight and control into the directionality of heat transport may still be exerted through considering the geometric characteristics of thermal conductivity. We show utilizing fundamental tenets, such as (a) that the net thermal conductivity of a composite can be altered through the arrangement of isotropic materials, and (b) the adaptation of the Fermat principle of least time, whereby heat flux choses the path of least resistance, may yield novel arrangements through which control over conductive heat flux may be achieved. Clearly, the motivation for such studies has been whether an alternative approach for understanding heat transport control, without resorting to EM methodologies, would be feasible and yield equivalent results. Another imperative has been the proposition of an effective thermal medium (ETM), equivalent to the analogous concept in electromagnetism, for considering the arrangement of the materials within the composite. While much of the discussion is on macroscale materials, principles of scale invariance may be applied to scale up or down the structures. While interfaces between the materials in the composite may indicate the necessity of nanoscale features, the ETM approach makes such issues unnecessary. We also hope to indicate whether truly nanoscale heat sources 30 could be feasible, e.g., as to what constitutes the minimum spot size to (/from) which heat may be generated.
We start by examining in detail the fundamentals and modalities of conductive heat flux transport in materials, which would then be extended to metamaterial-like structures. While the theoretical underpinnings will be first probed, this would lead the way to experimental verification of the analytical formulations.

II. UNDERSTANDING AND MANIPULATING THE DIRECTION OF THE CONDUCTIVE HEAT FLUX
We first discuss the fundamental aspects related to the tensorial form of the Fourier law of heat conduction and its manifestation in an effective thermal medium (which lays the foundation for a thermal metamaterial, generally). Subsequently, analytical and experimental investigations pertaining to the bending of heat flux in a given direction will be discussed.

A. The thermal conductivity tensor
It has been very well known, through the Fourier law of heat conduction, 31 that the thermal conductivity (κ ij ) is a second-order tensor, and is defined through the ratio of the flux in the i th direction (q i ) to the temperature gradient in the j th direction (∇T j ), using: In adopting the above formulation, we do not explicitly consider spatial and temporal effects, that may, for example, be pertinent to the nanoscale, 32 which may invalidate the Fourier law, 33,34 e.g., issues related to whether a temperature gradient may be defined at length scales smaller than a heat carrier mean free path 35,36 or the consideration of time scales 37 for heat flux transport, etc. Consequently, diffusive transport is the main consideration. It will be seen later that while such treatment is not complete, e.g., in terms of the coupling of the conductive heat to radiative heat transfer, it would provide an initial advance beyond present usage. Additionally, from a practical point of view, the thermal conductivity has been typically regarded as an engineering parameter, and as a scalar and isotropic materials property. Anisotropy in the thermal conductivity has not been generally considered even in confined nanostructures such as nanowires. 38 In the case of two-dimensional thin films, a few attempts have been made to correlate possible anisotropies to film and measurement geometry, 5 where for example, a conductivity variation dependent on the direction heat flow, whether cross-plane (κ cp ) or in-plane (κ ip ), was noted. For example, the possibilities and experimental demonstration of phonon focusing, [39][40][41] consequent to the elastic anisotropy in crystalline solids, implicitly considers the non-scalar nature of the thermal conductivity with as much as a factor of two variation. 40 While the assumptions of classical mechanics are clearly inadequate to explain electron motion, due to the spin degree of freedom, classical modeling has been typically considered adequate for thermal transport.
However, we will now show that considering the full tensorial nature of the thermal conductivity leads to a rich variety of behavior that may be exploited for manipulating the directionality of the heat flux. Such understanding could further lead to ideas into how an optimized arrangement of materials (i.e., as metamaterials) may facilitate the directing of heat flux, leading to a wide variety of possible uses including heat concentration as well as heat dispersion leading to both efficient thermal energy utilization as well as dissipation.
We first illustrate the anisotropic nature of the thermal conductivity with respect to a orthogonal, e.g., rectangular (x-y-z) coordinate system, i.e., In engineering practice, with the reporting on materials as if they possess isotropic thermal conductivity, 42,43 i.e., κ ij is equated to a single scalar value (≡ κ, say). In this case, the heat flux density vector (q i ) follows the respective temperature gradient, i.e., the heat flux in the horizontal/x-direction is only determined by the temperature gradient in that direction. The presence of off-diagonal terms, i.e., κ ij with i j, would induce cross-coupling and concomitant bending of the heat flux, e.g., with a substantial κ xy , the heat flux in the x-direction would be determined by the temperature gradient in both the x-and an orthogonal (i.e., y-/z-) direction. Such considerations immediately suggest controlled manipulation of the thermal flux and the placement of the materials, which facilitate such manipulation, would be highly desirable.

B. A two-dimensional composite thermal metamaterial
From a photonics viewpoint, a simple alternating arrangement of two materials (of distinct dielectric constants: ε) has been invoked as a prototypical example of a composite possessing a band structure along with a photonic bandgap 44 -the magnitude of which is proportional to the relative difference between the respective values of the ε. It would be interesting to make an analogy to the effects on the thermal conduction that could be caused due to two materials (of distinct thermal conductivity: κ). However, we first consider the diffusion equation when a steady state temperature distribution is prevalent. Consider a layering of two materials, e.g., when two sheets with nominally isotropic thermal conductivities κ 1 and κ 2 are alternatively stacked -as depicted in Figure 1(a). The structure represented in this figure is a prototypical example of layered media 45 that, for example, has laid the foundation for photonic crystals. 44 The κ 1 /κ 2 ratio is assumed to be greater than unity, and layers 1 and 2 are interchangeable. Consequently, the material with the higher thermal conductivity should be taken to have a value of κ 1 . The case of κ 1 /κ 2 = 1 corresponds to a homogeneous material. Considering one-dimensional parallel and perpendicular heat transport, between the surface on the right (maintained at a temperature: T h ) and the left surface (at a lower temperature: T c ), and assuming flux continuity and neglecting interfacial effects, it can be easily derived 22 that: FIG. 1. (a) Illustration of an anisotropic composite, constituting a metamaterial, of total length, l, which may be fabricated by alternatively stacking two thin sheets of thermal conductivities κ 1 and κ 2 . The left-hand side is maintained at a uniform temperature: T hot , while the right-hand side is at a lower uniform temperature, T cold (< T hot ). ) is always less than or equal to the transverse thermal conductivity, (κ y = κ 1 +κ 2 2 ), as the harmonic mean is less than or equal to the arithmetic mean.

C. The concept of an effective thermal medium (ETM)
A crucial criterion in the definition of a metamaterial is whether the composite fabricated, say through the alternating layering of materials as in Figure 1, may truly be considered a new material (e.g., possessing an anisotropic thermal conductivity, while being effectively homogeneous in its overall characteristics) for the relevant physical phenomena. Such an issue is relevant as interfacial thermal boundary resistance 46 has been widely acknowledged as one of the major impediments 47 to efficient heat transfer. Consequently, if interfaces play a major role then the composite may not ideally be considered to be a metamaterial, as interfacial variation could affect the characteristics. For the case of optical metamaterials, an arbitrary but widely accepted definition 48 is that the wavelength: structure (e.g., the layers) ratio should be 10:1, such that the specific composite features are not discernible and the metamaterial is truly a unique material. As thermal conduction is predominantly governed by diffusion, such a criteria seems difficult to define, given the lack of a characteristic wavelength/ length scale, e.g., as in the difficulty of defining a precise or relevant value of the mass (m) in the de-Broglie thermal wavelength: 32 Alternately, if it were assumed that the interfacial thickness was sufficiently small, this would be equivalent to the tenet that a linear temperature gradient can be defined. Considering a net temperature gradient T dev (= 2κ 2 κ 1 +κ 2 ∇T) along the total length of the unit comprising layers 1 and 2 (of Figure 1), it can be derived that the gradient across the individual layers would be ∇T 1 (= 2κ 2 κ 1 +κ 2 ∇T) and ∇T 2 (= 2κ 1 κ 1 +κ 2 ∇T), respectively. We have assumed that the n (an even integer) layers are constituted from a unit cell comprising a layer pair. Consequently, the temperature difference (T dev ), between that obtained from assuming a linear temperature gradient across the length, (equivalent to defining an effective thermal conductivity for the composite) and that considering temperature variation across the individual layers, would be: 22 It is to be noted that T dev decreases as n increases and is proportional to the κ 1 /κ 2 ratio, and consequently approximating the composite medium to possess an equivalent thermal conductivity is valid for a small layer thickness (=l/n) as well as for a higher thermal conductivity contrast, i.e., the approximation is exact for infinitesimally small thickness and when the composite is effectively constituted from a single material, say with κ 2 → 0. Eqn. (4) above is then a useful relation for understanding the basis for the effective medium approximation in the case of thermal metamaterials.

D. The re-direction of the heat flux
We now consider the application of a temperature gradient to a thermal metamaterial modeled as an ETM. When the orientation of the layers in the composite are rotated, in the plane, say, by an angle θ(−π/2 < θ < π/2, the layers are now oriented along new axes: x ′ (= x cos θ + y sin θ) and y ′ (= −x sin θ + y cos θ): see Figure 2.
It is to be noted that the angle θ is the amount of rotation necessary to align the coordinate axes fixed to the sample (x, y) to the coordinate axes along which the temperature gradient is applied (x ′ , y ′ ) and is considered to be positive in the counter-clockwise direction. The originally applied (horizontal) temperature gradient now acquires both horizontal and vertical components, with respect to the rotated layers. The resultant heat flux deviation may be interpreted in terms of rotating a thermal conductivity ellipsoid -the orientation of which in the unrotated guise is manifested through the elements of Eqn. (3). The degree of anisotropy would change from one composite sample rotation to another and the consequent change may be viewed in terms of tuning the κ i j . x, y and x ′ , y ′ correspond to original and rotated coordinate systems respectively. The κ 1 and κ 2 of the constituent layers were selected with representative values of 0.1 W/mK and 1 W/mK, respectively, and the surfaces on the right and the left hand side are maintained at 350 K and 300 K. The end effects and heat loss due to convection and radiation have been ignored. Square geometries were assumed. The temperature isotherms on the far right in Considering that the heat conduction equation, from Eqn. (1), can be written as ∇ · (κ ij ∇T j ) = 0, in the absence of heat sources/sinks in the ETM, the modified thermal conductivity (κ ij m ) in a changed coordinate system would be: 22 where J is the Jacobian matrix of the coordinate transformation between the new and the old coordinate systems, J T is the transpose of J and det (J) denotes the determinant. Consequently: A modified Fourier law could then be written as: The extent of heat flux bending, as inferred through the deviation (φ) of the q i would then be related to the ratio of the terms in off-diagonal to the diagonal terms in the κ ij m , and would be: In the expression above, c = κ y /κ x = (1+κ 1 /κ 2 ) 2 4×κ 1 /κ 2 . We note that φ (with a period of π) could be either positive, negative, or zero, which allows us to direct the heat flux density upwards or downwards through changing/tuning θ and/or varying κ 1 and κ 2 . Such upward/downward bending of the heat flux (corresponding to positive/negative φ) is shown in the computational simulations (using COMSOL ® Multiphysics) of Figure 2. This modulation of the heat flux bending has also been experimentally demonstrated (see Section II E), and found to be in accordance with the derived analytical relations. Figure 3 indicates that φ initially shows a linear variation with κ 1 /κ 2 , as a function of θ, and finally saturates at an angle equal to (90 • − θ), as apparent for larger values of θ. This should also be clear from Eqn. (8), which is increasing function of κ 1 /κ 2 and approaches a maximum value at κ 1 /κ 2 → ∞, with φ reaching a value = 90 − θ. At this limit, e.g., with κ 2 → 0, the heat flux would be forced to traverse layer 1, and be bent by an angle of 90 0 . A maximum contrast in the κ 1 /κ 2 ratio enabling such marked channeling may be considered analogous to the required dielectric constant contrast for observing significant band gaps between the propagating and non-propagating regimes, in photonic crystals. 49 Alternately, considering a minimum (say, with κ 1 /κ 2 → 1) would yield a φ = 0. At lower values of θ, the resolution of the horizontal temperature gradient along the rotated layers would be less well defined and higher order/cross-coupling effects may be important.
We infer that a larger contrast in the thermal conductivity values between the two layers enables a larger heat flux rotation at a smaller degree of composite rotation angle. As the range of thermal conductivity ratios, for well-understood, practical materials in vogue 50 could span five orders of magnitude, i.e., from ∼ 0.03 W/mK (for silica aerogel) to ∼ 2200 W/mK (for diamond) such characteristics may be easily demonstrated. The layers may be assembled in a composite-like form, for example, through an electrical discharge machining (EDM) process, 24 where grooves could be cut on a base plate of, e.g., copper (κ = 390 W/mK), and the spacing/channels in between filled with relatively lower (/higher) thermal conductivity material, e.g., PDMS (poly-dimethylsiloxane, κ = 0.17 W/mK) was used, to yield a high thermal conductivity contrast or a κ 1 /κ 2 ratio. A large length/width to thickness ratio (of >10, say) would ensure that heat flow can effectively be considered two-dimensional and simplify the κ i j (given by Eqn. (3)) effectively to a 2 X 2 matrix. Additionally, it has been assumed that the in-plane radiative heat transfer between the metal layers was negligible, based on the radiation flux calculated using the classical 51 Stefan-Boltzmann formulations or as modified through phonon radiative transfer mechanisms. 34 Such preliminary considerations clearly indicate the possibility of directing the heat flux in a different direction than that in which the thermal gradient is applied through carefully arranging the constituents in the composite. Such a notion is at the heart of constructing a metamaterial with novel functionality.

E. Experimental demonstrations
Considering that a major motivation for the design of thermal metamaterials was to implement optical equivalents such as lenses, we now indicate how the concomitant bending of heat flux could be accomplished through the use of such structures. The broad idea is to demonstrate the possible equivalent of Snell's law -now for the bending of the heat flux vector.
For the purpose of demonstrating proof of principle, thin layers (t = 2 mm) of commonly used engineering materials such as copper (κ 1 = 390 W/mK) and stainless steel (κ 2 = 42 W/mK) were alternately stacked (n = 30) at three different orientation angles, θ = 0 • , 45 • , and −45 • and placed in a temperature gradient of ∼ 80 K (T hot = 100 • C, and T cold = 22 • C) with stable heat sources enabled through hot and cold water baths, with large heat capacity. The composites (see inset to Figure 1(b), were fabricated in a square geometry with a length and width of 5 cm, and have a thickness of ∼ 0.3 cm. To minimize the perturbation of the heat flux at the interface with the background, the composite was embedded in brass (κ b ∼ 120 W/mK) to correspond to the geometric mean of the thermal conductivities 52 of the metamaterials' constituents. The large ratio of the length/width to that of the thickness, of ∼ 16, also enables us to treat the composites as an effective two-dimensional ETM, and simplifies the κ i j (given by Eqn. (3)) and experimental demonstration of the proof of concept. Thermally conducting silver epoxy (κ Ag ∼ 5 W/mK) was used to join the layers and in the ETM approximation (with a ∆T < 1 K, from Eqn. (4)) would not be involved in altering much the observed steady state and transient heat flux. The temperature profiles -both in the steady state and transient cases-were taken using an infra-red (IR) camera (FLIR 320), with a spatial resolution of ∼ 300 µm and a temperature resolution of 0.25 K. For obtaining an enhanced signal, the top surfaces of the composite were painted with a thin layer of black paint (of emissivity 0.95). Figure 4 reports the experimental results, in steady state, of the thermal flux through the metamaterials. For the case of θ = 0 • : Figure 4(a), the q i are seen to be perpendicular to the isothermal contours, as expected from Eqn. (1). However, with θ = 45 • , the effective thermal conductivity tensor can be calculated to be anisotropic, from Eqn. (6), as per: The q i are no longer perpendicular to the isothermal contours and were bent by an angle (φ ∼ 26 • ) with respect to the horizontal temperature gradient, due to the presence of a transverse component. The experimentally obtained φ was in excellent accord with the value computed through using Eqn. (8), where φ was computed to be ∼ 24 • (assuming the presence of a finite interface of thickness ∼100 µm) or ∼ 26 • (assuming negligible interfacial effects between the copper and the stainless steel). It is also to be noted that while for a material with isotropic thermal conductivity the temperature gradient tends to a constant value in the steady state, this would not be the case for a material with an anisotropic κ ij as observed in the Figures 4(b) and 4(c). Alternately, with θ = −45 • : Figure  4(c), negative off-diagonal components were evident in the thermal conductivity tensor, from Eqn. (6), per: The q i are now bent by an angle (φ ∼ −26 • ) with respect to the horizontal temperature gradient. The changing of the heat flux propagation directions in Figures 4(b) and 4(c) are reminiscent of the bending of light at the interface of two media. In the optical case, the propagation follows Snell's law, with propagation towards/away from the normal in the case of transmission from an optically less (/more) dense medium to a more (/less) dense medium. Considering an analogy of the thermal resistivity (the inverse of the thermal conductivity) to the optical density, a similar behavior was then observed, following the tenets of the thermal resistance minimization principle.
Indeed, an equivalent of Snell's law 53 had been previously derived 54 for heat propagation from one incident isotropic medium to another medium. This was done through analogy to the tangent law for electric field lines passing between two materials of different dielectric constants, 55 a tangent law for refraction for heat conduction was previously derived, 54 whereby for a given angle of incidence φ i , the angle of refraction φ r .
However, it was previously shown 22 that an alternative formalism considering the specific anisotropic character of the metamaterial captures the propagation behavior more accurately. Such detailed formulations also illuminate the possibility of negative refraction, 23 as will be considered in the next section, where the transmitted heat flux could bend anomalously with respect to the interface normal. An excellent correspondence between the simulations and the experiments (as in Figure 5) was observed, e.g., the experimental heat flux bending angle φ for the anisotropic composites oriented with θ = 45 • and −45 • were ∼ 26 • and ∼ −26 • , respectively while the simulations yielded values of φ ∼ 23 • and −23 • , closer to the computed value of ∼24 • (as discussed previously).
We have previously considered the notion of an ETM, where the individual constituents/ interfaces between the constituents would not matter in describing the passage of heat through the metamaterial. However, while experimental observations in Figures 4 & 5 seem to indicate such an idealization, simulations indicate wiggles in the heat flux -as indicated for example in the finite element method (FEM) based simulations of Figure 6.
Such modulations were ascribed to the successive and periodic refraction 23 of the heat flux at the interface of the two constituents in the metamaterial. The angle of refraction can be quantitatively estimated through considering the diagonal as well as the off-diagonal components of the thermal conductivity tensor of the composite as well as the oblique angle of incidence, 23  previously discussed. It would then be pertinent to investigate the circumstances leading to such heat flux refraction and the consequent applications.

F. Applications
The notion that thermal energy and the concomitant conductive heat flux may be re-oriented immediately suggests applications of the metamaterials discussed this far, such as thermal concentrators -where heat is a focused, or thermal cloak -where heat is avoided in a certain region: Figure 7.
The arrangement of designer bending blocks in a modular fashion: as depicted in Figure 7, would make it feasible to enhance (as for a thermal concentrator) or diminish (as for a thermal cloak or a diffuser) heat transfer to a particular region. There is also the possibility of thermal diffusers, which dissipates heat from a certain area, of use for example, in the cooling of electronic devices that reside on the region that would be affected. While we have thus far considered the application of one-dimensional temperature gradients to bulk objects, the principles may easily be extended to higher dimensions with similar functionality. Indeed, the thermal cloaking of bulk objects has been shown experimentally. 7,8 The ability to create such metamaterials would bring forth new opportunities for the creation of thermal elements and devices and forge new frontiers in the channeling of thermal energy, energy conservation as well as the harness of heat.

III. NEW DIRECTIONS IN CONDUCTIVE HEAT TRANSPORT IN METAMATERIALS
In addition to directing the heat flux, following the anisotropy in the thermal conductivity tensor, as discussed in the previous section, further probing leads to the possibilities of bending heat flux (either as expected or anomalously). Additionally, the anisotropy in the thermal conductivity may be related to the thermal diffusivity and intriguing issues related to the manipulation of heat including its stoppage. Such issues along with the possibility of constructing thermal lenses in analogy to optical lenses will be discussed in this section.

A. The refraction of heat flux: Positive and negative refraction
The phenomenon of negative refraction, 2 when first announced in the optical physics/EM context, sparked much excitement and intrigue in terms of the possibilities for a perfect lens. 2 The basic idea was that the amplitude and the phase of the evanescent waves could be resurrected through the agency of a negative refractive index, whereby a decaying wave-vector could be effectively reversed. However, the amplification of the evanescent waves (yielding a transmission as opposed to a decay) in the near field was more due to resonant coupling to surface plasmon-like excitations. 56,57 While it may be speculated that surface phonon polaritons (SPhPs) 58 -formed through the coupling of long-wavelength EM waves to optical phonons may serve the role of resonant excitations for the case of thermal energy transport, such hypotheses need to be much more carefully examined.
Nevertheless, in correspondence to the normally encountered refraction of light, we may define positive heat flux refraction to occur when the horizontal (x) and the vertical ( y) components of the incident and refracted heat flux vector points in the same direction and have the same sign. 59 In this context, we have shown in the previous section that metamaterial characteristics may be controlled to yield either a positive or negative heat flux bending angle, i.e., φ >/< 0. For example, if two thermal metamaterials -both with distinct and positive φ, were placed adjacent to each other, the AIP Advances 5, 053403 (2015) passage of heat flux would be as shown in Figure 8(i). Such a scenario would be termed as positive refraction. Consequently, the placement of metamaterials with contrary signs of the φ, may be thought to yield heat flux bending that is considered to correspond to negative/anomalous refraction, as indicated in Figures 8(ii) and 8(iii). However, the matching of the heat flux at the interface of the two materials must be carefully considered.
Assuming steady state conditions to prevail, we may derive the specific materials characteristics and arrangements that accomplish such heat flux passage through the composite. with no heat sources/sinks), along with (b) the continuity of temperature gradient along the interface (assuming no interfacial resistance). Consequently, an angle of refraction (φ r ) in the second media (subsequent to passage through the first) may be derived as: The ρ components can be derived using ρ ij = ( ρ xx ρ xy ρ yx ρ yy ij . It can then be deduced -from Eqn. (11) that the φ r may be positive/negative, according to whether: , respectively. The interesting aspect is that ordinary materials may be arranged to yield such non-intuitive attributes. We illustrate positive (φ r > 0) and negative refractions (φ r < 0) through considering two isotropic and commonly available materials, e.g., with κ 1 = 20 W/mK (e.g., carbon steel) and κ 2 = 0.1 W/mK (e.g., polystyrene) of equal thickness stacked alternately as indicated in Figures  1 or 8. Generally, the extent of bending or refraction seemed to be proportional to the thermal conductivity contrast between the constituents. The effective thermal conductivity of the composite in Figure 7 It was previously shown 22 that when the composite in Figure 1 was rotated in the x-y plane by θ (with −π/2 < θ < π/2) -see Figure 2, that the thermal conductivity would be modified as per Eqn. (8). We consider representative degrees of rotation of such a bilayered composite, e.g., by θ = π/4, π/3, and −π/6, yielding the following thermal conductivity matrices in the three cases respectively, with obvious anisotropy: We see straight away, e.g., in (C) above, that some of the elements in the matrix are negative. It should be cautioned that the signs (positive or negative) of the off-diagonal components are due to the operation of rotation and choice of coordinate system and may not be ascribed physical significance. However, metamaterials fabricated through stacking various combinations of the above cases yield diverse phenomena corresponding to the positive or negative bending of the thermal flux. For instance, the passage of heat through the metamaterial structures corresponding to stacking of combinations (A)/(B) and (A)/(C), and (B)/(C), are illustrated in Figure 8. The values of φ r were calculated from Eqn. (11). We have chosen these particular combinations as we found that φ r > 0 in the first case (representative of positive/normal refraction) while φ r < 0 (indicative of negative/anomalous refraction) in the latter two cases. Such conclusions were validated through simulations conducted using COMSOL ® Multiphysics software. For the simulations, the surfaces on the top and the bottom of a semi-infinite composite (of total height = 10 cm) were maintained at 300 K and 350 K respectively. Additionally, the material interfaces were assumed to be smooth without interfacial resistance and end effects, and heat loss due to convection and radiation has been ignored. We justify the assumption of negligible thermal interfacial resistance based on the rationale that (a) the interface is sufficiently thin, and (b) that no heat is stored or dissipated in the interfacial layer, in accordance with the boundary conditions specified earlier. Additionally, for such a two-dimensional semi-infinite composite, which is thin, the temperature gradient in the transverse direction tends to zero. It was also observed that the transmission phenomenon at the interfaces is reciprocal, i.e., reversing the direction of temperature gradients reverses the direction of the heat flux vectors preserving the incident and transmission angles.
It should be pointed out that the described situation for heat flux passage are yet distinct from the refraction related phenomena associated with a traditional optics viewpoint. Here, where ray propagation follows Snell's law, negative refraction (with θ r < 0) occurs for all incident angles (θ i ) -from θ r = Sin −1 ( n i n r Sinθ i ) if n r < 0 or n i /n r < 0, where n i and n r are the respective optical refractive indices of the incident and refracted media. As such effects are not always observed in our investigations, e.g., negative refraction (φ r < 0) only occurs for negative arguments of the Tan −1 function in Eqn. (11) when φ i < 0, the reverse bending may be termed as anomalous refraction. Our formulations can also be related to theoretical 60 and experimental 61 efforts in literature invoking "transformation thermodynamics," 60 which in essence relate to geometric transformations 5 for modulating the isotropic thermal conductivity in polar coordinate systems for obtaining circular thermal cloaks, where such reverse bending has been noted. The employed rectangular coordinate system is pertinent to considering the refraction of a heat flux at a given point of the ambient/metamaterial composite interface. A rotation of the thermal flux 62 would occur through a continuous change of the φ r . Incidentally, when φ r (= φ i − π) < 0, a backwards refracted/propagating heat flux could be obtained.

B. Thermal conduction in transient -The stopping of heat?
To further understand the ETM behavior of such anisotropic metamaterials, the transient propagation of the thermal flux was recorded and is shown in Figure 9. The development of the thermal flux profiles, as recorded at t = 90, 150, and 210 s, and showing the bending of the heat flux lines indicates the invariance to coordinate transformations 62 of the time dependent heat diffusion equation: ρ(r)C(r) ∂T (r ) ∂t = ∇ · [κ(r)∇T(r)], where ρ is the material density and C is the specific heat, and r is the spatial coordinate. Additionally, the thermal boundary resistance associated with the passage of heat flux from the background material (brass) to the metamaterial (treated as an ETM), due to dissimilar thermal impedances, 46 and the associated relative temperature drop at the given flux was seen. The thermal diffusivity: D th when formally defined through the ratio of the thermal conductivity to the product of the ρ and C could also be seen to be anisotropic given the constancy of the ρC product on geometric transformation. 60 The ρ and C were estimated from the average of the respective values of the constituent (stainless steel, with ρ SS = 7.8 g/cm 3 and C SS = 500 J/kg · K, and isotropic thermal conductivity, κ SS = 42 W/mK, and copper, with ρ Cu = 7.8 g/cm 3 and C Cu = 385 J/kgK, and isotropic thermal conductivity, κ Cu = 391 W/m K,) materials used for the layers in the metamaterials. If it is assumed that the diffusivity lies within an upper Wiener bound 63 assuming a parallel configuration for the layers (with temperature gradient along the layers) -with ρC = ( ρ Cu C Cu +ρ SS C SS 2 ) ∼ 3.7 X 10 6 J/m 3 K) and a lower Wiener bound 63 assuming a serial configuration for the layers (with temperature gradient perpendicular to the layers, with ρC = 1 2 ( ρ Cu C Cu κ SS +ρ SS C SS κ Cu κ SS +κ Cu ) ∼ 1.9 X 10 6 J/m 3 K). Simple averaging then results in a (ρC) av = 2.8 X 10 6 J/m 3 K) and a computed D ij equivalent to: 3.6 1.6 0 1.6 3.6 0 0 0 5.2 Such anisotropic diffusivity, hinging on the relative difference between the elements of the κ ij implies unequal times for the propagation of heat in such thermal metamaterials, e.g., with reference to the matrix Eqn. (13), the times for the heat to traverse a given length, l is of the order of l 2 /D th and would be in the ratio of 1:0.75:0.33 for propagation in the brass, the composite in the horizontal direction, and the composite in the vertical direction, respectively. While the D th for the brass was computed to be 3.6 X 10 −5 m 2 /s from κ,C, and ρ values taken from literature, 50 for the composite the κ ij was from Eqn.
(3) while C, and ρ were estimated from the average of the respective constituent (steel and copper) materials used for the layers in the metamaterials. The staggering of the temperature profiles, in Figure 8, over a time period of a few minutes at the background-metamaterial interface was also suggestive of a very low diffusivity as well as the possibility of nonlinear effects. [64][65][66]

C. The possibility of a thermal lens
The physical phenomena presented thus far, indicating the routing and the refraction of thermal flux, takes us into the realms of considering thermal lenses. We rwere motivated by the basic definition of a lens, as one "that reconfigures a transmitted energy distribution," 53 and its working principles, i.e., the geometric/layered arrangement of a material property such as the refractive index (the analogue of which in the case of heat, could be the thermal resistivity or its inverse, the thermal conductivity).
For the implementation of such a thermal lens, we adopt a variational methodology equivalent to Fermat's principle, 6 with the hypothesis that the thermal energy propagation would be such that the spatial integral of the thermal resistance (R th ) over a unit length of propagation (dl), i.e.,  R th dl, would be a minimum/extremum 54 (where R th = 1/κ A, with A as the heat flux pertinent area). The minimization would be based on the arrangement of media/materials through which the heat flux (Q) propagation would occur. It was proposed that the constitutive material character (through the κ) as well as the extensive parameters such as the thermal traversal length (L) as well as the A can be suitably adjusted to achieve the desired directionality of Q as well as the extent and magnitude of the focusing. Consequently, thermal energy re-orientation could be accomplished leading to functionalities such as energy focusing or dispersion akin to those paradigms familiar from basic optical lens phenomenology.
We will initially illustrate a prototype of a thermal concentrator, in almost exact analogy with a light-concentrating lens. Generally, when exposed to a light beam, an optical lens unit needs to perform the function of energy focusing, with the spatial extent of the focused spot related to the effective wavelength and the numerical aperture. 67 However, as incoherent heat transfer phenomena is generally 68 wavelength insensitive/diffusive the focal spot would be less well defined. The gradual increase in the optical path length,  ndl (where n is the optical refractive index) from the center to the edge of a typical optical converging lens was mimicked in the experiments on thermal lenses, through a similar variation of the κ. Ideally, such variation would need to be implemented with a very large number of layers and concomitant decrease in the layer thickness (t), and would be implicitly related to the assumption of an ETM (see Section II C) for the composite thermal lens. The thermal lens characteristics could then arise, e.g., (1) due to a change in the κ, keeping the Q traversal length constant, or (2) due to the variation of the geometry in the form of a change of the AIP Advances 5, 053403 (2015) FIG. 10. The variation of the output heat flux (Q) with θ (the angular position of a particular element of uniform κ along the periphery of the lens, reckoned with respect to the horizontal). For a given temperature gradient and geometry, the extent to which heat could be focused was found to be the largest for (ii) a CTL (concave thermal lens) formed using layers with varying κ from the center to the edge, in comparison to a homogeneous κ constituted (i) CTL or (iii) flat/planar geometry. L, at a constant κ. While the former case corresponds to a planar surface, the latter would imply a curved/lens-like geometry.
We indicate, in Figure 10, the modeled variation of the heat flux for the lens-like geometry, where a concave thermal lens (CTL), was composed of either (i) bulk material with isotropic/ homogeneous κ or (ii) layers with varying κ from the center to the edge, in comparison to a (iii) flat/planar geometry (with homogeneous κ). It was clear that, for a given temperature gradient, the Q enhancement was maximal for (ii). The area under the curves (proportional to the power) in the figure was identical in all the three cases, ensuring the conservation of Q. The symmetry of the CTL profiles may be related to the dependence of the Q to κ(θ)/L cos θ, where θ refers to the angular position of a particular element of uniform κ along the periphery of the lens, reckoned with respect to the horizontal, at the center of curvature of the CTL. While κ(θ) was variable for (ii), it was constant for (i).
For experimental demonstration of the proof of concept that spatial arrangement and geometry would indeed be relevant for heat flux focusing, polymeric nanocomposites (infiltrated with non-agglomerative 69 addition of carbon nanotubes to reactive ethylene terpolymer matrix 70 ) were chosen for the thermal lens materials. Such composites could relatively easily accommodate a gradual variation in the κ, which was necessary for creating an ETM. 71 The use of polymeric material implies low economic costs coupled with simple fabrication, and could at the very outset be related to the harness of low-grade waste heat.
We compare results from a planar thermal lens (PTL): Figure 11(a): Simulations, and Figure  11(b): Experiment, with κ decreasing linearly from the center to the edge, to that of a CTL: Figure  11 Figure 11(c), while a smoothly varying κ(θ) was desirable the experimental arrangement yielded more reliable results through the use of wedge-shaped samples. Conformal mapping techniques could be used for mapping the wedges onto circular geometries 72 implying similar geometric characteristics and energy focusing attributes as those of spherical lenses. 20 For the experimental investigations, the input Q to the thermal lens (at the bottom of the sample in Figure 11) was maintained at a uniform and constant value of ∼ 500 W/m 2 through a thick block of acrylic. Constant temperatures were maintained at the heat source and sink, employing a hot plate heater and a custom built Peltier cooler. A large block length ensured that entrance effects 73 could be ignored. Correspondingly, a copper block with a conforming shape was placed adjacent to the lens surfaces in both the PTL and the CTL arrangements and served as a heat collector. Both the Q source (acrylic block) and the Q sink (copper block) function in a manner similar to electrodes for electrical or electrochemical devices, serving as large heat capacity reservoirs. Isothermal conditions were ensured by holding the bottom (the hot side) of the acrylic at T hot = 380 K and the top (the cold side) of the copper block at T cold = 278 K. The heat conduction temperature profiles were non-invasively imaged, using an infra-red camera (FLIR A15 SC, with a spatial resolution of ∼ 0.3 mm and a temperature resolution of ∼ 0.25 K), with appropriate calibration. To achieve good signal fidelity through enhanced thermal emissivity from the surface of the composite, the top surfaces were painted with a thin layer of black paint, as previously indicated in Section II E.
The enhanced Q concentration through the use of thermal was evident through the FEA simulations (left panels) as well as in the experiments (right panels) of Figure 10. For the simulations, corresponding to the experimental situation, constant temperature boundary conditions were assumed at the entry and exit while the sides were insulated. The Q contours in Figures 11(b) and 11(d), were oriented to be perpendicular to the isothermal lines -obtained through IR imaging, and clearly indicate a focusing towards the center in both the PTL and the CTL cases. Such heat flux concentration is graphically indicated through the plots of Figure 12.
The extent of Q concentration was significantly larger through the use of the CTL, which was observed to induce a five-fold enhancement compared to the input/source flux and a 3.5 times increase (with a center heat flux: Q c ∼ 2460 W/m 2 ) compared to the PTL (with Q c ∼ 704 W/m 2 ). It was again verified, for both the lenses, that the Q averaged over the cross-section of the sample was ∼ 502 W/m 2 and ∼ 505 W/m 2 , respectively -in excellent accord with the flux through the preceding acrylic block (of ∼ 500 W/m 2 ), and as also computed through the ratio of the temperature drop across the acrylic of ∼ 75 K (with a T hot = 380 K) to the thermal resistance: L ac (=3 cm)/κ ac (=0.2 W/mK). The heat losses due to convection and radiation were estimated, through simulations, to be < 5% of the flux and have been ignored. Through ensuring energy conservation through the composite in the steady state, such Q values were intermediate between considering that the individual layers act as parallel or serial thermal resistors, falling within the respective upper and lower Wiener bounds 74 and signifying the ETM-like character 22 of the thermal lenses. It was evident that the dual use of varying the effective Q traversal length as well as the κ results in enhanced thermal energy concentration. Given the typical correlation of diffusive processes operating under gradients to Gaussian distributions, 75,76 the thermal energy concentration could be modeled appropriately.
In summary, this study has indicated the development of thermal lenses based on a thermal resistance minimization principle, leading to the concentration of heat flux. The subject area could also be further developed through analogies to concepts extant for diffusive light. 76 However, many issues yet remain to be resolved including (a) minimum spot size for the focused thermal energy, which could be related to a mean phonon wavelength, in analogy with the Rayleigh criterion 53 used in optics. In this context, the mean free path is the minimum distance over which a temperature gradient (and a heat flux) could be defined, 32 and may also be considered a practical limit for heat focusing. Moreover, at low heat flux the efficiency of focusing would depend on the convective and radiative heat transfer losses. While the convective loss could be drastically reduced, e.g., if the lenses were placed in a vacuum-like environment, the radiative losses could be substantial. Consequently, adequate radiation shielding should be provided for high efficiencies. The efficiency may also be improved through alternate materials choice and arrangement as well, e.g., a large thermal conductivity contrast ratio (κ center /κ edge ) would also indicate a greater concentration of the heat flux, following the discussion in Sections III D and III E, and further improve the flux concentration efficiency. 24

IV. ENERGY HARVESTING IN METAMATERIALS
We have, so far in the paper, focused on how thermal energy orientation and harness could be accomplished through the use of layered materials. The regulation of the temperature gradient that is evident through such metamaterials oriented design also yields ideas for other types of energy transduction mechanisms. At the very outset, the flow of heat current across a material has ben known to induce a voltage (through the Seebeck effect 77 ) -the magnitude of which is inversely proportional to the electrical carrier concentration in the material. There are also concomitant thermoelectric and electrothermal effects 78 that can be manifested in metamaterials. Such phenomena could be used for electrical voltage and power generation through thermal energy harvesting. 79,80 We give a brief introduction to such phenomena, with close reference to the fundamental formulation previously discussed in Section II A -II D.

A. Thermoelectric effects: The anisotropy of the Seebeck coefficient tensor
As an example of the relatively unique attributes of metamaterials, we discuss the possibility of tuning the Seebeck coefficient (S), through positive, zero, and negative values, through the arrangement of the constituent media in metamaterials. The S has been defined as the ratio of the electrostatic potential difference to the temperature gradient, in the absence of electrical current. 81 It has been shown that the magnitude of the coefficient, |S| varies inversely with the reduced Fermi energy, η (= E F /k B T), where E F is the Fermi energy, k B is the Boltzmann constant and T is the temperature. 82 The implication is that a small (/large) S is obtained for a metallic (/insulating) substance. 81 For example, the S values are of the order of 1-10 µV/K for metals, while they range from 100 µV/K for degenerate semiconductors to 1 mV/K for non-degenerate semiconductors. The sign of S depends on whether electrons or holes constitute the current, in which case the S < 0, or S > 0. For superconductors, the S has been defined to be zero. 81 However, a major issue in the traditional usage of the S is that it is typically considered to be a scalar quantity and ignores the tensorial characteristics, which relates the temperature gradient ∇T j to the electric field, E i through: The Seebeck coefficient tensor (S ij ) is a relation between the i th -component of the electric field to the j th component of the temperature gradient vectors, quite similar in form to Eqn. (1). For a rectangular coordinate system, e.g., i, j = x, y, and z, the full form, cf. Eqn.
(2), would be: S xx S xy S xz S yx S yy S yz S zx S zy S zz (15) When the S ij is equated to a single scalar value, it is being assumed that the cross-diagonal terms are all equal to zero. However, it would be interesting to understand whether S ij , with i j,would be feasible and the resultant consequences. Such possibilities are evident through the consideration of the anisotropy of the S ij in metamaterial architectures. For example, it is implicit from Eqns. (14) and (15), that with a given temperature gradient, say in the horizontal (x-) direction, off-diagonal values for the S ij imply cross-components of the electric field/voltage, i.e., a finite S x y indicates voltages in the x-direction as well as in the perpendicular y-direction. It can also be shown that arbitrary values of the S may be obtained just through arrangement and orientation of thin layers in a given composite.
We first show that S ij is invariant to coordinate transformation. The implication is that for a geometric transformation, say, from one coordinate system, [x 1 , x 2 , x 3 ] to another, say [x 1 in rectangular coordinates, Eqn. (14) would be modified to: ) denotes the Seebeck coefficient tensor in the transformed coordinate system. J is the Jacobian matrix of the coordinate transformation between the new and the old coordinate systems, J T is the transpose of J, and det (J) denotes the determinant.
We indicate the principles of such a transformation through application on a multilayered system, consisting of two materials: see Figure 1(a), with isotropic values of the Seebeck coefficient (S 1 and S 2 , respectively), electrical conductivity (σ 1 and σ 2 , respectively), as well as the thermal conductivity (κ 1 and κ 2 , respectively). Consider a horizontal temperature gradient (say, in the x-direction) applied to the composite depicted in the figure. Under one-dimensional parallel and perpendicular heat transport, between the surface on the left (maintained at a temperature: T hot ) and the right surface (at a lower temperature: T cold ), we can derive (in a manner similar to that reported earlier for Eqn. (3): At the very outset, this is a remarkable relation in that it involves the parameters S and κ for the longitudinal effects while S and σ are involved in the transverse dimensions. The rationale for the specific form of S x is that the temperature drops (and the consequent measured voltages, from V = S ∆T) are additive across the two materials and are related to the respective thermal conductivity values. Alternately, for the S y /S z , a single temperature gradient (and a voltage) is obtained. While it was previously noted that the Seebeck coefficient was defined in the absence of an externally measured current, the definition does not preclude internal currents. Indeed, for a semiconductor with two types of carriers, i.e., the electrons (e: in the conduction band) and holes (h: in the valence band), 83 the overall S = S e σ e +S h σ h σ e +σ h similar in form to the S yy and the S z z components of Eqn. (17). The electronic and hole currents circulate in the semiconductor and equilibrate to yield a net current of zero. A similar rationale may be used for the case of a parallel arrangement of individual layers, where addition of the internal currents (= V/R) through layers 1 and 2, constituting the composite.
While the weighted averaging ensures that the net S in all the three dimensions lies in between the isotropic Seebeck coefficients of the constituent materials, it would be possible to obtain any designed intermediate value through the appropriate use of the layers' σ and κ. Additionally, we observe that the S may be zero: (1) the conditions for ensuring a zero S x are that κ 1 κ 2 = − S 1 S 2 while S y /S z can both tend to zero when σ 1 σ 2 = − S 2 S 1 , in which case the S x would be very large), (2) for S 1 and S 2 of opposite signs, S y /S z > S x . As previously noted, the S has been defined to be zero for superconductors. Metamaterials offer the possibility of a zero S at room temperature! Now, when the layer orientation in the composite is rotated, in the plane, say, by an angle θ (−π/2 < θ < π/2, with θ considered positive in the counter-clockwise direction) the layers would now be oriented along new axes: x ′ (= x cos θ + y sin θ) and y ′ (= −x sin θ + ycos θ); see Figure  2, for the definitions. Such a geometrical transformation is identical to what has been indicated earlier through Eqn. (6) for understanding the re-direction of the heat flux vector. The originally applied (horizontal) temperature gradient now acquires both horizontal and vertical components, with respect to the rotated layers. The degree of anisotropy would change from one composite sample rotation to another and the consequent change may be viewed in terms of tuning the S ij . The resulting Seebeck coefficient tensor of the composite would be given by: The modified electric field would then be: An immediate implication of the above relation would be that depending on the angle of rotation (θ), the electric fields and the measured voltages could be positive, zero, or negative.

B. Experimental investigations
Such off-diagonal Seebeck related thermopower (through the ratio of the transverse voltage generated in response to a longitudinal temperature gradient/heat current) was measured in "artificially created tilted metallic multilayer structures" constituted from sintered alternating constantanchromel foils, where the heat was produced through laser irradiation. 84 A large thermoelectric anisotropy (related to the difference of the parallel and the perpendicular values of the S) of ∼ 1.5 mV/K, was also noted 85 through the use of an Al/Si layered system. It is remarkable that such a large effect was reported in this composite, at room temperature, given that the S values are of the order of 1-10 µV/K for metals, and around 100 µV/K for semiconductors, and may indicate problems with measurement and/or interpretation. A concomitant Peltier effect (an electrothermal phenomena, 81 where electrical current could cause a cooling in the transverse direction) was also considered in an alternately layered Pb/Bi 2 Te 3 composite with a recorded temperature difference of ∼ 22 K, transverse to the direction of the applied electrical currents. 86 Thickness ratio (of the Pb thickness to that of the Bi 2 Te 3 ) dependence was noted for the S anisotropy (parallel vs. perpendicular) with a maximum for a ratio close to unity.
Another relatively recent effort 87 in the direction of energy harvesting, concerns the use of such layered metamaterial-like features to enhance the thermoelectric figure of merit 88,89 (ZT = S 2 σ /κ). Such new methods are to be greatly welcomed given that most recent progress in the ZT seems to have been achieved 82 at the expense of reducing the quality of the material through diminishing the κ, with little attention being paid to the power factor (= S 2 σ). 83,90 The broad idea involves simultaneously regulating the (i) thermal conductivity, (ii) electrical conductivity, and the (iii) Seebeck coefficient, through alternating thermoelectric (TE) and metal layers, and inducing thermal and electrical anisotropies provides a new methodology for tuning composite thermoelectric characteristics. The principles mainly involve the relations in Eqns. (3) and (17) for the conductivities (thermal or electrical) and the Seebeck coefficients, respectively. The experimental arrangement (comprised of metallic Ni layers, with isotropic values of S ∼ −20 µV/K, ρ ∼ 17 µΩ − cm, and κ ∼ 51 W/mK and TE Bi 0,5 Sb 1.5 Te 3 layers, with S ∼ 210 µV/K, ρ ∼ 1.2 mΩ − cm, and κ ∼ 1.1 W/mK) is indicated in Figure 13(a). It is to be noted, at the very outset, that the electrical field gradients (in the horizontal direction) and the temperature gradient (in the vertical direction) are orthogonal.
When the layer orientation is rotated by an angle (θ), we have seen using the precepts previously employed for the thermal conductivity (e.g., Figures 2 and 3) that there would be a considerable modulation of the composite properties. From an intuitive point, for the case of the Ni/TE composite under discussion, a low value of the θ would favor passage of the electrical current (J e ) along the metallic layers (implying a large σ) while presenting a relatively large thermal impedance (implying a low κ) : Figure 13(b). Such variations should of course be balanced with the constituent change of the S with respect to the overall ZT, per the relations in Eqns. (17) and (18). It is interesting to now consider, from Figure 14(a), that the ZT as well as the S are now a function of the angle of orientation of the layers in the thermal metamaterial. This is in much distinction to the usual mode of thought where these parameters are typically considered to be functions of the carrier concentration. 82 Additionally, from Eqn. (19), it may be inferred that the electric field orientation and the sign of the voltage may also be changed due to varying θ. It should be noted that, while Figure 14(a) indicates that the S is always positive and would be maximized at θ = 45 • (per S xz = (−S x + S y ) cos θ sin θ, according to the off-diagonal term in Eqn. (18)), it may be inferred -also from Eqn. (18) that the S may be negative or zero as well, depending on the choice of the θ. As discussed previously, the notion of S = 0 is quite novel, and may be used to reduce parasitic thermoelectric voltages that often creep up due to subtle temperature differences.
Preliminary experimental verification of the above principles (with a Ni/BST thickness ratio of unity) was done in samples fabricated by the SPS (Spark plasma sintering) methodology and seems to accord well with the theoretical estimates: Figure 14(b), for the generated power (right axis) and the voltage (right axis). An increase of the σ concomitant with a decreased κ (by as much as 90%) was reported. Additional modifications, such as the (i) selective introduction of porosity into either of the layers to either increase or decrease the figure of merit, and (ii) the power generation (where the individual layers act in concert to serve as a solid state heat exchanger) could both substantially enhance the utility of such layered thermal metamaterials.
In summary, it was interesting to consider that the layered arrangement of materials can yield desired voltages and electric fields in addition to thermal energy orientation. While the efficiency of such schemes still remains to be explored, e.g., the averaging procedure (either with the layers in parallel or in the perpendicular orientation) would imply that only a mean of the constituent  Figure 12). A low (/high) θ corresponds to electrical conduction through the metal (/thermoelectric) respectively. The xz orientation implies a diagonal orientation: see Figure 12 properties would be obtained, the notion of multi-functionality of thermal metamaterials is very attractive and could yield new applications. It is also important to note that losses due to convective and radiative transfer have not been adequately considered and may pose a problem for low temperature heat recovery. While interfaces involved in joining the two layers, have been mostly neglected in the ETM approximation (Section II C) implicitly assumed, they may yet be important in determining the reproducibility of the observed metamaterial characteristics due to the inevitable variability of the interfacial roughness and constituent thermal conductivity.

V. FRONTIERS IN LAYERED THERMAL METAMATERIALS
The relevance of metamaterials to thermal energy re-orientation and focusing as well as the possibilities of thermal energy harvesting has only involved very preliminary investigations. The notion of layered arrangement has served as the basis and the details of the efficiency of such arrangements have to be much elucidated. For instance, the validity of even fundamental considerations, such as the ETM, need to be investigated for designating a thermal metamaterial. In contrast to electromagnetic metamaterials, where an order of magnitude difference in the wavelength: structure ratio is taken as a criterion for an effective medium, 48 the lack of a characteristic wavelength in thermal conduction makes such a criteria more difficult to define in the thermal case. Such a foundation is important in that interfacial effects are often considered to dominate heat conduction processes between disparate materials. We discuss the influence of such possible effects and speculate on future directions for unifying heat transfer phenomena (i.e., conduction and radiation) based on metamaterials.

A. The role of interfaces: Probing the validity of the ETM
Interfaces may play a major role in modulating heat transfer 91,29 in thermal metamaterials. At the very outset, the interfacial thermal resistance 46 (R int ) is defined for an area, A, through: R int A = l int /κ int , where l int and κ int are the effective thickness and thermal conductivity of the interface, respectively. In practical construction of a multilayered metamaterial/composite, such as that considered in this paper, there is an intrinsic variability in both the l int and the κ int , and could play a major role in establishing a true metric of the thermal boundary resistance (TBR). 47 Generally, there could be a substantial variation of the heat flux bending, due to such physical parameter (i.e., the l int and κ int ) fluctuations intrinsic to an interface. This may be understood through the joining of Al layers (κ 1 = κ Al ) together with epoxy (κ 2 = κ int ), which serves as the interfacial (int) material, following the concepts laid out in Section II. When both l int and κ int are considered the heat flux bending angle would be modified to φ int ( see cf. Eqn. (8), to a value: The change in the φ int (= dφ int ) due to the variation in the l int (= ∂l int ) and the κ int (= ∂κ int ) may be understood through a relation of the form: Consequently, the experimental observation of the heat flux could yield insights into the interfacial characteristics and has been implemented, through a non-contact methodology. 25 A higher κ Al /κ int ratio contributes to the smaller deviation of the heat flux bending angle. It was pertinent that the interfacial conductivity variation may play a larger role compared to the thickness variation. As there would always be such thickness and thermal conductivity variation, the degree to which a thermal metamaterial can be approximated as an effective thermal medium (ETM) should be carefully considered. One would expect that in an ideal thermal metamaterial, that the constituent TBR should be zero/very low.

B. Transduction of the conducted heat: The interplay of conduction and radiation
Initial metamaterials related efforts applied to the thermal regime, incorporated designs directly from those initially postulated for electromagnetics 92 -with fabricated structures of a particular absorption/emission resonance frequency in the infrared regime 93 corresponding to thermal wavelengths. In contrast, the layered materials approach discussed in this paper provides an alternative take through g non-resonant structures, as appropriate to diffusive transport, and conductive heat flux in an effective thermal medium. However, it would be relevant to relate in more detail the heat conduction mechanisms from the atomic scale upwards.
In this context, thermal conduction through the agency of phonons has been extensively investigated both from a theoretical and experimental point of view in the bulk as well as nanostructures. 39,[94][95][96][97][98][99] It would be desirable to understand a correspondence between the phonons and radiated photons in order to trace a pathway of heat from the source to a sink. Moreover, the photons may then be put to useful purpose, e.g., for photovoltaic or thermo-photovoltaic (TPV) energy harness. 79,100-103 A mode of thought may consider a quasi-classical argument that a given temperature corresponds to an equivalence, through k B T, of the phonon and photon energy, i.e., hc/λ photon = hv/λ phonon , where h is the Planck constant, c is the velocity of photon/EM wave transport (=3 X 10 8 m/s), v is the thermal transport velocity and λ photon and λ phonon are the respective photon and phonon wavelengths. Typically, v is of the order of 10 3 -10 4 m/s. From a Brillouin zone construction, it is well known that the range of wave-vectors ranges in between zero and ∼ 10 10 m −1 , implying a λ phonon range from infinity to ∼ 10 −10 m (the assumption being that the atomic spacing is of the order of 0.1 nm). Correspondingly, the range would be ∼(6-60) µm < λ photon < ∞, which marks the range of the thermal wavelengths. It would then be interesting to consider the correspondence between the phonons to the emitted radiation. Moreover, the role of the optical phonons (which have a much smaller group velocity compared to the acoustic phonons) in such phonon-photon coupling is far from clear. The coupling of long-wavelength EM waves to optical phonons has been implicated in the formation of polaritons (SPhPs), 58 which may serve the role of resonant excitations and have yielded evidence for coherent thermal radiation very close to the surface. 68 Such studies may also shed further light on the possibilities of enhanced radiative transfer between closely spaced layers or bodies. 51,104

VI. CONCLUSIONS
Much of the initial work in thermal metamaterials for heat flux control had been motivated 18 from principles formulated for electromagnetic waves and subsequently adapted to heat transport, e.g., through adaptation of geometric and transformation 18 optics based approaches. 6,20,105 However, over the past two years, the field has developed alternative methodologies where metamaterials could be used, e.g., for thermal cloaking, 7,8 but which may be yet based on Maxwell equations. 106 It was the aim of the present paper to provide an alternate viewpoint from the consideration that the simplest composite/metamaterial could consist of two constituents, e.g., arranged as layers: Figure 1. With this basis and on the foundation of the Fourier law of heat conduction, various intriguing physical phenomena that could be pertinent to heat flux conduction, such as (a) anisotropy and re-directing of heat flux in a direction other than in which the temperature gradient was applied: Section II, (b) positive and negative refraction of the heat flux; Section III A, (c) the stoppage of heat due to the concomitant anisotropic diffusivity: Section III B, (d) the possibility of thermal lenses: Section III C, and (d) thermal energy harvesting through metamaterial based thermoelectric effects: Section IV, were all probed through analytical and computational methods and experimentally verified in many cases.
The possibility that heat flux can be controlled and transduced to useful energy (e.g., through coherent thermal radiation) provides a large imperative for the continued study of such phenomena.
For the further progress and vitality of the field of thermal metamaterials, there is a continual need for new insights as well as for the re-examination (and renunciation) of the earlier concepts. Given the rapid progress of scientific investigations in the present day, there seems to be an urgent need to demonstrate unique functionalities, such as the extremely efficient concentration or dissipation of thermal energy as well as to determine the limits of such functionalities. The merger of thermal conduction with radiation, in the process achieving an insight into the relationship between diffusive and wave-like transport, would be a worthwhile outcome of such investigations with the potential of yielding hitherto unsought of applications.