Free-form thermal cloaks in three dimensions

Cloaking oneself from the environment to evade predators or confuse prey is widely seen in nature, such as cuttlefish, octopus, and squid^1,2. Emulating this intriguing bio-function, physicists have developed the celebrated transformation theory^3,4,5,6,7,8, which predicts the existence of artificial thermal cloaks capable of hiding three-dimensional (3D) objects omnidirectionally from thermal detection^6,8. The cloaks, however, require spatially graded and highly anisotropic thermal conductivity, posing a significant challenge for engineering realization^9,10,11. Indeed, most realizations of thermal meta-devices thus far remain two-dimensional (2D), mostly with regular shapes (cf. circles)^{12,13,14,15,16,17,18,19,20,21,22,23}. Extension to general 2D geometries is a recent progress, achieved by inversely designed microstructures and their assembly into thermal meta-device^{9,23,24,25,26}. A recent development proposes integrating analytical rank-2 laminates with a 2D de-homogenization technique, leading to concise and well-connected 2D thermal meta-devices¹⁰.

Many real-world applications are three-dimensional, and how to realize the transformation theory in a general 3D domain and engineer free-form 3D omnidirectional thermal meta-devices remains far from being understood. The fundamental and long-standing challenge is how to create functionally robust, well-connected, and practically fabricable 3D microstructures that deliver the required graded and highly anisotropic 3D graded conductivity^9,10,11. The majority of existing studies in 3D are theoretical and numerical ^{27,28,29,30,31,32,33,34,35,36,37}, most of which assume continuous and physically unrealized materials. A few studies have realized 3D meta-devices with simple geometries (i.e., sphere or ellipsoid) ^38,39,40 or 2D-extruded domains^41,42, where conductive shells or regular laminates may produce the needed properties. A recent attempt to achieve irregular 3D cloaks suggests circumventing the transformation theory and using a numerical conformality tracing technique¹¹. While producing fabricable 3D structures with isotropic properties, the cloaking is not omnidirectional, and achieving directional independence demands strong anisotropy⁴³. Indeed, the structural realization of 3D graded anisotropy can realize the general 3D transformation theory and enable the engineering of wide-ranging 3D omnidirectional thermal metamaterials. In 3D, otherwise popular numerical inverse design widely used in 2D studies (such as topology optimization^44,45) becomes prohibitively costly, results in overly complex microstructures, and would have difficulties in ensuring 3D connectivity.

Here, we physically realize the transformation theory in 3D space and create 3D omnidirectional thermal cloaks with arbitrary internal and external geometries. This is achieved by leveraging special spectral properties of the transformed 3D conductivity and developing a new 3D de-homogenization method^{46,47,48,49,50} with a functionally robust, yet geometrically simple, 3D lattice microstructure. Featuring complex shapes, the created 3D thermal cloaks exhibit concise structures with guaranteed global connectivity, successfully concealing 3D objects with a wide range of conductivities under multi-directional thermal detection. Moreover, the created 3D cloaks can be efficiently fabricated by metal 3D printing and mold casting, and the 3D omnidirectional cloaking is accurately reproduced in experiments. With a proposed integration of 3D de-homogenization with spherical harmonics—the Fourier transform on a sphere—and a transverse isotropy regularization technique, we achieve 3D thermal cloaks with unprecedented geometric complexity, such as the shape of a human face, while still exhibiting excellent thermal cloaking. The findings resolve the fundamental challenge of realizing 3D transformation thermotics, creating a practical route for engineering wide-ranging omnidirectional 3D thermal meta-devices as well as other 3D conductive systems, such as electrical conduction and pressure flow.

3D thermal cloak by de-homogenization

As illustrated in Fig. 1A, we aim to conceal an object (apple-shaped) in a region (pear-shaped) from omnidirectional thermal detection. The apple surface is defined by surfaces R₁(θ, ψ) and the pear by R₂(θ, ψ), as shown in Fig. 1B (expressions for R₁ and R₂ are given in SI 1). The required 3D anisotropic thermal conductivity κ is obtained by the transformation thermotics formula \({\boldsymbol{\kappa }}=\frac{{\kappa }_{0}}{\det ({\bf{J}})}{\bf{J}}{{\bf{J}}}^{{\rm{T}}}\)⁸, where κ₀ is the (isotropic) thermal conductivity of the background, and J is the Jacobian of the coordinate transformation (see detailed expression of κ in SI 2), which depends solely on the domain geometry. This study uses thermal conductive encapsulant (DOW DOWSIL^TM TC-6020) as the background material with κ₀ = 2.72 Wm⁻¹ K⁻¹. The three locally ordered principal values κ_i, i = 1, 2, 3 are shown in Fig. 1C. The magnitudes and distribution patterns of the principal values κ_i differ significantly, demonstrating strong anisotropy and spatial gradient. Although complex, the κ possesses special spectral properties that can be exploited for the subsequent design process (see detailed analysis of κ’s spectral properties in SI 3).

