Designing wideband thermal cloaks remains a challenge, especially at high frequencies. We propose an optimization approach for the design of a thermal cloak for an arbitrary object with large thermal conductivity (copper), in a given frequency band and for a specific diffusion direction. Cloaking performance is assessed as a function of frequency (in the optimal direction) and as a function of angle (at the optimal frequency). Near-perfect cloaking is achieved over a finite frequency band, and, moreover, the thermal cloak performs well in the time domain, including in the transient regime, irrespective of the initial temperature distribution. Interestingly, this specially optimized cloak also works fairly well for other objects with large thermal conductivity but breaks down for those of low thermal conductivity.

Different routes have been proposed to achieve thermal cloaking in a transient regime. The transformation thermodynamics approach leads to a highly heterogeneous and anisotropic diffusivity.^{1} In the steady state, it has been proposed that an incomplete cloak employing naturally available materials can also function very well.^{2} In addition, shielding, concentrating, and inverting heat current have been experimentally demonstrated with engineered materials in Ref. 3. However, thermal cloaking can only be achieved after certain laps of time in transient regimes.^{4} Transient thermal cloaking has been validated experimentally in the time domain with a ten layer cloak.^{5} Unfortunately, short times correspond to high frequencies for which long wave homogenization techniques require a very fine layering of the cloak; typically, ten thousand layers are required to achieve a good cloaking efficiency.^{6} One way to reduce the required number of layers in the thermal cloak would be to implement an impedance matching concept proposed for electromagnetic cloaks.^{7} However, an alternative route is to consider an active thermal cloaking approach,^{8,9} based on Peltier elements surrounding the object to cloak. A single thermal dipole suffices to achieve some cloaking in the steady state regime, as theoretically and experimentally demonstrated in Ref. 10. Here again, achieving cloaking at short transient times remains a challenge. In this article, we would like to propose a third route based on an optimization approach,^{11} which has been implemented for designing microwave cloaks over a finite bandwidth.^{12}

In the past two decades, topology optimization (TO),^{13} originally well developed in mechanics, has become a widely used tool in other areas of computational physics,^{14} such as photonics^{15} or thermal conduction, convection, and heat transfer,^{16} and has allowed the inverse design of a broad range of devices. In essence, density based TO is an inverse design procedure that can produce highly optimized structures serving a dedicated objective. One of its main advantages is that it offers unparalleled design freedom since the material distribution is updated locally (at the pixel or voxel level) inside the domain of interest. On the other hand, fabrication constraints often limit this versatility and several auxiliary tools can be included to tackle those issues,^{17} for instance, imposing minimal length-scales or ensuring the connectivity of the resulting layout. Since the number of degrees of freedom is usually prohibitively high to obtain the gradient using naive finite differences, adjoint sensitivity analysis^{18} is an indispensable part of all inverse design algorithms.

^{19,20}level set methods combined with the covariance matrix adaptation evolution strategy (CMA-ES) for the design of thermal cloaks,

^{21}bifunctional cloaks in electrostatics and static thermal diffusion,

^{22}or devices with simultaneous thermal cloaking and concentrating capabilities.

^{23}However, no study on the topology optimized thermal cloak in the frequency domain has been reported in the literature. We choose to work in the frequency regime, and thus, adapting the strategy of Ref. 11 requires turning the real wavenumber of the Helmholtz equation for propagating waves into a purely imaginary one for diffusive pseudo-waves.

^{24}Thus, we consider the time-harmonic heat equation (in distributional sense) for the temperature

*T*,

*κ*is the conductivity (W m

^{−1}K

^{−1}, i.e., watt per meter kelvin in SI units),

*ρ*is the density (kg m

^{−3}, i.e., kilogram per cubic meter in SI units),

*c*is the heat capacity (J K

^{−1}kg

^{−1}, i.e., joule per kilogram kelvin in SI units), and

*ω*is the angular frequency (rad s

^{−1}, i.e., radian per second in SI units) of a periodic heat source. A solution of (1) in a homogeneous background with conductivity

*κ*

_{b}, density

*ρ*

_{b}, and heat capacity

*c*

_{b}is given by a pseudo-plane wave

*T*

_{b}= exp(

*j*

*k*_{b}

**·**

**) with wavevector**

*r*

*k*_{b}= (

*k*

_{b}cos

*θ*,

*k*

_{b}sin

*θ*), where

*θ*is the angle of incidence and

*k*

_{b}= (1 +

*j*)/

*δ*, in which $\delta =2\kappa b/\omega \rho bcb$ is the thermal diffusion length. Our formulation then considers the diffused temperature field

*T*

_{d}=

*T*−

*T*

