Design of Two-Dimensional Transient Circular Thermal Cloaks with Imperfect Interfaces

In this paper, analytic modeling for the design of a transient thermal invisibility cloak with imperfect interfaces is presented together with numerical simulations. In contrast to steady-state conditions, it is shown that an object can only be made partially invisible under a transient-state condition with either ideal or imperfect interfaces. The thermal visibility of an object to the external region can be optimally suppressed under certain conditions referred to as the “weak invisibility conditions” for the transient response, which are different from the “strong invisibility conditions” that can completely conceal an object in a steady state. In the formulation, a homogeneous metamaterial with constant volumetric heat capacity and constant anisotropic conductivity tensor is employed. It can be demonstrated that the interface’s bonding conditions will have a significant effect on the design of metamaterials. Two typical types of imperfect interfaces, referred to as low-conductivity- and high-conductivity-type interfaces, are considered. Conditions, that render an object mostly undetectable, are analytically found and expressed in simple forms under quasi-static approximations. Within the quasi-static limit, the thermal localization in the target region can be tuned with the anisotropy of the conductivity tensor. Thermal shielding or concentrating effects in the target region are exemplified based on finite element simulations to demonstrate the manipulation of heat flux in the target region. The present findings make new advances in theoretical fundamentals and numerical simulations on the effect of the imperfect interface in the transient regime and can serve as guidelines in the design of thermal metamaterials through the entire conduction process.


Introduction
Inspired by the novel innovation of an electromagnetic invisibility cloak proposed by Pendry et al. [1], thermal invisibility cloaks with tunable functionalities have attracted substantial attention in recent years [2][3][4][5][6][7][8][9][10][11][12][13] for various applications in advanced thermal management. Thermal metamaterials are often composed of different materials with internal structures. The bonding conditions between different constituent materials will in general affect the heat-transferring behavior across the interface, causing diminishing or even breaking down the functionalities of the metamaterials. For example, Zheng and Li [14] demonstrated that, in the presence of interfacial thermal resistance, the temperature field outside the cloak could be greatly distorted, if the interfaces are inappropriately considered to be perfect. In practical situations, imperfectness always exists across the interface between two different solids in contact owing to the interfacial thermal resistance [15] or roughness at their common boundaries, or caused by the mismatch of physical properties between the adjacent phases [14]. The effect of imperfect interfaces will become particularly significant and non-negligible while the thermal device is in the nano-scaled [15,16]. However, most previous studies mainly considered ideal interfaces including steady-state thermal cloak made of inhomogeneous conductivity and constant volumetric heat capacity, which works perfectly in the steady state, is always detectable with a small amount of scattered heat in the transient state. Ji et al. [26] explored the cloaking of a complex shape on the basis of either the concept of neutral inclusion or the transformation method in both steady-state and transient regimes and showed that the neutral inclusion method is more flexible and easier in practical implementation. Zhang et al. [33] reported that the transient heat flow velocity can be controlled with the regulation of effective density and specific heat capacity via a proposed thermal architected metamaterial. To our knowledge, the subject of the present context is new in the literature.

Results
Based on the first principle of thermodynamics, the temperature field T is governed by Fourier's relation ∇·k∇T = ρc p (∂T/∂t) in the absence of a heat source, in which k is the conductivity tensor, ρ the mass density, and c p the specific heat capacity at constant pressure. The total response of the temperature field consists of the transient and the steadystate responses. The former will decay with time while the latter governs the long-term response. Apart from the steady-state response, the volumetric heat capacity s = ρc p plays an essential role in a transient response during the heat conduction process. When the volumetric heat capacity s or ∂T/∂t is very small, the state change of the system will be negligibly small and thus the heat conduction process will asymptotically proceed under the "quasi-static process". Incorporating the effect of imperfect interfaces, here we will extend the framework from our previous work [27] to explore how a thermal cloak can be properly designed in the transient regime under the quasi-static approach.
A two-dimensional model is illustrated in Figure 1. Region I is the target region for heat flux control, made of the same material as the surrounding host medium Region III with constant isotropic conductivity k 0 and volumetric heat capacity s 0 . Region II is the designed metamaterial with constant volumetric heat capacity s 2 , and constant anisotropic conductivity tensor, expressed as where k r and k θ are respectively the constant radial and circumferential conductivities. Both of the interfaces at r = a and r = b are considered as either LC-or HC-type imperfect interfaces. It is mentioned that, by eliminating the effect of imperfect interfaces, the analytic result for the case of perfect interfaces will be readily obtained. The transient thermal excitation T inc = U inc e −iωt = e i(K 0 x−ωt) is applied horizontally along x axis from the lefthand side of the thermal device, where K 0 = √ iωs 0 /k 0 refers to the wave number of the pseudo diffusion plane wave in Region III. The scalar temperature field U inc can be expanded in cylindrical coordinates by the Jacobi-Anger expansion [34], given as U inc (r, θ) = e iK 0 x = e iK 0 rcos θ = ∞ ∑ n=−∞ i n J n (K 0 r)e inθ , (2) in which J n refers to the Bessel function of the first kind. Taking the boundary conditions and geometric symmetry into account, the temperature field for the transient response can be assumed as T(r, θ, t) = U(r, θ)e −iωt . As such, the scalar temperature fields that comply with Fourier's relation can be expressed in cylindrical coordinates as [34] U 1 (r, θ) = ∞ ∑ n=−∞ A n J n (K 0 r)e inθ , r ≤ a, where H (1) n is Hankel function of the first kind with order n, K 2 = √ iωs 2 /k r is the wave number in Region II, and the n th -order anisotropy parameter is equal to λ = √ k θ /k r multiplied by n, given as λ n = n √ k θ /k r = nλ. The value of λ represents the anisotropy of the conductivities in Region II. A value of λ = 1 stands for an isotropic material in Region II. The coefficients A n , B n , C n and S n will be determined by the interface conditions at the interfaces. To make the thermal cloak undetectable, the scattering coefficients S n have to be canceled out. However, for the transient response, it is not possible of eliminating the infinite scattering modes S n by a limited number of boundary conditions and interface conditions. As such, one can only get rid of the primary modes of scattering fields to suppress the entire scattering fields to the maximum extent such that the thermal cloak can be made mostly invisible. The scattering width σ 2D , which is introduced to measure the overall visibility of the inclusion to external observers, can be defined as [35] where T sc = T 3 − T inc denotes the scattering fields in Region III. The subscript 2D refers to the two-dimensional or cylindrical cases. With the substitution of (3), it will be shown that interfaces in Sections 2.1-2.3. The main results will be given in Equations (7) and (8) for ideal interfaces, Equations (10) and (11) for LC-interfaces, and Equations (15) and (16) for HC-interfaces.