To realize the complex and anisotropic κ, while achieving high efficiency and global connectivity, we develop a 3D de-homogenization strategy that consists of local microstructural realization and global formation of the meta-device. For local microstructure, we propose a simple yet highly robust 3D lattice as shown in Fig. 1D. It is made from two isotropic constituents, Materials A and B, with conductivity κ_A ≥ κ_B. This study uses AlSi10Mg (κ_Al = 113.0 Wm⁻¹ K⁻¹) as Material A and PDMS (κ_PDMS = 0.16 Wm⁻¹ K⁻¹) as Material B. The Material A part consists of three square-cross-section bars characterized by the cross-section widths b_i, i = 1, 2, 3, and the rest of the cell is filled with Material B. The geometrically simple three-parameter microstructure is highly robust as it achieves a wide range of the physically admissible 3D conductivity. As demonstrated in Fig. 1D, the achievable κ_i of the 3D lattice (obtained by numerical homogenization, indicated by green circles) fills most space inside the theoretical bound (“G-closure”)^51,52,53 of two-constituent composites, indicated by the pink surface, with the exception of the extreme points of the material space. The lattice’s coverage of the conductivity parameter space is also much greater than that of some other simple microstructures, such as rank-3 laminates, which achieve more extreme properties but in a smaller parameter range (see SI 4.1, 4.2 and Supplementary Fig. 1 for detailed analysis). The latter makes rank-3 or closed-walled microstructures impractical for cloaking purposes. Such important insight about the achievable parameter space was not revealed previously, even though the same form of lattice microstructure was adopted in a previous numerical study³³. The broad coverage of property space rules out the need for much more cumbersome and computationally expensive free-form numerical inverse design^9,24. We note that the geometric simplicity of the microstructure is not only a great advantage but also a requirement for the later global de-homogenization technique.

The parameter identification of the b_i triplets that produce the required κ_i is carried out based on a precomputed surrogate model of the 3D lattice, where the parameter-property database (b_i → κ_j) is constructed using numerical homogenization (see SI 4.1). The parameter identification results in a three-variable least-squares problem solved at all locations of the cloak in a parallel fashion with low computational cost (see SI 4.3 for details of the parameter identification). The distributions of b_i are qualitatively similar to those of κ_i and also show a strong spatial gradient (see SI 4.3 and Supplementary Fig. 2 for the plot of b_i distributions).

With the microstructural local realization, the next task is to connect the lattice smoothly to form the global structure where the local members’ principal directions are aligned with the eigenvectors of κ (denoted as v_i), which are shown in Fig. 1E. The distributions of v_i are also irregular and spatially graded, which is a major challenge for the smooth alignment of the 3D lattices. To achieve this, we develop a 3D de-homogenization approach that represents the three global lattice groups through wave functions and Boolean operations. The basic idea has been used in producing high-resolution 3D structures with optimized mechanical performance^{46,47,48,49,50}. In this approach, we solve for a pseudo-conformal mapping ϕ in the 3D domain such that its gradients align with the eigenvector fields, i.e.,  ∇ ϕ_i = α_iv_i where α_i is a scalar field. The ϕ_i is obtained by solving a constrained least-square problem¹⁰ (see SI 5 and Supplementary Figs. 3 and 4 for details). The obtained ϕ_i, together with b_i are input to cosine wave functions that generate three mutually perpendicular surface groups (see SI 6 and Supplementary Fig. 5 for the detailed procedure), based on which we apply Boolean operations to obtain the three lattice groups. Finally, taking the union of the lattices leads to the global lattice structure shown in Fig. 1F (minor post-processing is applied at the top and bottom center, see SI 6.1 and Supplementary Fig. 6 for details). The complete procedure for generating the 3D cloaking structure is also illustrated in Supplementary Movie 1, and the computational cost is documented in SI 7. A qualitative comparison of different design approaches is provided in SI 8 and Supplementary Table 1.