_{b}solution of a problem with non-homogeneous materials distribution,

*z*direction, and we choose the object to be cloaked as a cylindrical rod Ω

_{rod}centered at the origin of a Cartesian coordinate system, made of copper (

*κ*

_{Cu}= 394 W m

^{−1}K

^{−1},

*ρ*

_{Cu}

*c*

_{Cu}= 3.4 × 10

^{6}J kg

^{−1}K

^{−1}), and of radius

*R*

_{o}=

*δ*

_{0}/6, where

*δ*

_{0}= 6.8 mm corresponds to the thermal diffusion length in the background for a reference frequency

*ω*

_{0}= 1 rad s

^{−1}. It is embedded in a square computational domain Ω

_{b}of size

*l*

_{x}= 2

*δ*

_{0}and filled with iron (

*κ*

_{Fe}= 80 W m

^{−1}K

^{−1},

*ρ*

_{Fe}

*c*

_{Fe}= 3.4 × 10

^{6}J kg

^{−1}K

^{−1}), with homogeneous Neumann boundary conditions (i.e., perfect insulator with zero normal flux) applied on its outer boundaries. In the absence of cloak, the temperature field is distorted by the presence of this object as seen in Fig. 1: in (b), isothermal lines are straight parallel lines, unlike in (c). The cloak to be optimized is a shell design domain Ω

_{des}around the object of radius

*R*

_{c}= 1.5

*R*

_{o}, in which the material distribution is parameterized by a continuous scalar density function

*s*∈ [0, 1]. A filtered density $s\u0303$ is used to avoid smaller features and pathological “chessboard” patterns by solving the following Helmholtz equation

^{25}with homogeneous Neumann boundary conditions on Ω

_{des}:

*R*

_{f}being a real positive parameter characterizing the filter radius [

*R*

_{f}= (

*Rc*−

*R*

_{o})/2 in the rest of the paper]. Next, we employ a threshold projection to obtain a binary design consisting of two materials only,

^{26}

*ν*= 1/2 and

*β*is a real positive parameter and is increased during the course of the optimization. Finally, the material properties inside the cloak are defined by the solid isotropic material with penalization (SIMP) method

^{27}as

*κ*

_{Ti}= 22 W m

^{−1}K

^{−1},

*ρ*

_{Ti}

*c*

_{Ti}= 2.35 × 10

^{6}J kg

^{−1}K

^{−1}).

The gradient based optimization is initialized with a uniform density *s*_{0} = 0.5 and performed for 20 iterations or until convergence on the objective or design variables. This step is then repeated for *n* global iterations setting *β* = 2^{n}, where *n* is an integer between 0 and 7, restarting the algorithm with the optimized density obtained at the previous step as an initial guess.

*T*

_{cloak}and

*T*

_{obj}are the fields for the cloaked and uncloaked cases, respectively, and $I(T)=\u222b\Omega b|T|2dr$. This objective function essentially quantifies the reduction in the amount of energy diffused in the presence of the cloak. All the numerical results are obtained using open source libraries with bindings for the python programming language. Equations (2) and (3) are solved by the finite element method (fenics

^{28}) using second order Lagrange shape functions, the sensitivity analysis is then obtained with the adjoint method

^{29}(with automatic differentiation through dolfin-adjoint

^{30}), and the method of moving asymptotes (MMA

^{31}) is used for the optimization (nlopt package

^{32}).

The layout of the optimized cloak is a well-defined distribution of copper and titanium as seen in Fig. 1(a) and is approximately mirror symmetric along the *x*-axis (i.e., in the diffusion direction) although we did not impose any symmetry constraints. The temperature distribution is shown in Fig. 1(d), and one clearly sees the cloaking effect achieved, as we recover quasi-straight and parallel isothermal lines in the background. The final objective is *ϕ* = 2.03 × 10^{−2}; that is, the diffused thermal energy is reduced by ∼98% with respect to the object alone.

Since our target was set with a preferred direction and at a fixed frequency, we now explore the behavior of the proposed device when those parameters vary. Figure 2 shows the objective function Φ as a function of the pulsation of the incident wave. Unsurprisingly, it is minimal at the design target *ω* = 1 rad/s, and we remark that it increases slower for higher frequencies: this could be explained by the fact that at higher frequencies, the decay is stronger so that the field is weaker at the location of the cloaked object as can be seen in the insets of Fig. 2, showing the temperature distribution for *ω* = 0.1 rad/s (left) and *ω* = 6 rad/s (right). Overall, on the studied frequency band, heat diffusion is still reduced compared with the uncloaked case. This is in stark contrast to propagating waves where the invisibility effect usually relies on resonant phenomena and is, therefore, narrowband,^{11,33} although this can be mitigated by optimizing over a frequency band but with a trade-off on the performances. Furthermore, the relatively wideband behavior of the proposed design lets us envision that the cloak should perform well in the time domain, as we will investigate later in this study.

