Independent Manipulation of Heat and Electrical Current via Bifunctional Metamaterials

Spatial tailoring of the material constitutive properties is a well-known strategy to mold the local flow of given observables in different physical domains. Coordinate-transformation-based methods (e.g., transformation optics) offer a powerful and systematic approach to design anisotropic, spatially-inhomogeneous artificial materials ("metamaterials") capable of precisely manipulating wave-based (electromagnetic, acoustic, elastic) as well as diffusion-based (heat) phenomena in a desired fashion. However versatile these approaches have been, most designs have so far been limited to serving single-target functionalities in a given physical domain. Here we present a step towards a"transformation multiphysics"framework that allows independent and simultaneous manipulation of multiple physical phenomena. As a proof of principle of this new scheme, we design and synthesize (in terms of realistic material constituents) a metamaterial shell that simultaneously behaves as a thermal concentrator and an electrical"invisibility cloak". Our numerical results open up intriguing possibilities in the largely unexplored phase space of multi-functional metadevices, with a wide variety of potential applications to electrical, magnetic, acoustic, and thermal scenarios.


I. INTRODUCTION
Traditionally, conventional materials have been devised and engineered to serve only singletarget applications. In an integrated circuit, for example, each component is designed to play a specific role: metallic interconnection lines carry electric currents, while a separated block works as a heat sink for dissipating heat. If a single building block could be designed to perform multiple functions in different physical domains, independently but at the same time, this could lead to a completely new way to design complex systems. Natural media are not conceived to accomplish multiple functionalities at the same time, and for this reason taming different physical phenomena at will is a tough proposition.
A new avenue could be paved with the employment of properly engineered artificial materials. Driven by the ability to induce physical responses absent in Nature, the field of "metamaterials" has seen a tremendous growth in recent years. One of the catalysts for the progress made in this field, theoretically as well as experimentally, has been the so-called "transformation optics" theory [1,2]. Viewing the rerouting of energy flow as a distortion of space from a coordinate transformation, the correspondence between constitutive material parameters and geometric transformations can serve as a powerful recipe for designing and fabricating artificial structures. This approach has been utilized not only for the manipulation of electromagnetic waves [3][4][5], but also for acoustics [6][7][8][9][10], elastodynamics [11][12][13][14][15][16], electrostatic [17][18][19] and magnetostatic [20][21][22][23][24] fields, as well as liquid surface waves [25] and diffusive heat flow [26][27][28]. The reader is referred to [29][30][31] for recent perspectives and reviews of metamaterial applications to diverse fields. Within this framework, also worth of mention are some recent multiphysics studies aimed at exciting surface plasmon polaritons in graphene via the interplay of light and sound waves [32,33].
From the mounting experimental applications in various physical branches, it is clear that the strength of the transformation-optics theory is first and foremost its unconventional versatility.
Taking advantage of it, one may envision applying the theory to simultaneously manipulate multiple physical phenomena in independent fashions. For example, a material may be designed to exhibit a particular thermal functionality while its electrical functionality is made drastically different via separate but intertwined coordinate transformations.
Through the example of designing a metamaterial shell that behaves as a thermal concentrator and an electrical "invisibility cloak" at the same time, we present here a framework that allows access to the phase space of multi-functionality with metastructures. Utilizing coordinate transformations while effectively linking phenomena in multiple physical domains, we demonstrate a step towards a general platform that can be called the "transformation multiphysics".
Accordingly, the rest of the paper is organized as follows. In Sec. II, with specific reference to the thermal and electrical scenarios, we outline the modeling aspects pertaining to the transformation media, their effective-medium implementation, and the numerical simulations (with details relegated in Appendices A and B). In Sec. III, we illustrate a proof-of-principle example of synthesis in terms of realistic material constituents. In Sec. IV, we provide further insight into the response exhibited by our metastructure, as well as some bounds dictated by practical limitations. Finally, in Sec. V, we provide some brief conclusions and perspectives.