The global lattice is naturally well-connected without post-processing, and its relative widths comply with the b_i fields while its orientations align with the v_i fields. As shown in the sectional view in Fig. 1F, the radial members (index 1) are generally much smaller than the other two lattice groups, reflecting the relative magnitudes of κ_i in Fig. 1C. The zoomed-in views demonstrate the high resolution enabled by the 3D de-homogenization approach, which is another advantage over numerical inverse designs. The local asymptotic properties of the 3D structure match those of the local 3D lattice. The structure thus achieves 3D omnidirectional thermal cloaking as illustrated in Fig. 1G and as elaborated in the following discussion.

We post-evaluate the performance of the 3D cloak using finite element analysis (FEA, see SI 9.1 for details). As a reference, the temperature profile of the bare core (no cloaking structure) under Z-direction applied heat is shown in Fig. 1H. The isotherms in the background are non-uniform and curved toward the less conductive core (filled with κ_B), exposing the core to thermal detection. For the cloaking case where the structure covers the core, the temperature distributions under Z- and X-direction (Y direction same as X due to rotational symmetry) are shown in Fig. 1I. The background temperature is mostly uniform as indicated by the equidistant and straight isotherms, and the temperature gradient at the core is diminished, showing successful concealing of the apple-shaped core. Quantitatively, we use a relative temperature difference (RTD) value that measures the relative difference of the background temperature from the reference linear temperature profile, with RTD = 1 corresponding to the bare core case without cloaking (Fig. 1H) and RTD = 0 representing a perfect cloaking (see SI 9.2 for RTD’s definition). The cloak’s RTD values under Z- and X-direction applied heat are 0.037 and 0.079, showing little difference from the reference profile and hence great cloaking performance. More discussions on the FEA results are provided in SI 9.1 and Supplementary Fig. 7. We note that the cloaking performance is insensitive to the core material (the cloaked object), the background material (see SI 10.1, 10.2, and Supplementary Fig. 8), and the presence of a heat source in the background (see SI 10.3 and Supplementary Fig. 9). However, it is influenced by the fluctuation of material properties (see investigation in SI 10.4 and Supplementary Fig. 9). Also, the current microstructure feature size is a balance between the separation of scale and fabricability but achieves great cloaking performance, while remaining practically fabricable, as demonstrated next.

Fabrication and validation

A great strength of the proposed strategy is the relatively concise lattice structures that can be practically fabricable despite their graded and highly anisotropic homogenized properties. We fabricate the aluminum lattice structure using metal 3D printing (Direct Metal Laser Sintering, Proto Labs Inc.) with an imposed minimal feature size limit of 0.5 mm. The PDMS parts and the background encapsulant are fabricated with mold casting (see SI 11 and Supplementary Fig. 10 for detailed fabrication process). The left and middle panels of Fig. 2A show the CAD model and the 3D-printed lattice, respectively, featuring identical geometries with many details. The right panel of Fig. 2A shows the complete specimen with casted PDMS parts and encapsulant background.

The setup of the thermal conduction test is shown in Fig. 2A. The specimen is placed between two aluminum plates for smooth heat paths. The top plate is connected to an electric heat source, and the bottom plate is soaked in iced water. The four side faces are cast with a 2 mm-thick PDMS layer to realize the adiabatic boundary condition. The hot end’s temperature is set to 40 °C. The temperature distribution is recorded using an infrared camera (Telops FAST M100hd). Details of the experiment setup are given in SI 11.

The steady-state temperature distributions (recorded after 60 min of heating) of the front and back faces under vertically and horizontally applied heat are shown in Fig. 2B. In the front face, the isotherms in the background remain largely straight and equidistant, demonstrating the unperturbed temperature field as in the homogeneous medium case. Also, the temperature gradients at the core are diminished, with little heat transfer. On the back face, the isotherms are also straight and uniform, a direct display of the 3D cloaking where the apple-shaped object is concealed from thermal detection. The temperature profiles at various sections (S1–S6) are shown in Fig. 2C together with their numerical counterparts. The profiles in the background sections (S1, S3, S4, and S6) are close to linear, and those crossing the cloak (S2 and S5) show a clear plateau at the core region. The experimental results match well with the numerical simulation, validating the proposed 3D de-homogenization strategy. The design, fabrication, and experiment process of the apple-pear cloak are illustrated in Supplementary Movie 2. We note that the choice of the base materials imposes certain practical requirements on the application conditions, such as the operational temperature at which the material properties are stable and the potential influence of interfacial thermal resistance (see SI 12.1 and 12.2 for investigations). Also, a discussion on the cloaks’ transient performance and characteristic time is provided in SI 12.3.