Ideal Interface
The case of ideal interfaces will be examined first. In a steady state, the design of a thermal invisibility cloak only involves the connection that balances the equivalence in effective conductivities between the surrounding background material and the cloaked inclusion. Additionally, in the transient regime, the volumetric heat capacity will take effect on the dynamic balance as well. In the quasi-static limit, it is analytically found that when the universal connections = , for = 0 (monopole mode), are fulfilled, the cloak can be made least visible to the exterior observers. It is of interest that these two constraint conditions confine separately the balance in volumetric heat ca- is applied on the left-hand side, where ω is the angular frequency, and K 0 = √ iωs 0 /k 0 the wave number. Region I and Region III are made of the same isotropic material with constant conductivity k 0 and constant volumetric heat capacity s 0 . Region II is the designed anisotropic metamaterial with constant volumetric heat capacity s 2 , constant conductivities k r and k θ .
The detailed derivation of (5) is given in Appendix A. It is remarkable that in the quasi-static approximation, ∂T/∂t = −iωT shall be very small. This requires a very small frequency ω → 0 as well as a long diffusion length in each region, namely K 0 r 1 in Region I and K 2 r 1 in Region II. In this scenario, the scattering fields are mostly contributed by the primary scattering orders, i.e., the first two orders (n = 0 for the monopole and n = 1 for the dipole mode). Consequently, the scattering width will be asymptotically equal to It turns out that, on the basis of the concept of SCT [31], the scattering fields can be mostly suppressed when S 0 and S 1 are zero. This will give the "weak invisibility conditions" for the transient thermal cloak that can mostly conceal the inclusion. It is remarkable to note that this also coincides with the results for the spherical case that the perturbation around the thermal cloak is primarily governed by the first two modes of the scattered field as well in the quasi-static limit, as illustrated in [24]. We mention that the thermal cloak that fulfills weak invisibility physically implies that the effective material parameters of the thermal device are equivalent to those of the background medium in an averaging sense [32]. This concept also corresponds to that proposed by [36] for determining the effective shear modulus of a composite material based on a generalized self-consistent method. In the following, these exact connections for different types of bonding conditions will be derived in a simple form. In the following the universal connections as the criteria for the design of homogeneous metamaterial to optimally achieve thermal invisibility will be analytically examined respectively for the cases of ideal interfaces, LC-and HC-type interfaces in Sections 2.1-2.3. The main results will be given in Equations (7) and (8) for ideal interfaces, Equations (10) and (11) for LC-interfaces, and Equations (15) and (16) for HC-interfaces.