Figure 3(a) shows the diffusion reduction objective *ϕ* as a function of the angle of incidence *θ* of the pseudo-plane wave. It is minimum as expected at *θ* (corresponding to diffusion from left to right for which we performed the minimization). Interestingly, for the opposite direction (*θ* = 180°), this metric has a second local minimum at Φ = 0.072 and the diffusion is still reduced as seen in the field map in Fig. 3(b). Outside those two extremal configurations, the performances degrade and even reach a regime of superdiffusion (i.e., Φ > 1) where the energy diffused in the background is higher than for the object alone, as confirmed by the temperature map in Fig. 3(c) for *θ* = 90°, where one can observe strongly distorted isotherms. This angular dependence has been observed in electromagnetism,^{33,34} and enforcing more and more stringent symmetries on the design has been shown to increase the angular tolerance at the expense of poorer scattering reduction.

*t*= 0.01 s. In addition, one needs to choose an initial condition for

*T*.

We start by setting a homogeneous initial distribution *T*_{0} = 0 and perform the time domain simulation up to *t* = 5 s. At very short times (*t* < 0.1 s), the object and cloak have a very low diffusion as illustrated by the coefficient *K*(*t*) [see the green lines in Fig. 4(a)] since the perturbation of the temperature profile did not reach the central region where the cloak/object is placed. We can, therefore, question the meaning and accuracy of assessing the cloaking effect at short times. Cloaking is maximal (Ψ = 0.025) at *t* ≃ 0.8 s and reaches an equilibrium value of 0.05 for long times [see the green line in the inset of Fig. 4(a)]. The observation of a snapshot of the temperature field at *t* = 0.8 s confirms that the cloak performs well [see Fig. 4(d)] as illustrated by straight isotherm lines and the similarity with the free space situation [Fig. 4(b)], in contrast to the bare object [Fig. 4(c)].

*K*(

*t*). Eventually, the cloak reduces diffusion compared to the object alone and reaches a similar thermal equilibrium as in the case of homogeneous zero initial temperature, as the metric Ψ reaches the same asymptotic value for large

*t*[see the blue curve in the inset of Fig. 4(a)]. This is confirmed by the maps in Figs. 4(e)–4(g), where the field for the uncloaked case shows disturbed temperature gradients, whereas a very similar distribution is obtained for the background and cloak. Although it is impossible to confirm numerically for every initial distribution, this feature shows the robustness of the proposed design toward variations in the ambient conditions.

Finally, we investigate whether the cloak still functions for objects that are no longer made of copper. The left panel in Fig. 5 shows the diffusion reduction Ψ as a function of time for objects made of different materials (see the legend for the values of material properties). One can see that we have two cases: the performances are poor for objects with a lower conductivity such as titanium, Polyvinyl Chloride (PVC), and Polydimethylsiloxane (PDMS), and the diffusion still reduces significantly for highly conducting materials (aluminum, diamond, and graphene). Field maps on the right of Fig. 5 illustrate that insulating objects disturb the field strongly, while for conducting ones, the isotherms are straighter, demonstrating that the cloaking somehow retains its efficiency. The results also suggest that the cutoff for satisfying performances is when the conductivity of the object is higher than the conductivity of the background (iron), which puts a cutoff at *κ* = 80 W m^{−1} K^{−1}, and that the influence of *ρc* is minimal. These results suggest that our technique is versatile and would be robust to impurities or fabrication uncertainties.

In conclusion, we have proposed and numerically designed a unidirectional thermal cloak based on inverse design minimizing the diffusion. Although the optimization was carried out in the frequency domain, we have assessed its performances in the time domain and found out that it performs well, even for different initial temperature distributions. In addition, although it was designed for a given copper object, the cloak seems to perform well for objects with a high conductivity but would break down for insulating materials. These results pave the way for controlling the flow of heat with structured materials, complementing other design techniques based, for instance, on transformation thermodynamics.

Sébastien Guenneau thanks Marc Holderied and Richard Craster for insightful discussions and their support through the UK EPSRC Grant No. EP/T002654/1. Benjamin Vial and Yang Hao acknowledge the support from UK EPSRC through Grant Nos. EP/R035393/1 and EP/N010493/1.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Benjamin Vial**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Sébastien Guenneau**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Yang Hao**: Funding acquisition (equal); Project administration (equal); Supervision (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Topology Optimization: Theory, Methods, and Applications*