Geometric versatility

The proposed methodology works well for more general and complex 3D geometry. Figure 3A, B shows five 3D cloaks (Dsgs. 1–5) with diverse geometries, where those of Dsgs. 1–3 are defined by surfaces of revolution (similar to that of the apple-pear cloak), and those of Dsgs. 4 and 5 come from general surfaces (parameters are given in SI 1.1 and 1.2). The lattice structures of Dsgs. 1–3 feature relatively small radial lattice members compared to the other two directions, as the radial κ₁ is the smallest. The cloaking performances of Dsgs. 1–3 are shown in Fig. 3C, where the RTD values in all directions are smaller than 0.08, a significant reduction from the bare case of 1.00. Importantly, the cloaking performance is insensitive to the core material’s conductivity as demonstrated in the top-left plot of Fig. 3A, which shows the RTD of Dsg. 1 with wide-ranging κ_c from κ_PDMS (0.16 Wm⁻¹ K⁻¹) to κ_Al (113 Wm⁻¹ K⁻¹). All RTD values are below 0.08 and show further reduction with a higher κ_c. The cloaking performance is also insensitive to the change of background material, as shown in the bottom left of Fig. 3A, where the background’s conductivity κ₀ is increased from 2.72 Wm⁻¹ K⁻¹ to 12.5 Wm⁻¹ K⁻¹. The rise in κ₀ leads to a decrease in RTD in Z direction and some increase in RTD in X and Y directions, with all smaller than 0.12. The increase of κ₀ also results in a close-to-linear rise of Material A’s usage, as shown in the Material A’s volume fraction (VF) v.s. κ₀ relation in the same plot.

Dsgs. 4 and 5 are defined by more general geometries with the surfaces dependent on both θ and ψ, as shown in the cross-sectional views of Fig. 3B. The orientations of the 3D lattice structures are more complex and irregular than the surface of revolution cases (Dsgs. 1–3) and do not lie in a plane in general. This is seen in the lattice group 2 of Dsg. 4 and lattice group 3 of Dsg. 5, plotted in the inset of Fig. 3B. The orientation of the lattice groups aligns with the associated principal vector fields of κ. The lattices are generated through the union of two iso-surface level sets defined by ϕ_i and b_i as described earlier, which are also plotted in Fig. 3B. Despite the more complex structures, Dsgs. 4 and 5 produce great 3D cloaking performance, with RTD in all three directions smaller than 0.15 as shown in Fig. 3C, outperforming most of the 2D cloaks in one of the latest studies²⁵. For demonstration, we fabricate Dsgs. 1, 2, 4, and 5 using stereolithography (SLA) 3D printing with plastic, as shown in Fig. 3D, accurately reproducing the details.

As further validation, we fabricate and test Dsg. 3, as shown in Fig. 3E, where the small lattice members are well-preserved in the metal 3D-printed specimen. The steady-state experimental temperature distributions under Z- and X-direction applied heat are shown in Fig. 3E, which demonstrates successful 3D cloaking with uniform and straight isotherms in the background and diminished temperature gradient at the core. A more comprehensive analysis of the experimental results is provided in SI 13 and Supplementary Fig. 11.

Human-face 3D thermal cloak

The developed strategy allows us to engineer 3D thermal meta-devices with record-breaking geometric complexity. We demonstrate this via a 3D thermal cloak with the shape of human faces, as illustrated in Fig. 4C. The inner boundary R₁ is a male face with glasses, and the outer boundary R₂ consists of two distinct male faces on opposite sides. Due to the surfaces’ high complexity and the requirement of analytical parameterization for transformation thermotics, we propose using spherical harmonics to fit the surfaces based on discrete point clouds from their STL models. Spherical harmonics can be seen as the Fourier transform on a sphere, which analytically represents a complex closed surface. It is particularly suited for 3D thermal metastructures for their closed-surface domains (see SI 1.3 for detailed procedure). The surfaces fitted by spherical harmonics are shown in Fig. 4C, which preserves most facial details. The complex geometry results in highly irregular and somewhat disordered distributions of v_i, which can cause computational issues in the numerical solution of ϕ_i. To resolve it, we adopt an approximated κ to the original one by exploiting its near transverse isotropy in the v₂–v₃ subspace. The near transverse isotropy is evidenced in Fig. 4A. The distribution of the three principal values demonstrates overall close κ₂ and κ₃, while κ₁ is significantly smaller. The proportion of all locations satisfying (κ₃ − κ₂)/κ₃ ≤ η for a range of relative difference η is shown in the right of Fig. 4A. Approximately 60% of all locations have a relative difference smaller than 5%, and 85% of the domain is smaller than 15%. This observation inspires us to approximate κ with one that is transversely isotropic in the v₂–v₃ subspace. This allows us to modify v₂ and v₃ to a more regular and numerically friendly distribution. Details of the approximation and modifications are provided in SI 5.2.