Ideal Interface
The case of ideal interfaces will be examined first. In a steady state, the design of a thermal invisibility cloak only involves the connection that balances the equivalence in effective conductivities between the surrounding background material and the cloaked inclusion. Additionally, in the transient regime, the volumetric heat capacity will take effect on the dynamic balance as well. In the quasi-static limit, it is analytically found that when the universal connections are fulfilled, the cloak can be made least visible to the exterior observers. It is of interest that these two constraint conditions confine separately the balance in volumetric heat capacities and conductivities. In the quasi-static limit, we find that the volumetric heat capacity is necessary to be homogeneous throughout based on the monopole constraint (7), whereas the dipole mode (8) gives an identical connection to the invisibility condition for the steadystate response [37]. This implies that the design of thermal invisibility for the transient response will work perfectly in a steady state. As will be mentioned later, depending on whether the interfaces are of LC-or HC-type, the effective conductivities of the cloaked region will be diminished or enhanced. To signify the effect of the imperfectness at the interfaces, let us define an imperfectness parameter g = k 0 /k G = k 0 /(λk r ). Note that k G = λk r = √ k r k θ refers to the geometric mean of the conductivities of the metamaterial in Region II. Thus (8) can be simplified as g = 1. This also indicates that, in the case of ideal interfaces, the connections for thermal invisibility only depend on material properties s 2 and k G for the metamaterial, regardless of geometric size or the anisotropy ratio. When the imperfectness of the interfaces starts to increase, the value of g will begin to deviate from unity to achieve thermal invisibility. A value of g < 1 and g > 1 will be invoked for the cases of LC-and HC-type interfaces to compensate for the reduction or enhancement of effective conductivities respectively. The size effect can be observed when g = 1.

Low Conductivity-Type Interface
The LC-type interface exhibits thermal contact resistance impeding heat flowing across the interfaces, also known as Kapitza resistance [15]. It is mentioned that the interface conditions remain the same for steady and transient states for LC-type interfaces [38].
Specifically, the interface conditions for LC-type interfaces at r = a and r = b can be expressed as [39,40] here T 1 , T 2 , and T 3 are temperature fields in Region I, II, and III, respectively, and ∂/∂n denotes the normal derivative at interfaces. The interface parameter β is defined as the ratio between the temperature jump and the normal heat flux across the interface, with a dimension of W·m −2 K −1 . The subscripts a and b in β refer to the interfaces at r = a and r = b respectively. A value of β → ∞ represents an ideal interface, while in contrast, the case that β → 0 stands for the adiabatic contact. With the substitution of (3) into (9), the weak invisibility conditions for the LC-type can be obtained by eliminating S 0 and S 1 simultaneously, given as where c = (a/b) 2 is the area fraction of Region I. The non-dimensional interface parameters are defined the same as those in Reference [27] for consistency. The detailed derivation is given in Appendix B. Regarding the effect of imperfectness at the interfaces, it is remarkable to note that, when β a → ∞ and β b → ∞ , that is ∼ β a = ∼ β b = 0, the case of LC-type interfaces will reduce back to that of ideal interfaces. When the values of ∼ β a and ∼ β b increase as the thermal resistance grows at the interfaces, the effective conductivities of the metamaterial have to increase correspondingly to compensate for the effect of impedance across the interfaces. As such, the value of k G = √ k r k θ shall be larger than the conductivity in the background medium k 0 , which leads to g = k 0 /k G < 1 for the LC-type cases. Unlike the perfect bonding cases, the weak thermal invisibility relates not only to the material properties s 2 and k G , but also the imperfectness parameter g, area fraction c = (a/b) 2 , and λ = √ k θ /k r . In addition, the LC-type parameters ∼ β a = k 0 /(aβ a ) and ∼ β b = k 0 /(bβ b ) exhibit size effect for the design of a thermal cloak. The effect of thermal resistance on the thermal cloak will become more significant with scaling down the thermal device. For a detailed exposition of the effect of interfaces, one can also refer to our previous work [27]. In addition, under the quasi-static approximation, we can find these two universal conditions are decoupled, namely, the monopole connection (10) constrains the volumetric heat capacity, while the dipole mode (11) rules the conductivities in Region II separately. We mention that the monopole and dipole constraints will be coupled, related to both volumetric heat capacities and conductivities simultaneously if the process of heat conduction does not proceed slowly enough to fulfill the quasi-static assumption. Only when the rate of temperature change is relatively small during the conduction process will the constraint conditions for these two modes be decoupled. In addition, it can be seen that (11) is exactly the same with the strong invisibility condition for a steady-state response for the case of LC-type interfaces as shown in Reference [27]. These results indicate that, in the quasi-static limit, the system is nearly in static equilibrium at each instant of time during the heat diffusion process. In this scenario, the dipole mode solely connects the steady-state behavior while the time effect of the transient response reflects entirely on the monopole constraint. When the volumetric heat capacity is negligibly small, the dynamic response will asymptotically approach the steady-state response. In this case, one can observe that the monopole constraint (10) becomes trivial, whereas the dipole mode alone governs the universal connection for invisibility entirely. Physically this implies that the inclusion, which can be least revealed with the design via the connections (10) and (11) in a transient state, will be completely concealed as a steady state is established after a certain period of time. This phenomenon also agrees well with the numerical results to be discussed later in Section 5.
Additionally, it is worthwhile to mention that, if the volumetric heat capacity in Region I, denoted as s 1 , is considered different from that in Region III, the monopole condition will become The above equation can be rewritten as This result is analogous to that found in Reference [24] for a three-dimensional problem. When s 1 = s 0 , (12) will give c = 1 or s 2 = s 0 . The former is a trivial result indicating that the whole medium is homogeneous throughout. To our concern, the latter s 2 = s 0 is selected as the constraint condition for monopole mode, as shown in (10).