A. Thermal and Electrical Transformation Media
Although, in principle, our approach could be applied to different physical domains, our focus here is on the thermal and electrical responses. As illustrated in Fig. 1(a), we begin by considering an auxiliary space ′ r = ′ x , ′ y , ′ z ( ), filled with an isotropic medium of thermal and electrical conductivities ′ κ and ′ σ . At equilibrium, the stationary heat and electrical conduction equations in the absence of sources are given by with ′ T and ′ V denoting the temperature and electrical potential, respectively. In the homogeneous case (i.e., ′ κ and ′ σ constant), if temperature and potential differences exist at the two boundaries, the heat-flux and electrical current-density would be directed along straight, parallel paths, as schematically depicted in Fig. 1(a). This is the typical behavior of natural materials.
Next, we introduce two coordinate transformations to a new curved-coordinate space r , namely, with the subscripts "t" and "e" denoting the thermal and electrical domains, which induce different local metric distortions in the two physical domains. For instance, as shown in Fig. 1(a), we consider a concentrator-type transformation in the thermal domain, and an invisibility-cloaktype transformation in the electrical domain. By exploiting the form-invariance properties of Eqs.
(1), the temperature and potential distributions in the transformed domains can be readily related to the original quantities as [31]: Moreover, the distortion effects induced by the coordinate transformations can be equivalently obtained in a flat, Cartesian space r = x, y, z ( ) filled with an inhomogeneous, anisotropic "transformation medium" [cf. Fig. 1 with denoting the Jacobian matrices associated with the two coordinate transformations, and the superscripts " −1" and " −T " indicating the inverse and inverse-transpose, respectively. In such a medium, the heat-flux and electrical current-density would follow markedly different paths. For instance, in the concentrator/cloak example chosen, the heat-flux would tend to concentrate in the inner region, whereas the current-density would tend to circumvent that region, as schematically depicted in Fig. 1(b).

B. Effective-Medium Modeling and Synthesis
Although it is generally impossible to find a natural material exhibiting the desired constitutive relationships in Eqs. (4), these can be approximated to a certain extent by means of metamaterials. Results available in the literature [18,27,28] deal with the design of a single functionality (e.g., cloak or concentrator) in a single domain (e.g., thermal or electrical), and the only example of bifunctional device implements the same functionality in both thermal and electrical domains [34]. Here, the task requires us to prescribe different functionalities in multiple domains, and we proceed by following a synthesis approach based on the mixture of N different types of material inclusions embedded in a host medium [ Fig. 1(c)]. The host and inclusions are characterized by their thermal and electrical conductivities κ n and σ n , respectively, and filling fractions f n , with n = 0,1,…, N , with the subscript "0" denoting the host medium. Each inclusion is also characterized by a depolarization tensor Γ n , which depends on its shape and orientation [35]. We are therefore led to finding the material and structural compound parameters where κ nom and σ nom are the desired nominal constitutive tensors [arising from Eqs. (4)], whereas κ eff and σ eff are the effective constitutive tensors characterizing the mixture, which can be related to the host and inclusion parameters via approximate mixing formulae [35] (see also Appendix A for details). We highlight the nonlinear character of Eqs. (6) (stemming from the mixing formulae), and the coupling between the thermal and electrical domains via the structural compound parameters f and Γ . Moreover, the search space is constrained by the passivity requirements κ n ≥ 0 and σ n ≥ 0 , as well as by the self-consistency conditions 0 < f n < 1 , The synthesis is significantly simplified if the same functionality is required in both domains, as in [34]. In this case, F t = F e and [from Eqs. (4)] κ nom ′ κ = σ nom ′ σ , which implies that the problems in Eqs. (6) are decoupled, with only one synthesis needed. Here, the scenario is more complex, and to induce two distinct functionalities we exploit the concept of "neutral" inclusions from the theory of composites [36,37], i.e., inclusions that are matched with the host medium in one physical domain, so that they are effectively "visible" only in the other domain. This assumption too decouples the thermal and electrical syntheses in Eqs. (6), but it does not constrain the two functionalities to be identical. Clearly, working with natural material constituents, the required neutrality conditions may only be fulfilled approximately.
Nevertheless, in principle, such inclusions may be properly engineered via multilayered composites, e.g., along the lines of [38,39].
Focusing on a two-dimensional scenario in the associated ρ,φ, z ( ) cylindrical coordinate system, and letting κ ρ,nom , κ φ,nom , σ ρ,nom , σ φ,nom the nominal values of relevant constitutivetensor components to be synthesized, and κ ρ ,eff , κ φ,eff , σ ρ ,eff , σ φ,eff the corresponding effective parameters, we consider a three-phase mixture featuring two types of elliptic cylindrical inclusions with axes locally oriented along the cylindrical coordinates ρ and φ . The search parameter space therefore comprises the constitutive parameters κ = κ 0 ,κ 1 ,κ 2 { } and , and relevant depolarization-tensor components Γ 1ρ = 1− Γ 1φ and Γ 2 ρ = 1− Γ 2φ . These latter components, for an elliptical inclusion with axes A ρ (along the ρ direction) and A φ (along the φ direction), are given by [35] Assuming that the type-1 inclusions are thermally neutral κ 1 = κ 0 ( ) and the type-2 inclusions are electrically neutral σ 2 = σ 0 ( ), and considering standard Maxwell-Garnett mixing formulae [35], the effective parameters can be written as follows (see Appendix A for details): inclusions, whereas the electrical parameters in Eqs. (10) and (11) do not depend on the type-2 inclusions. Under these conditions, the synthesis problem can be solved analytically in closed form. Referring to Appendix A for the general solution, we consider here the limiting case κ 2 κ 0 and σ 1 σ 0 , which yields the particularly simple results where σ nom = σ φ,nom σ ρ,nom , κ nom = κ φ,nom κ ρ,nom , and the filling fractions f 1 and f 2 appear as free parameters. It can be readily verified that the results in Eqs. (12) and (13) are inherently feasible, as they yield κ 0 ≥ 0 , σ 0 ≥ 0 , 0 < Γ 1φ < 1, and 0 < Γ 2φ < 1, for arbitrary values of the nominal anisotropy ratios κ nom and σ nom . However, practical considerations (related to the spatial arrangement of the inclusions) as well as model-consistency issues effectively restrict the attainable anisotropy ratios to moderate values. These aspects are discussed in more detail in Sec. IV-B below.

C. Numerical Modeling
All numerical simulations of the thermal and electrical responses in our study are carried out by means of COMSOL Multiphysics 4.2, a finite-element-based commercial software package that allows multiphysics simulations in the presence of anisotropic, inhomogeneous constitutive parameters [40]. In particular, for our simulations, we utilize the "Heat Transfer" and "AC/DC" modules [40] in order to solve the stationary, sourceless heat and electrical conduction More specifically, the magnitudes of these (vector) observables are represented in false-color scale, while their local directions are indicated by the superimposed streamlines.

A. Thermal Concentrator and Electrical Cloak
The above synthesis procedure can be applied to the scenario illustrated in Fig. 1 by introducing two (scalar) radial coordinate transformations for which Eqs. (4) can be particularized in terms of the relevant components [34] with the overdot denoting differentiation with respect to the argument. As schematically illustrated in Fig. 3(a), the transformations in the thermal and electrical domains map an annular cylinder of radii R 1 = 2cm and R 2 = 12cm in the transformed space r onto an annular cylinder of radii R c > R 1 and R 2 and a cylinder of radius R 2 , respectively, in the auxiliary space ′ r . From the functional viewpoint, F t yields a concentration effect (with c = R c R 1 > 1 denoting the concentration factor), whereas F e yields an invisibility cloaking effect. In order to achieve these effects, only the boundary values are prescribed, whereas the function behaviors in between are only partially constrained (see Appendix B for more details). In our example below, we exploited this degree of freedom by selecting the two mapping functions so that Though not strictly necessary, this choice allows us to utilize two types of inclusions with identical shape (just rotated of 90°) and filling fractions, which arguably facilitates their spatial arrangement (see Sec. III-B and Appendix B for details).

B. Preliminary Ideal-Parameter Metamaterial Synthesis
Starting from the six-layer nominal profiles in Fig. 3(b), for each sampled value, we extract the scaled conductivities κ ρ ′ κ , κ φ ′ κ , σ ρ ′ σ and σ φ ′ σ (with ′ κ and ′ σ denoting the background parameters), and compute the nominal anisotropy ratios κ nom and σ nom . Next, we choose the filling fractions f 1 = f 2 taking into account the aforementioned assumptions and limitations (see also Appendices A and B). As a rule of thumb, taking into account that the transformation media to synthesize tend to become more isotropic towards the exterior layers (see Appendix B), we assume a gradually decreasing law (from interior to exterior layers) fulfilling the bound f 1 = f 2 ≤ 0.2 (see also our discussion in Sec. IV-B below). Assuming neutral inclusions (κ 1 = κ 0 and σ 2 = σ 0 ) with κ 2 κ 0 and σ 1 σ 0 (assumed, for simplicity, κ 2 = 0 and σ 1 = 0 ), we now have all the entries in Eqs. (12) and (13) to compute the unknown Table I shows, for each layer, the computed parameters. As anticipated (see also Appendix B), we observe that, in view of the particular choice in Eq. (17), i.e., the two types of inclusions have identical shape (just rotated of 90°). As a consequence, from Eqs. (9) and (10), we also obtain that κ 0 ′ κ = σ 0 ′ σ . Assuming elliptical inclusions, also shown in Table I are the axis ratios calculated from the depolarization factors [35], Table I provides the geometrical, structural, and constitutive parameters for the host medium and the two types of inclusions in each layer, which constitutes all the information needed for an inclusion-based implementation.
Based on this information, we generated the geometry in Fig. 4  In particular, in the inner region, the enforced heat-flux is enhanced by a factor 1.51 (concentrator), while the current-density is reduced by a factor 0.55 (cloak).

C. Realistic-Parameter Metamaterial Synthesis
The above design is idealized in the sense that assumes the availability of host materials with strictly prescribed constitutive parameters (cf. Table I) and neutral inclusions, which, in practice, may only be somehow approximated.
While maintaining the same geometrical and structural parameters as those in Table I and Figs. 4(a)-4(c), we tried to approximate this ideal-parameter configuration by utilizing only five realistic material constituents (detailed in Table II)

A. Comparison with Conventional Material Shell
To better understand the effects of our bifunctional metamaterial shell, it is insightful to

B. Realistic Anisotropy Bounds
As anticipated in Sec. II-B, although our synthesis procedure in Eqs. (12) and (13) that may be difficult to arrange in a spatially efficient fashion.
The latter requirement (relatively high filling fractions), on the other hand, entails significantly dense mixtures, for which the assumed Maxwell-Garnett mixing formulae may not represent an adequate model [35].
It is therefore important to estimate some realistic bounds on the anisotropy ratios that arise from these limitations, and the consequent constraints in the coordinate-transformations that may be implemented. In our design procedure above, we found that values of the depolarization factors 0.1 ≤ Γ 1,2φ ≤ 0.9 and of the filling fractions f 1,2 ≤ 0.2 usually allow spatially-efficient arrangements of the inclusions, which are also adequately modeled by the Maxwell-Garnett mixing formulae (see Appendix A). For the case of neutral inclusions (κ 1 = κ 0 and σ 2 = σ 0 ) with κ 2 κ 0 and σ 1 σ 0 , Fig. 7(a) shows the anisotropy-ratio κ ρ,eff κ φ,eff or σ ρ,eff σ φ,eff that can be attained from Eqs. (8)-(11) by letting the depolarization factors and filling fractions vary within the above mentioned ranges. We emphasize that the neutral-inclusion assumption decouples the syntheses in the thermal and electrical domains, so that the results in Fig. 7 It makes therefore sense to represent these bounds in the two-dimensional Fig. 7(b). In such space, a given anisotropy ratio corresponds to a straight line passing through the origin, with the slope decreasing with increasing values of the anisotropy-ratio κ ρ,eff κ φ,eff (or σ ρ,eff σ φ,eff ). In particular, the (dashed) bisector represents the identity transformation F ρ ( ) = ρ . Thus, for fixed filling fractions ( , by drawing the lines corresponding to the maximum and minimum anisotropy ratios attainable [extracted from Fig. 7(a)], we obtain an angular sector [cyan shaded in Fig. 7 which contains all the possible combinations between F ρ ( ) and F ρ ( ) ρ that can be implemented within the assumed parameter constraints. To give an idea, also shown in Fig. 7(b) are the curves pertaining to the ideal concentrator (red) and cloak (blue) coordinatetransformations in Fig. 3(a). It can be observed that only a portion of these curves actually falls within the allowed angular sector. The inset shows a magnified view of these portions, with the Γ φ markers corresponding to the discretized samples in Fig. 3(b), which were purposely chosen so as to fall within the allowed region.
The above analysis provides useful indications for synthesizing more general functionalities, different from the concentrator and cloak in the chosen example.

V. CONCLUSIONS AND PERSPECTIVES
The characteristics that we have illustrated in this study are a vivid example of artificial

APPENDIX A: DETAILS ON THE EFFECTIVE-MEDIUM FORMULATION
In our approach, the effective thermal and electrical constitutive parameters of a multiphase mixture composed of N types of inclusions embedded in a host material are modeled via simple Maxwell-Garnett mixing formulae [35], With specific reference to the cylindrical geometry of interest for our study, with inclusions aligned along the local ρ (and φ ) direction [cf. Fig. 4(d)], we obtain for the relevant constitutive-tensor components [35] κ ρ ,eff The effective-medium model above relies on the calculation of the static polarizabilities of the inclusions, assuming that each inclusion is embedded in an infinite host medium. While this may be an acceptable assumption for sparse mixtures, it becomes inaccurate for densely packed inclusions, for which the medium effectively "seen" outside the generic inclusion is different from the host medium. In this latter scenario, more refined models can be applied such as, e.g., the Polder -van Santen mixing formulae, which approximate the "apparent" medium outside the inclusions as something in between the host medium and the effective medium [35]. However, this yields implicit equations that need to be solved numerically. Acknowledging these limitations, in our approach, we restrict the structural parameters of the mixtures so as to avoid the dense-packing conditions, and remain within the range of applicability of the Maxwell-Garnett mixing formulae in Eqs. (A1) and (A2). As anticipated in Sec. II-B (and better quantified in Sec. IV-B), this inherently limits the attainable anisotropy ratios to moderate values.
The three-phase mixture considered in Sec. II-B is the simplest reduction ( N = 2 ) of Eqs.
(A1) and (A2) that still allows the joint and independent synthesis of different functionalities in the thermal and electrical domains. To better understand this aspect, we first consider the simpler two-phase mixture, i.e., with only one type of inclusions embedded in a host medium. By particularizing Eqs. (A1) to this case ( N = 1 ), we obtain for the thermal parameters with κ 1 = κ 1 κ 0 , and the second equalities following from the consistency conditions f 0 + f 1 = 1 and Γ 1ρ + Γ 1φ = 1 . Similarly, for the electrical parameters, we obtain from Eqs. (A2) with σ 1 = σ 1 σ 0 . We note that, for the extreme values Γ 1φ = 0 and Γ 1φ = 1 , Eqs. (A3)-(A6) reduce to the well-known expressions pertaining to radial and angular multilayers [35], respectively, which have been widely utilized to design coordinate-transformation-inspired metamaterial structures implementing single functionalities (e.g., cloak, concentrator) in the thermal, electrical, or magnetic domains [17][18][19][20][21][22][23][24][26][27][28]. However, it can be verified that the mixing formulae in Eqs. (A3)-(A6) do not provide enough degrees of freedom to design different anisotropic behaviors in the thermal and electrical domains. For instance, assuming that the parameters κ 1 , f 0 and Γ 1φ in Eqs. (A3) and (A4) are chosen so as to guarantee that κ ρ,eff κ φ,eff > 1, it will be then impossible to achieve an anisotropy ratio σ ρ,eff σ φ,eff < 1 from Eqs. (A5) and (A6), for any choice (subject to the passivity condition) of the material parameter σ 1 . This can be verified in a rather straightforward fashion for the limit (multilayer) cases Γ 1φ = 0 or Γ 1φ = 1 , and in a more cumbersome fashion (we relied on the "Reduce" symbolicalgebra tool of Mathematica™ [43]) for general values of Γ 1φ . Clearly, it represents a significant curtail of the capabilities to independently manipulate the phenomena in the two physical domains. For instance, it is clear from Fig. 3(b) that the joint synthesis of a thermal concentrator (which requires κ ρ κ φ > 1 ) and an electrical cloak (which requires σ ρ σ φ < 1 ) would not be possible with this type of mixtures. It is worth highlighting that these constraints may be relaxed in the presence of negative-conductivity material constituents. While these materials are not available in Nature, they may be in turn synthesized as metamaterials. For instance, Fang et. al [19] recently demonstrated experimentally the possibility of synthesizing an artificial material exhibiting negative electrical conductivity, by means of active devices together with resistor networks.
An easier way to overcome the above limitations, while maintaining the passivity requirements, entails considering a three-phase mixture, featuring two types of inclusions embedded in a host medium. For this scenario ( N = 2 ), we now obtain from Eqs. (A1) and (A2): While it is now possible, in principle, to achieve different anisotropy ratios in the thermal and electrical domains, the general synthesis in Eqs. (6) remains a formidable task, which, in the general case, can be addressed in a weak fashion, i.e., by minimizing a suitable cost function parameterizing the mismatch of the effective and nominal parameters, and possibly including regularization terms.

APPENDIX B: DETAILS ON THE COORDINATE TRANSFORMATIONS
As mentioned in Sec. III-A, the thermal-concentration and electrical-cloak effects are essentially established by the boundary values of the mapping functions in Eqs. (16), whereas the function behaviors in between are only partially constrained by the continuity requirement (in order to avoid additional boundary conditions) as well as by the nonnegative character of their logarithmic derivatives [ F t F t ≥ 0 , F e F e ≥ 0 , in order to guarantee passivity, cf. Eqs. (15)]. The cloak transformation utilized in our study [blue curve in Fig. 3(a)] belongs to the general class of algebraic transformations which satisfy the required boundary conditions (16), and map a cylinder of radius R 2 in the auxiliary space onto an annular cylinder of radii R 1 and R 2 in the transformed space, thereby creating a "hole" of radius R 1 that admits no image in the auxiliary space. For γ = 1 , Eq. (B1) reduces to the standard (linear) cloak transformation introduced by Pendry et al. [2]. Here, we consider instead which yields Recalling Eqs. (15), this ensures that i.e, that the arising transformation medium tends (for ρ → R 2 ) to an isotropic material matched with the background medium. As mentioned in Sec. III, while not strictly necessary, this assumption simplifies the inclusion-based implementation.
For the concentrator transformation, we exploit the degrees of freedom in the choice of the mapping function, by enforcing the condition in Eq. (17). Looking at Eqs. (12) and (13), it can be observed that with this assumption (i.e., κ nom = 1 σ nom ), together with f 1 = f 2 , we obtain Γ 2φ = 1− Γ 1φ = Γ 1ρ , which means that the two types of inclusions have identical shape (just rotated of 90°). Once again, while not strictly necessary, the above assumption may facilitate the inclusion-based implementation, allowing more efficient packing strategies. Recalling Eqs. (6), the condition in Eq. (17) yields the differential equation with the boundary condition . By substituting Eq. (B1) in Eq. (B5), we then obtain which admits the simple analytical solution considered in our study [red curve in Fig. 3(a)]. From Eq. (B7), we can easily calculate the concentration factor, thereby verifying the concentrator functionality.  Table I. Geometrical, structural, and constitutive parameters of the ideal-parameter synthesis, from the piecewise-constant nominal profiles in Fig. 3(b), assuming neutral inclusions (κ 1 = κ 0 , σ 2 = σ 0 ) with κ 2 = 0 and σ 1 = 0 .  Corresponding relevant constitutive-tensor components κ ρ ′ κ = ′ κ κ φ (red) and σ ρ ′ σ = ′ σ σ φ (blue). Outside the annulus R 1 < ρ < R 2 , the coordinate transformations reduce to the identity,  Table II. represents the allowed region, and the dashed bisector the identity transformation. Also shown are the curves pertaining to the concentrator (red) and cloak (blue) coordinate-transformations in Fig. 3(a). The inset shows a magnified view of the allowed portions of these curves, with the markers corresponding to the discretized samples in Fig. 3(b).