The resulting three lattice groups from the de-homogenization are shown in Fig. 4B, and their sizes reflect the relative magnitude of κ_i in Fig. 4A. Taking the union of the lattice groups forms the 3D cloaking structure shown in Fig. 4D, where the sectional view by the cut plane is overlapped with the outer and inner surfaces. The lattice structures show high anisotropy and spatial gradient, with orientations aligning with the principal vectors. The SLA 3D-printed samples are shown in Fig. 4E, preserving all the details in the CAD model. Despite the complex geometry, the structure produces satisfactory 3D thermal cloaking performance as demonstrated by the FEA-post evaluated temperature profiles in Fig. 4F, where the external and internal boundaries and the cross-sections of the 3D cloak are overlaid. Under applied heat in all three directions, the isotherms in the background are uniform and curved around the human-face core, inside which the temperature gradient is small. The RTD values under the Z-, X-, and Y-direction heat are 0.229, 0.278, and 0.094, respectively, which are lower than most of the 2D thermal cloaks in ref. ²⁵ despite the much more complex 3D geometry.

This study demonstrates a precise, efficient, and practical design-fabrication strategy for realizing omnidirectional, arbitrary-shaped 3D thermal meta-devices predicted by transformation theory. This is made possible by (1) the discovery of special spectral properties of the transformed 3D conductivity that significantly accelerates the subsequent design process; (2) the proposed geometrically simple 3D lattice microstructure with the required large effective property space; (3) a new 3D de-homogenization technique featuring multi-step Boolean operations tailored to the 3D lattice; and (4) the integration of transformation thermotics with spherical harmonics aided by a proposed transverse isotropy regularization.

The combined method enables the creation of 3D omnidirectional thermal meta-device structures that are locally concise, naturally well-connected, capable of adapting to highly complex domains, and practically fabricable. Regardless of the geometry, the 3D thermal cloaks demonstrate great performance in concealing objects under multi-directional thermal detection, as evidenced by their small RTD indexes. With metal 3D printing, we accurately fabricate the 3D thermal cloak and successfully reproduce the omnidirectional cloaking in experiments.

The findings provide a competitive and general solution to the long-standing challenge of accurately realizing 3D graded and anisotropic material properties. They also advance the field of thermal metamaterial from idealized 2D geometries toward practical applications, where objects have arbitrary and complex 3D geometries. The developed methodology is general and transformable, readily extensible to other 3D thermal and conduction-dominated meta-devices across various disciplines.

Definition of a complex 3D domain

This study adopts irregular 3D domains defined by two closed surfaces via spherical coordinates. We use parametric surfaces for most 3D domains in this study, and use spherical harmonics to represent the complex human-face domain. Detailed definition of the surfaces is given in SI Section 1.

3D transformation thermotics and its spectral properties

The study uses 3D transformation thermotics on complex domains. Its detailed procedures are given in SI Section 2. This study discovers and analyzes several special spectral properties of the transformed 3D anisotropic thermal conductivity tensor, some of which can greatly simplify the design process. The analysis and its implications for design are discussed in detail in SI Section 3.

3D lattice microstructure: property space and parameter identification

This study uses a simple and robust three-parameter 3D lattice defined by the cross-sectional areas of the three mutually perpendicular square-sectioned bars. The lattice’s effective 3D thermal conductivity tensor is obtained via numerical homogenization, based on which its property space is obtained. The study reveals that the lattice’s property space fills most regions inside the theoretical bound. Based on the homogenized properties, the study develops a surrogate model that maps parameters to the principal values of the 3D thermal conductivity, which facilitates the efficient parameter identification. Details of these methods and analysis are provided in SI Section 4.

3D de-homogenization

This study develops a new 3D de-homogenization method that smoothly maps the spatially graded 3D lattice to any global complex 3D domains while effectively preserving their asymptotic properties. The 3D de-homogenization is built upon a 3D pseudo-conformal mapping and a multi-step Boolean operation. Details of these techniques are provided in SI Sections 5 and 6. The overall computational cost of the design process is documented in SI Section 7.