High Conductivity-Type Interface
The HC-type interface, in contrast to the LC-type interface, can be modeled in terms of a highly conducting interphase layer [30,41]. The interface conditions of an HC-type interface are derived based on the dynamic balance within the interphase layer with an asymptotic approach. The volumetric heat capacity of the interphase s int have to be taken into account to balance the discrepancy resulting from the discontinuous drop of heat flux in a transient state. A detailed exposition of the derivation of the HC-type interface condition with transient effect is given in Appendix C. The HC-type interface conditions at the interfaces can be expressed as [27] ( where the surface volumetric heat capacity of the interphase is defined as s s = s int h, and h the thickness of the interphase. Nevertheless, s s is negligible as s int is finite and h → 0 . The HC-type interface conditions in the transient regime can be reduced back to that in a steady state [29]. The parameters of the HC-type interfaces are represented by α a and α b respectively for the interface at r = a and at r = b, with a dimension of W·K −1 . In the quasi-static limit, we find that, when the material parameters and geometric parameters fulfill the universal connections the visibility to the exterior observers can be minimized. As the transient effect diminishes and only the steady-state response remains, the condition for thermal invisibility will solely rely on the dipole mode constraint (16), which is indeed the strong invisibility condition as introduced in [27]. Again, it is noted that both monopole and dipole conditions will depend on volumetric heat capacities and conductivities simultaneously if the quasi-static approximation, namely K 0 r 1 in Region I and K 2 r 1 in Region II, is not applicable. In (16), the dimensionless parametersα a = α a /(ak 0 ),α b = α b /(bk 0 ) are defined to indicate the bonding situations for the HC-type interfaces. It will give ideal interfaces when α a → 0 , α b → 0 , namelyα a → 0 andα b → 0 . The connections for weak thermal invisibility will Materials 2023, 16, 2297 8 of 20 reduce to that of ideal interfaces, which leads to g = 1 for the dipole mode condition. The value ofα a,b will increase accordingly with the increase of g, which implies that k G have to be lowered down to balance the high conducting interfaces to achieve thermal invisibility. The case that α → ∞ represents the limiting situation that the interfaces are infinitely superconducting [27,42]. Together with the LC-and HC-type interfaces, the relation between the interface parameters ∼ β a andα a versus g can be illustrated in Figure 2 in [27]. Similar to the case of LC-type interfaces, the thermal invisibility for the case of HC-type interfaces depends on volumetric heat capacity s II , as well as g, c, and λ. One can observe thatα a = α a /(ak 0 ) andα b = α b /(bk 0 ) will become noticeable in magnitude while scaling down the value of a and b. As a result, the effect of HC-type interfaces will be non-negligible with small-scale thermal devices. static approximation, namely ≪ 1 in Region I and ≪ 1 in Region II, is not applicable. In (16), the dimensionless parameters = /( ), = /( ) are defined to indicate the bonding situations for the HC-type interfaces. It will give ideal interfaces when → 0, → 0, namely → 0 and → 0. The connections for weak thermal invisibility will reduce to that of ideal interfaces, which leads to = 1 for the dipole mode condition. The value of , will increase accordingly with the increase of , which implies that have to be lowered down to balance the high conducting interfaces to achieve thermal invisibility. The case that → ∞ represents the limiting situation that the interfaces are infinitely superconducting [27,42]. Together with the LC-and HC-type interfaces, the relation between the interface parameters and versus can be illustrated in Figure 2 in [27]. Similar to the case of LC-type interfaces, the thermal invisibility for the case of HC-type interfaces depends on volumetric heat capacity Ⅱ , as well as , , and . One can observe that = /( ) and = /( ) will become noticeable in magnitude while scaling down the value of and . As a result, the effect of HC-type interfaces will be non-negligible with small-scale thermal devices.  Figure 2. Temperature contours, together with temperature and heat flux profiles based on finite element simulations (COMSOL) for isotropic Region II with LC-and HC-type interfaces. The imperfectness parameter is specified as = 2/3 in (a) for LC-type, and = 3/2 in (b) for HC-type interfaces. The temporal frequency is = 10 . For numerical simulations, we consider = 1 μm and = 1/4. We plot the contour of the temperature difference between the concerned cases and the reference case illustrated in Figure A2. The temperature contours are shown to demonstrate the temperature field for the transient responses for various time periods in (i). The temperature difference between our concerned cases and Figure A2 are plotted in (ii) to highlight the scattering field in Region III. The temperature profile as well as the heat flux profile along the -axis are shown in (iii) and (iv). The dashed lines in the profile figures are the results in Figure A2. Note that the value of (in W⋅K −1 ⋅m −1 ) is listed in each case for reference. The temperature is described in Celsius scale (°C), while the heat flux is in MW/m 2 .

Theoretical Analysis
The design of a thermal cloak in the transient regime, distinct from that in the steady regime, relates not only to the geometric parameters and conductivities but also to the volumetric heat capacity in each region. However, based on that we have assumed a constant volumetric heat capacity in both Region I and Region III, the volumetric heat The imperfectness parameter is specified as g = 2/3 in (a) for LC-type, and g = 3/2 in (b) for HC-type interfaces. The temporal frequency is ω = 10 −4 . For numerical simulations, we consider b = 1 µm and c = 1/4. We plot the contour of the temperature difference between the concerned cases and the reference case illustrated in Figure A2. The temperature contours are shown to demonstrate the temperature field for the transient responses for various time periods in (i). The temperature difference between our concerned cases and Figure A2 are plotted in (ii) to highlight the scattering field in Region III. The temperature profile as well as the heat flux profile along the x-axis are shown in (iii) and (iv). The dashed lines in the profile figures are the results in Figure A2. Note that the value of k C (in W·K −1 ·m −1 ) is listed in each case for reference. The temperature is described in Celsius scale ( • C), while the heat flux is in MW/m 2 .

Theoretical Analysis
The design of a thermal cloak in the transient regime, distinct from that in the steady regime, relates not only to the geometric parameters and conductivities but also to the volumetric heat capacity in each region. However, based on that we have assumed a constant volumetric heat capacity s 0 in both Region I and Region III, the volumetric heat capacity in Region II have to be identical to the neighboring regions such that the monopole constraint for thermal invisibility can be fulfilled in the quasi-static limit, irrespective of the bonding conditions. The dipole mode invisibility condition will prevail while the temporal frequency ω ≈ 0, thereby, governing the long-term response as time elapses. Different from those for steady-state regimes, the thermal localization in Region I is difficult to be precisely manipulated at every instant of time for a transient response. However, within the quasi-static limit, we find that the manipulation of thermal localization in Region I is the same as that in the steady-state regime [27], still mainly controlled by the anisotropy ratio λ = √ k θ /k r . The perturbation of the controlled heat flux due to the transient response in the target region will eventually goes out with time. The value of volumetric heat capacity will affect the amount of perturbation over the transient response, and govern how long it will take for a steady state established.
When the value of k θ , or more precisely the value of λ, is greater than a critical value λ cr , the heat flow will tend to be guided around Region II circumferentially. This will impede heat flux from entering Region I, exhibiting a shielding effect. Concentrating effect is opposite to the shielding behavior. When the radial conductivity k r in Region II is increased to a level such that λ < λ cr , the heat flux will be directed toward the core region, leading to a concentrating effect in Region I. A value of λ = λ cr indicates the case without thermal localization in Region I. As mentioned in [27], the value of λ cr is highly pertinent to the bonding conditions at the interfaces. In the case of ideal interfaces, one has λ cr = 1, indicating that the thermal shielding (λ > 1) or concentrating effect (λ < 1) can be simply achieved by the value of λ. When there is imperfect bonding at the interfaces, the value of λ cr will be reduced in magnitude. This gives λ cr < 1 for either LC-or HC-type interfaces, as illustrated in [27]. This phenomenon physically explains by the fact that there shall be some of the heat energy dissipated while heat flow passes across the imperfect interface. As such, a shielding effect will be much more achievable than a concentrating effect in the presence of imperfect interfaces. In addition, one can find that the amplification ratio of concentrating effect is always bounded by the geometric parameter b/a as the steady state is nearly established, while a perfect shielding effect can be always achievable regardless of the bonding conditions at interfaces.
In Table 1 we list the weak invisibility conditions for a transient-state response in comparison with the strong invisibility conditions for a steady-state response [27]. Obviously, it can be seen that the weak invisibility condition in regard to the dipole mode for a transient response is identical to the strong invisibility condition for the steady-state case. It is of interest to mention that, mathematically the eigenfunctions of the dipole mode (n = 1) in (3) will asymptotically approach the polynomials composed of r ±λ , which correspond to the eigenfunctions in the steady-state regime. This also offers an interpretation, explaining why the dipole mode constraint in the transient regime will be identical to the thermal invisibility condition in a steady state.

Methods
In addition to a homogeneous volumetric heat capacity, namely s 2 = s 0 , the design of a thermal cloak in the transient regime under the quasi-static limit can be constructed via a strategy as that proposed for the steady-state response [27]. According to the design methods, the design of a thermal cloak can be constructed via the concept of neutral inclusions [30,43] through a series of homogenization procedures based on a composite cylindrical assemblage (CCA) [44] model. Incorporating the effect of an imperfect interface, this method provides a systematic way to decompose the derivations into three separate steps, greatly simplifying the algebraic complexity. The first step will involve the homogenization of Region I coated with an imperfect interface at r = a. Secondly, the homogenization for the composite cylinder of the coated core region enclosed by the designed metamaterial, Region II, is determined. To optimally achieve invisibility, the effective material properties of the entire inclusion that concurrently concern the outer imperfect effect at r = b will be equivalent to that of the background medium under weak invisibility conditions. It is notable that only the second step is required for the cases of ideal interfaces. In Section 5, the numerical simulations based on the finite element method will be illustrated to validate our theoretical methodology.
This idea can be extended to design a multi-layered thermal metamaterial via a series of homogenization by applying two fundamental techniques recursively [27]. This can be utilized to realize a suitably designed thermal metamaterial with tailored orthotropic material properties and offers an excellent alternative in comparison with the classical method [45]. For the realization of a thermal metamaterial, the classical method suggests that a metamaterial with cylindrically anisotropic conductivities k r = 2 κ −1 and k θ = (k A + k B )/2 can be fabricated as a layered heterogeneous composite material consisting of two alternating materials with isotropic conductivities k A and k B [45]. Nevertheless, this asymptotic approach requires a sufficient number of layers to render the cloak asymptotically equivalent to that made of a homogeneous material with cylindrically orthotropic conductivities k r and k θ , especially when the values of k A and k B differ greatly. In contrast, it is of our interest that, for the design of thermal cloaks via the concept of neutral inclusions, the thermal invisibility conditions not only can be always fulfilled irrespective of the number of annular layers but also be capable of counting into the effect of bonding imperfections at each interface.

Numerical Simulations
The thermal cloak based on the universal connections (10) and (11) for LC-type interfaces, and (15) and (16) for HC-type interfaces will be demonstrated to exhibit optimal cloaks that can suppress the disturbance in Region III to the maximum extent. In Figure A2 the reference case in which the whole medium is a homogeneous material is illustrated. To meet the requirement for the quasi-static assumption, a relatively small angular frequency of the incident excitation is taken as ω = 10 −4 . Additionally, it is mentioned that the requirement of the quasi-static approximation, which is confined by K 0 b = √ iωs 0 /k 0 b 1, can be guaranteed with a small size of the radius b, low value of s 0 , or large k 0 as well. When the size of the thermal device is much smaller than the diffusion length, the scattered field will be very small and decay with time quickly. Thus a relatively small thermal device with a radius a = 0.5 µm and area fraction c = 1/4 is considered. On the other hand, a value of small volumetric heat capacity physically will barely store energy with time and converge to a steady state in a very short time. For numerical demonstration, a small value of s 0 = 1 MJ·K −1 ·m −3 is specified with appropriate thermal conductivity of k 0 = 1 W·K −1 ·m −1 . Note that the volumetric heat capacities of natural solid materials vary roughly from about unity to 4 MJ·K −1 ·m −3 [46], for example, lead at 1.44 MJ·K −1 ·m −3 , and iron at 3.569 MJ·K −1 ·m −3 . The value of thermal conductivity of natural solids varies a great deal ranging from 0.02 W·K −1 ·m −1 (aerogel, recognized currently as the solid material with the lowest thermal conductivity) to about 2000 W·K −1 ·m −1 (for example diamond).
In Figure 2, the case of λ = 10 0 is considered which refers to an isotropic Region II. For numerical simulations, the imperfect interface can be simulated as an ultrathin interphase layer with thickness h and isotropic conductivity k C since the interfaces between two connected materials are always considered perfectly bonded in the finite element analysis software COMSOL. The value of k C is related to the imperfect interface parameters β and α in such a way that k C = βh [39] and k C = α/h [29] for the LC-and HC-type interfaces respectively. As a result, we consider two interphase layers ranging from a to a + h and b to b + h respectively. For the cases of LC-type interfaces, the value of k C shall be extremely small compared with the conductivities of neighboring phase materials; while in contrast, k C will reach a value much higher than that of the adjacent contact materials for HC-type interfaces. Theoretically k C will approach zero or goes to unbounded respectively for the cases of LC-or HC-type interfaces while the thickness h → 0 . Here a sufficiently small value of h = 10 −3 b = 10 −3 µm is employed to give a convergent solution. The corresponding value of k C for each case is listed in Figure 2. The parameter that signifies the bonding situations is specified as g = 2/3 in Figure 2a and g = 3/2 in (b) for the cases of LC-and HC-type interfaces respectively. The transient excitation T inc = U inc e −iωt is applied on the left-hand side, and −T inc on the right-hand side to possess the symmetric condition. In the first row in Figure 2, the temperature contours are plotted to demonstrate the temperature field for the transient responses at various periods of time. According to the concept of SCT, the perturbation in Region III shall be minimized while the weak invisibility conditions are fulfilled. The temperature difference between our concerned cases and Figure A2 below the temperature contours, denoted as (ii), are plotted to highlight the scattering field in Region III and demonstrate the shielding effect or concentrating effect in the core region. The temperature and heat flux profiles along x-axis are illustrated respectively in (iii) and (iv) below the contours to highlight the thermal localization as well as the jump across the interfaces on temperature or heat flux for LC-or HC-type cases respectively. The dashed lines in the profile diagrams represent the results in Figure A2. When the profiles of the studied case (solid lines) overlap with the referenced case (dashed lines) in Region III, the anisotropic region will be invisible to the background medium. From the results presented in [27], one can observe that an isotropic metamaterial in Region II with imperfect interfaces shall exhibit a shielding effect since λ cr < 1 for either type of imperfect interface.
It can be seen from the profile diagram in Figure 2a,b that the heat flux in Region I is reduced in amount compared with the dashed lines, exhibiting a shielding effect as it should be. On the other hand, the dashed lines in the temperature and flux profiles are nearly overlapped with the solid lines in Region III, resulting in excellent invisibility over time. As time elapses, the temperature difference in Region III will decrease. A steady-state response has been asymptotically established at the time t = 6 s, at which the thermal cloak is almost invisible, and thus the temperature, as well as heat flux in Region I and III, will be linearly distributed along the horizontal direction, as can be seen in the profile diagrams.
To demonstrate the thermal localization in Region I manipulated with the anisotropy ratio in Region II, the cases of λ = 10 −1 and λ = 10 1 are presented in Figure 3. In Figure 3a,c, the radial conductivity is one hundred times the circumferential conductivity, namely k r = 10 2 k θ . As can be seen in the temperature contour diagram, the heat flux is directed toward the core region, and an evident concentrating effect is observed. In contrast, in Figure 3b,d we have k θ = 100k r for the case of λ = 10. The heat flux diffusing from Region III is guided around Region II, exhibiting an excellent shielding effect in Region I.

Conclusions
This work presents a detailed design of a thermal cloak with imperfect bonding interfaces focusing particularly on the transient state. It is shown analytically and demonstrated numerically that an object can only be made partially concealed under a transientstate condition with or without the effect of imperfect interfaces. The universal connections, also referred to as "weak invisibility conditions", that can optimally suppress the visibility for the design of a thermal cloak made of homogeneous anisotropic material and constant volumetric heat capacity, were analytically found. It is remarkable to see that the For numerical illustration, the imperfectness parameter g and anisotropy ratios λ are specified as (a) g = 2/3 and λ = 10 −1 ; (b) g = 2/3 and λ = 10 1 ; (c) g = 3/2 and λ = 10 −1 ; and (d) g = 3/2 with λ = 10 1 . Note that the cases of g = 2/3 will invoke the LC-type interfaces; whereas g = 3/2 refers to the cases of HC-type interfaces. The geometric size and applied boundary conditions are identical to those described in Figure 2.

Conclusions
This work presents a detailed design of a thermal cloak with imperfect bonding interfaces focusing particularly on the transient state. It is shown analytically and demonstrated numerically that an object can only be made partially concealed under a transient-state condition with or without the effect of imperfect interfaces. The universal connections, also referred to as "weak invisibility conditions", that can optimally suppress the visibility for the design of a thermal cloak made of homogeneous anisotropic material and constant volumetric heat capacity, were analytically found. It is remarkable to see that the weak invisibility conditions comprise two constraint conditions that correspond to two distinct modes, monopole and dipole modes. In the quasi-static limit, the monopole constraint condition, which confines the volumetric heat capacity of the metamaterial, characterizes the time effect on the invisibility condition. On the other hand, the dipole mode constraint condition, which only relates to the connections of conductivities, corresponds to the strong invisibility condition in the steady-state regime. This finding suggests that a thermal cloak designed via weak invisibility conditions will work excellently as an optimal solution that can minimize the visibility to the observers in the transient regime, and will turn into an ideal cloak whenever the steady state is achieved. Additionally, we also show that similar to those in the steady-state regime, the thermal localization in the quasi-static limit can be mainly controlled by the anisotropy of conductivities, and will be greatly influenced by the bonding conditions at the interfaces. A shielding effect is found to be always achievable. While the effect of concentrating will tend to slim down, especially when the imperfection arises at the interfaces. It is mentioned that the design concept could be applied to more complex geometric shapes, such as ellipses or ellipsoids. However, the mathematics involved will be much more complicated.
The present findings make new advances in theoretical fundamentals and numerical simulations on the effect of the imperfect interface in the transient regime and can serve as a guideline in the design of thermal metamaterials through the entire conduction process.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Scattering Width for a Two-Dimensional System
In this Appendix A, we provide a detailed derivation for Equation (5). The scattering width σ 2D , defined in (4), is utilized to measure the visibility of the scattered field away from the inclusion [35]. Under the quasi-static approximation, with the substitution of T(r, θ, t) = U(r, θ)e −iωt and (3) 3 into (4), the scattering width can be expanded as As dΩ = rdθ, the integration with respect to θ simply yields the value of 2π. As such, the value of σ 2D can be simplified as With consideration of the asymptotic expansion for Hankel function of the first kind with large argument [34] H (1) the scattering width in (A2) can be expressed as in which K 0 = ωs 0 /(2k 0 ) is the real part of the wave number K 0 , such that K 0 = (1 + i)K 0 as introduced in our main text. In the quasi-static limit, the value of K 0 is relatively small, thereby e −K 0 r → 1 . Consequently, we have Note that |K 0 | = √ 2K 0 in (A5). Since the values of S n are independent of the radial coordinate r, the formula for determining the scattering width can be reduced to (5).

Appendix B. Derivation of Weak Thermal Invisibility Conditions
In this Appendix B, we provide a derivation for the exact connections of thermal invisibility in the transient regime with imperfect interface effects. Specifically, a detailed derivation for the universal connections for thermal invisibility (10), (11), (15), and (16) in the transient regime under the quasi-static limit assumption will be given in detail. In cylindrical coordinates, the heat conduction equation ∇·k∇T = s(∂T/∂t) can be written as With the assumption of T(r, θ, t) = U(r, θ)e −iωt , (A6) can be rewritten as where K = √ iωs/k r refers to the wave number, k θ and k r denote the circumferential conductivity and radial conductivity respectively. Note that for Region I and Region III, k θ = k r = k 0 . In Region II the anisotropy parameter is denoted as λ = √ k θ /k r . With the variable decomposition of U(r, θ) = U(r)e inθ , (A7) will become a Bessel differential equation whose solutions consist of two series of independent Bessel functions. As such, the scalar temperature fields that fulfill Fourier's relation can be expressed in cylindrical coordinates as described in (3).

Appendix B.1. LC-Type Interfaces
For LC-type interfaces, the interface conditions (9) can be rewritten in matrix form with the substitution of (3) as where x = A n B n C n S n T is the undetermined coefficient vector. The matrix ∼ A and vector ∼ b are given as Note that the prime notation stands for the derivative with respect to the argument Kr. To our concern, the scattering coefficient S n can be derived from Cramer's rule as S n = ∆ S n /∆, in which ∆ = det ∼ A . The value of ∆ S n is given by λ n (K 2 a) 0 −k 0 K 0 J n (K 0 a) k r K 2 J λ n (K 2 a) k r K 2 H (1) The boxed terms in (A9)-(A11) are marked to highlight the difference between the case of LC-type interfaces and that of perfect interfaces. It will recover back to the case of perfect interfaces when all the boxed terms vanish. Based on the concept of the scattering cancellation technique [31], we are aiming to eliminate the first two orders of scattering fields, namely S 0 and S 1 . Equivalently, the weak invisibility conditions can be derived by letting ∆ S 0 = 0 and ∆ S 1 = 0 simultaneously. This will give the results shown in (10) and (11) in the main text with the quasi-static approximation.

Appendix B.2. HC-Type Interfaces
With the substitution of the admissible temperature fields (3), the interface conditions for the HC-type cases (14) can be rewritten in matrix form Ax =b.
The matrixÂ and vectorb are given as Figure A1. A schematic illustration of a small element on the curved interface subjected to the heat flux on the top face which is in contact with the matrix (denoted by the superscript ), on the bottom face which is in contact with the inclusion (denoted by the superscript ), and along the four edges (with superscript ).

Appendix D. Numerical Simulations of a Homogeneous Background Medium
In Figure A2, we illustrate the numerical results for the reference case in which the whole medium is homogeneously made of a material with isotropic conductivity k 0 and constant volumetric heat capacity s 0 . Hence in this case Region III 2 will remain thoroughly undisturbed during the process of conduction. The geometric configuration is the same as that in Figure 2. Temperature and flux profiles are illustrated below the contours.

Appendix D. Numerical Simulations of a Homogeneous Background Medium
In Figure A2, we illustrate the numerical results for the reference case in which the whole medium is homogeneously made of a material with isotropic conductivity and constant volumetric heat capacity . Hence in this case Region III 2 will remain thoroughly undisturbed during the process of conduction. The geometric configuration is the same as that in Figure 2. Temperature and flux profiles are illustrated below the contours.  Figure A2 is a homogeneous medium throughout which is taken as the reference medium in comparison with the designed cases. For numerical simulations, we consider = 1 μm and = 1/4. The temporal frequency is = 10 . The temperature profile as well as the heat flux profile along the -axis are shown at the bottom. As a reference case, all the profile diagrams are illustrated with dashed lines to be compared with others. The temperature is described as Celsius (℃), while the heat flux is in MW/m 2 . Figure A2. Temperature contours, together with heat flux profiles based on finite element simulations (COMSOL) for a background medium. Note that the configuration of Figure A2 is a homogeneous medium throughout which is taken as the reference medium in comparison with the designed cases. For numerical simulations, we consider b = 1 µm and c = 1/4. The temporal frequency is ω = 10 −4 . The temperature profile as well as the heat flux profile along the x-axis are shown at the bottom. As a reference case, all the profile diagrams are illustrated with dashed lines to be compared with others. The temperature is described as Celsius (°C), while the heat flux is in MW/m 2 .