Comparison with other design methods

The comparison of the proposed design method for thermal metamaterial with some other existing strategies is provided in SI Section 8.

3D finite element analysis and parametric studies

This study uses 3D FEA to post-evaluate the 3D cloaks’ performance. Details of the FEA and performance measures are provided in SI Section 9. This study provides some parametric studies on the 3D cloaks’ performance, including their sensitivity to changes in core material, background material, presence of heat source, and material property fluctuations. Details of the parametric studies are given in SI Section 10.

Fabrication and experiment

This study uses metal 3D printing to fabricate the Material A lattice structure (AlSi10Mg) and adopts PDMS (DOW DOWSIL^TM 184) casting with plastic 3D-printed molds to manufacture the Material B parts. For the background, this study uses thermal conductive encapsulant (DOW DOWSIL^TM TC-6020) with mold casting. For the experiment, we use electric heating for the hot end and iced water for the cold end of the specimen. Details of the fabrication and experimental setup are provided in SI Section 11. Additional analysis of the experimental result of the heart-shaped 3D cloak is provided in Section 13.

Practical application considerations

This study provides some discussions on practical application aspects in SI Section 12, which include the range of operational temperature, potential influence of the interfacial thermal resistance, estimate of transient performance, and characteristic time scale.

1. Hanlon, R. Cephalopod dynamic camouflage. Curr. Biol. 17, R400–R404 (2007).

   Article  CAS  PubMed  Google Scholar 

2. Woo, T. et al. The dynamics of pattern matching in camouflaging cuttlefish. Nature 619, 122–128 (2023).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

3. Dolin, L. To the possibility of comparison of three-dimensional electromagnetic systems with nonuniform anisotropic filling. Izv. Vyssh. Uchebn. Zaved. Radiofiz. 4, 964–967 (1961).

   Google Scholar 

4. Pendry, J. B., Schurig, D. & Smith, D. R. Controlling electromagnetic fields. Science 312, 1780–1782 (2006).

   Article  ADS  MathSciNet  CAS  PubMed  Google Scholar 

5. Leonhardt, U. Optical conformal mapping. Science 312, 1777–1780 (2006).

   Article  ADS  MathSciNet  CAS  PubMed  Google Scholar 

6. Fan, C. Z., Gao, Y. & Huang, J. P. Shaped graded materials with an apparent negative thermal conductivity. Appl. Phys. Lett. 92, 251907 (2008).

   Article  ADS  Google Scholar 

7. Gersborg, A. R. & Sigmund, O. Extreme non-linear elasticity and transformation optics. Opt. Express 18, 19020–19031 (2010).

   Article  ADS  PubMed  Google Scholar 

8. Huang, J.-P. Theoretical Thermotics: Transformation Thermotics and Extended Theories for Thermal Metamaterials (Springer, 2020).

9. Sha, W. et al. Robustly printable freeform thermal metamaterials. Nat. Commun. 12, 7228 (2021).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

10. Li, W., Sigmund, O. & Zhang, X. S. Analytical realization of complex thermal meta-devices. Nat. Commun. 15, 5527 (2024).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

11. Xu, L. et al. Free-form and multi-physical metamaterials with forward conformality-assisted tracing. Nat. Comput. Sci. 4, 532–541 (2024).

    Article  PubMed  Google Scholar 

12. Han, T. et al. Experimental demonstration of a bilayer thermal cloak. Phys. Rev. Lett. 112, 054302 (2014).

    Article  ADS  PubMed  Google Scholar 

13. Shen, X., Li, Y., Jiang, C., Ni, Y. & Huang, J. Thermal cloak-concentrator. Appl. Phys. Lett. 109, 031907 (2016).

    Article  ADS  Google Scholar 

14. Qu, Y. et al. Thermal camouflage based on the phase-changing material GST. Light 7, 26 (2018).

    Article  Google Scholar 

15. Hu, R. et al. Illusion thermotics. Adv. Mater. 30, 1707237 (2018).

    Article  ADS  Google Scholar 

16. Han, T. et al. Full-parameter omnidirectional thermal metadevices of anisotropic geometry. Adv. Mater. 30, 1804019 (2018).

    Article  Google Scholar 

17. Hong, S., Shin, S. & Chen, R. An adaptive and wearable thermal camouflage device. Adv. Funct. Mater. 30, 1909788 (2020).

    Article  CAS  Google Scholar 

18. Yang, F., Tian, B., Xu, L. & Huang, J. Experimental demonstration of thermal chameleonlike rotators with transformation-invariant metamaterials. Phys. Rev. Appl. 14, 054024 (2020).

    Article  ADS  CAS  Google Scholar 

19. Zhu, Z. et al. Inverse design of rotating metadevice for adaptive thermal cloaking. Int. J. Heat. Mass Transf. 176, 121417 (2021).

    Article  Google Scholar 

20. Ji, Q. et al. Designing thermal energy harvesting devices with natural materials through optimized microstructures. Int. J. Heat. Mass Transf. 169, 120948 (2021).

    Article  Google Scholar 

21. Ji, Q. et al. Deep learning based design of thermal metadevices. Int. J. Heat. Mass Transf. 196, 123149 (6) (2022).

    Article  Google Scholar 

22. Hirasawa, K., Nakami, I., Ooinoue, T., Asaoka, T. & Fujii, G. Experimental demonstration of thermal cloaking metastructures designed by topology optimization. Int. J. Heat. Mass Transf. 194, 123093 (2022).

    Article  CAS  Google Scholar 

23. Zhu, Z. et al. Field-coupling topology design of general transformation multiphysics metamaterials with different functions and arbitrary shapes. Cell Rep. Phys. Sci. 4, 101540 (2023).

    Article  Google Scholar 

24. Sha, W., Xiao, M., Huang, M. & Gao, L. Topology-optimized freeform thermal metamaterials for omnidirectionally cloaking sensors. Mater. Today Phys. 28, 100880 (2022).

    Article  Google Scholar 

25. Wang, Y., Sha, W., Xiao, M., Qiu, C.-W. & Gao, L. Deep-learning-enabled intelligent design of thermal metamaterials. Adv. Mater. 35, 2302387 (2023).

    Article  CAS  Google Scholar 

26. Zhu, Z. et al. Arbitrary-shape transformation multiphysics cloak by topology optimization. Int. J. Heat. Mass Transf. 222, 125205 (2024).

    Article  Google Scholar 

27. Han, T. et al. Theoretical realization of an ultra-efficient thermal-energy harvesting cell made of natural materials. Energy Environ. Sci. 6, 3537–3541 (2013).

    Article  CAS  Google Scholar 

28. Chen, T., Weng, C.-N. & Tsai, Y.-L. Materials with constant anisotropic conductivity as a thermal cloak or concentrator. J. Appl. Phys. 117, 054904 (2015).

    Article  ADS  Google Scholar 

29. Li, T.-H. et al. Design of diamond-shaped transient thermal cloaks with homogeneous isotropic materials. Front. Phys. 11, 110503 (2016).

    Article  ADS  Google Scholar 

30. Wang, J., Bi, Y. & Hou, Q. Three-dimensional illusion thermal device for location camouflage. Sci. Rep. 7, 7541 (2017).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

31. Li, Q. & Vipperman, J. S. Non-singular three-dimensional arbitrarily shaped acoustic cloaks composed of homogeneous parts. J. Appl. Phys. 124, 035103 (2018).

    Article  ADS  Google Scholar 

32. Zhang, X., He, X. & Wu, L. Three-dimensional thermal illusion devices with arbitrary shape*. Chin. Phys. B 28, 064403 (2019).

    Article  ADS  Google Scholar 

33. Ji, Q. et al. Non-singular homogeneous polyhedral heat cloak and its realization. ES Energy Environ. 7, 29–39 (2020).

    Google Scholar 

34. Xiao, Y., Chen, Q. & Hao, Q. Inverse thermal design of nanoporous thin films for thermal cloaking. Mater. Today Phys. 21, 100477 (2021).

    Article  CAS  Google Scholar 

35. Chen, T. & Lin, J.-H. Exact thermal invisibility for spherical cloaks with imperfect interfaces. AIP Adv. 12, 075214 (2022).

    Article  ADS  Google Scholar 

36. Zhuang, P., Wang, J., Yang, S. & Huang, J. Nonlinear thermal responses in geometrically anisotropic metamaterials. Phys. Rev. E 106, 044203 (2022).

    Article  ADS  CAS  PubMed  Google Scholar 

37. Shan, Q., Shao, C., Wang, J. & Xia, G. Enhanced thermal invisibility effect in an isotropic thermal cloak with bulk materials. Chin. Phys. Lett. 40, 104401 (2023).

    Article  ADS  Google Scholar 

38. Xu, H., Shi, X., Gao, F., Sun, H. & Zhang, B. Ultrathin three-dimensional thermal cloak. Phys. Rev. Lett. 112, 054301 (2014).

    Article  ADS  PubMed  Google Scholar 

39. Chen, T. et al. Three-dimensional direct current invisibility cloak produced with bulk materials. Opt. Express 29, 42580–42588 (2021).

    Article  ADS  CAS  Google Scholar 

40. Zhou, Y. et al. Reconfigurable three-dimensional thermal dome. Engineering 46, 236–244 (2025).

    Article  Google Scholar 

41. Li, Y., Bai, X., Yang, T., Luo, H. & Qiu, C.-W. Structured thermal surface for radiative camouflage. Nat. Commun. 9, 273 (2018).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

42. Qin, F., Wang, L., Jiang, W., Zhang, Y. & Liu, H. Three-dimensional patterned distribution of thermal conductivity in the volume for effective thermal concentration. J. Appl. Phys. 129, 055102 (2021).

    Article  ADS  CAS  Google Scholar 

43. Xu, L. et al. Heat diffusion invariant. Phys. Rev. Lett. 135, 067103 (2025).

    Article  ADS  CAS  PubMed  Google Scholar 

44. Bendsøe, M. P. & Kikuchi, N. Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Eng. 71, 197 – 224 (1988).

    Article  MathSciNet  Google Scholar 

45. Bendsøe, M. P. & Sigmund, O. Topology Optimization: Theory, Methods and Applications (Springer, Berlin, Heidelberg, 2003).

46. Pantz, O. & Trabelsi, K. A post-treatment of the homogenization method for shape optimization. SIAM J. Control Optim. 47, 1380–1398 (2008).

    Article  MathSciNet  Google Scholar 

47. Groen, J. P., Stutz, F. C., Aage, N., Bærentzen, J. A. & Sigmund, O. De-homogenization of optimal multi-scale 3d topologies. Comput. Methods Appl. Mech. Eng. 364, 112979 (2020).

    Article  MathSciNet  Google Scholar 

48. Geoffroy-Donders, P., Allaire, G. & Pantz, O. 3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method. J. Comput. Phys. 401, 108994 (2020).

    Article  MathSciNet  Google Scholar 

49. Groen, J. P., Thomsen, C. R. & Sigmund, O. Multi-scale topology optimization for stiffness and de-homogenization using implicit geometry modeling. Struct. Multidiscip. Optim. 63, 2919–2934 (2021).

    Article  Google Scholar 

50. Stutz, F. C. et al. Synthesis of frame field-aligned multi-laminar structures. ACM Trans. Graph. https://doi.org/10.1145/3516522 (2022).

51. Lurie, K. A. & Cherkaev, A. G-closure of a set of anisotropic conducting media in the case of two dimensions. Dokl. Akad. Nauk SSSR 259, 328–331 (1981).

    MathSciNet  Google Scholar 

52. Lurie, K. A. & Cherkaev, A. V. Effective Characteristics of Composite Materials and the Optimal Design of Structural Elements, pp 175–271 (Birkhäuser Boston, Boston, MA, 1997).

53. Cherkaev, A.Variational Methods for Structural Optimization (Springer-Verlag, Berlin/Heidelberg/New York/etc., 2000).

Download references

Authors and Affiliations

1. Department of Civil and Environmental Engineering, University of Illinois Urbana-Champaign, Urbana, IL, USA

   Weichen Li, Yibo Wang & Xiaojia Shelly Zhang

2. Department of Civil and Mechanical Engineering, Technical University of Denmark, Koppels Allé, Kongens Lyngby, Denmark

   Ole Sigmund

3. Department of Mechanical Science and Engineering, University of Illinois Urbana-Champaign, Urbana, IL, USA

   Xiaojia Shelly Zhang

4. National Center for Supercomputing Applications, University of Illinois Urbana-Champaign, Urbana, IL, USA

   Xiaojia Shelly Zhang

Contributions

O.S. and X.S.Z. conceived the research and supervised the research project. W.L., Y.W., O.S., and X.S.Z. developed the methodology, carried out the theoretical and numerical investigation, and drafted and revised the paper. W.L. and Y.W. performed the experiment. X.S.Z. and O.S. secured the funding support. X.S.Z. provided the resources.

Corresponding author

Correspondence to Xiaojia Shelly Zhang.

Competing interests

The authors declare no competing interests.

Peer review information

Nature Communications thanks the anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions