Directed Thermal Diffusions through Metamaterial Source Illusion with Homogeneous Natural Media

Owing to the utilization of transformation optics, many significant research and development achievements have expanded the applications of illusion devices into thermal fields. However, most of the current studies on relevant thermal illusions used to reshape the thermal fields are dependent of certain pre-designed geometric profiles with complicated conductivity configurations. In this paper, we propose a methodology for designing a new class of thermal source illusion devices for achieving directed thermal diffusions with natural homogeneous media. The employments of the space rotations in the linear transformation processes allow the directed thermal diffusions to be independent of the geometric profiles, and the utilization of natural homogeneous media improve the feasibility. Four schemes, with fewer types of homogeneous media filling the functional regions, are demonstrated in transient states. The expected performances are observed in each scheme. The related performance are analyzed by comparing the thermal distribution characteristics and the illusion effectiveness on the measured lines. The findings obtained in this paper see applications in the development of directed diffusions with minimal thermal loss, used in novel “multi-beam” thermal generation, thermal lenses, solar receivers, and waveguide.


Introduction
The theories of transformation optics (TO) [1] and conformal mapping [2] provide a wide platform for manipulating waves/fluxes under pre-designed methods in various fields, such as electromagnetics [3], optics [4], mechanics [5], acoustics [6], and thermodynamics [7]. Many significant achievements in the above fields, such as the novel invisible cloaks, efficient concentrators, lenses, illusions, etc., have been promoted by employing topological structures which contain form invariances of the governing functions following coordinate transformations. In this case, the tailored topological structures were achieved with near-zero parameters (such as near-zero relative permittivity, relative permeability, and relative thermal conductivity) to reveal the novel manipulative behaviors [8]. Using the pre-designed spatial transformations and the well-established effective medium approximation (EMA) [9], the electromagnetic and optic waves can be precisely controlled and directed through the tailored configurations of relevant material properties. In order to realize the parallel beam shifter and the beam splitter, a modified class of optical transformation was proposed, which employs finite, embedded coordinate transformations to steer the electromagnetic wave [10]. Further studies in embedded optical transformations were carried out, and conversion from the cylindrical waves to plane waves was achieved over a short range, using layered metamaterials [11]. To enhance the feasibility of realizing the tailored material properties, the discrete elements with

Theoretical Derivation of Conductivities Tensor
For realizing the proposed thermal source illusions, which are able to restore the outside thermal flux (temperature distributions) generated by the sources to make thermal illusions with varying characteristics of the existing sources as shown in Figure 1a,b, the schematic of the transformation process is shown in Figure 1c,d.

Theoretical Derivation of Conductivities Tensor
For realizing the proposed thermal source illusions, which are able to restore the outside thermal flux (temperature distributions) generated by the sources to make thermal illusions with varying characteristics of the existing sources as shown in Figure 1a,b, the schematic of the transformation process is shown in Figure 1c,d. For simplicity and generality, the entire system is restricted to 2D case, and one vertex of the inscribed polygon of the outer circle is fixed on the line of x = 0, i.e., the y direction, of the original and transformational domains. Obviously, the proposed source illusions are achieved with topological structures by mapping the inner circles onto the inscribed polygons of the outer circles. To demonstrate the transformation process, the entire system is divided into 2N parts through attaching each vertex and midpoint of the inscribed polygons with the origin O (0, 0). Furthermore, the separated parts are classified into two types, i.e., the parts with two vertices of (sin ((2n − 2) ⁄ ), cos ((2n − 2) ⁄ )), and ((sin ((2n − 1) ) ⁄ , cos ((2n − 1) ⁄ ))) are denoted as Type I, in which n ∈ [0, N], and the rests are Type II. Note that we adopt the rotation angles between the configurations and y directions in this paper. To keep the description of the transformation process simple, the transformations of regions in Type I are taken as an example. It is clear that the linear mapping function [13,14,23] can be employed to respectively map the lines of OAn and OCn into OAn' and OCn'. The linear mapping can be carried out by validating the relations of the original and transformational domains along the directions of the principle axes. That is, the relevant expansive/compressive coefficients should be solved through the following expressions:  For simplicity and generality, the entire system is restricted to 2D case, and one vertex of the inscribed polygon of the outer circle is fixed on the line of x = 0, i.e., the y direction, of the original and transformational domains. Obviously, the proposed source illusions are achieved with topological structures by mapping the inner circles onto the inscribed polygons of the outer circles. To demonstrate the transformation process, the entire system is divided into 2N parts through attaching each vertex and midpoint of the inscribed polygons with the origin O (0, 0). Furthermore, the separated parts are classified into two types, i.e., the parts with two vertices of (sin((2n − 2)π/N), cos((2n − 2)π/N)), and ((sin((2n − 1)π/N), cos((2n − 1)π/N))) are denoted as Type I, in which n ∈ [0, N], and the rests are Type II. Note that we adopt the rotation angles between the configurations and y directions in this paper. To keep the description of the transformation process simple, the transformations of regions in Type I are taken as an example. It is clear that the linear mapping function [13,14,23] can be employed to respectively map the lines of OA n and OC n into OA n and OC n . The linear mapping can be carried out by validating the relations of the original and transformational domains along the directions of the principle axes. That is, the relevant expansive/compressive coefficients should be solved through the following expressions: For the map of arc A n C n to segment A n C n , it can be directly realized in a simplified way by approximating the length of arc A n C n to segment A n C n . However, such a simplification would introduce major errors into the transformations, once the inscribed polygon with fewer edge numbers (N), i.e., the length of arc A n C n is not small enough compared with the radius of the inner regions (r) [14]. Hence, we further refine the configurations inside the part A n A n C n C n of Type I into m parts through the bisecting angle of A n OC n , in order to reduce the errors due to the simplifications. Each refined configuration is assigned a serial number p (A~D) , ranging from 1 to m. That is, the general vertices of an arbitrary element inside Type I can be obtained considering the positions of azimuths: where r and r 1 are the radii of the inner and outer circles. n denotes the serial number of each adjacent element inside Type I considering the positions of azimuth, ranging from 0 to N − 1. m are the multiples of equal dividing for the central angle and, p A and p C are the serial number of each characteristic component with equal central angle which is counting from the side of OA. According to the linear mapping function, the approximate inner triangle OA n,p A C n,p C can be mapped into the outer OA n,p A C n,p C with the following transformation processes upon the linearly expansive processes.
Based on the previous studies [16,17], the coefficients of a I,n , b I,n , d I,n , and e I,n are the elements of ∂x /∂x, ∂x /∂y, ∂y /∂x, and ∂y /∂y in the related Jacobi matrix used for combining such transformation processes. Owing to the 2D domain, the elements of ∂z /∂x and ∂z /∂y are zero, i.e., c I,n,p , and f I,n,p are independent of the changes in the z direction. Hence, the related transformed Jacobi matrix can be derived upon solving the coefficients of a I,n , b I,n , d I,n , and e I,n . The above equations can be transformed into matrix forms through using the correspondences on the both sides of the simultaneous equations in Equation (3).
Obviously, the coefficients of a I,n , b I,n , d I,n , and e I,n can be obtained through multiplying the inverse matrices of the corresponding parameters in the original region at the both sides of Equation (4). Considering the positions of characteristic points, the general expansive components in the refined elements of Type I can be achieved as follows.
For the elements in the parts of Type II, the conductivities can be obtained using the similar mapping function. The general positions of the vertices inside the elements of Type II are presented as follows: In Equation (7), p B and p D are the corresponding serial number of each characteristic component counting from the side of OB. Considering the mapping relations and the characteristic points, the expansive coefficients of the elements inside Type II can be derived by taking Equation (7) into Equation (1). Hence, the general expansive coefficients of Type II can be obtained.
x B n,p B y B n,p D n,p D 1 Taking Equations (5), (6), (8), and (9) together, the Jacobi matrix for realizing such 2D expansive actions are expressed as: where i denotes the types of the functional regions, i.e., Type I or II. Using the above Jacobi matrix, the conductivity components of each element, for achieving the proposed source illusions in the expansive regions, are obtained through solving the expressions as: κ = J·κ·J det(J) with det (J) = a i,n ·e i,n − b i,n ·f i,n , in which κ denotes the original conductivity between the inner and outer circles. Hence, the conductivity components can be expressed as follows: The components of the conductivities of Types I and II for the source illusions can be derived through Equations (11)-(13), by separately considering the transformations in Equations (5), (6), (8), and (9). Note that the values det (J) both in Types I and II are r 2 1 r 2 cos π N , owing to the linear mapping functions along the principle axes of the relevant curves. Hence, the components of conductivities of the proposed source illusions are independent of the range of the transformed regions. That is, the thermal illusions would theoretically be observed outside the outer circles, once the conductivity components in the characteristic direction are in accordance with the values calculated by Equations (11)-(13), without considering the sizes of the inner and outer circles.
For the sake of comparisons with those methods "isotropic cloaking formalism" and "anisotropic cloaking formalism" [37,38] employed in manipulating surface waves, the proposed method with local rotational and radial components can be used to achieve varied thermal behaviors by employing different geometric components. Owing to the linear maps between characteristic lines, the rotationally-symmetric profiles are not required, and the limit parameters can be significantly reduced.

Selections of Materials and Geometric Profiles of Thermal Source Illusions
The theoretically general components of the relevant thermal conductivity tensors of the elements inside Types I and II of the proposed source illusions are presented above. The next work is proposing the method for establishing such devices. As [26] pointed out, the relative conductivities for such composited systems, considering the deflections of heat flux [26][27][28], i.e., rotations of the local coordinate systems, could be obtained through alternative configurations of layered materials. Hence, the relative conductivities along the principle axes for functional regions of the source illusions can be deduced based on the effective medium theory.
In Equation (14), θ denotes relative rotation angle between the local configurations of the layered materials and corresponding principle axes (y axes). R(θ) and R(θ) are the unimodular rotation matrices. Furthermore, the relative conductivities along the principle axes can be obtained based on effective medium theory with the fact of θ = π 2 − 1 2 arctan 2(a i,n d i,n + b i,n e i,n ) Considering the findings from studies in bending heat flux with layered material configurations, the requirements of the relative conductivity components along the principle axes could be achieved by alternative configurations of negative and positive layers, with certain area fractions [29,32,34].
where, κ P and κ N are the relative conductivities for the positive and negative layers. γ P,l and γ N,l denote the composited fractions of the two layers in certain regions and the sum of γ P,l and γ N,l is 1.
In order to validate the derivations and demonstrate the source illusions, four schemes with varying profiles of the inscribed polygons, including triangle (N = 3), square (N = 4), pentagon (N = 5), and hexagon (N = 6), are proposed. As the illusion performances are independent of the physical dimensions, i.e., the radii of the inner circular spaces (r) where the thermal source located, and the function regions between the inner and outer circles (r 1 ), the radii of the outer (inscribed polygons) and inner circles are set as 0.1 m and 0.02 m, i.e., r 1 = 0.1 m and r = 0.02 m. The entire systems are placed at the center of the square plates with dimensions of 400 mm × 400 mm, filled with nickel steel (50% Ni) with a thermal conductivity of κ = 19.6 W·m −1 ·K −1 . In addition, the inner circles are used for diffusing the thermal flux generated by the sources. Hence, the nickel steel (50% Ni) is also employed to fabricate the inner circular regions as it has small effect on the restoring function. For the parts between the outer circle and corresponding inscribed polygons, nickel steel is also employed, as the related spaces are unchanged during the transformation process. Considering the characteristics of the designed distribution of thermal fields under related polygonal profiles, each area inside Types I and II of all the schemes is divided into two parts: 10% portions of the area at the positions of the geometric deformations are considered as one part, while the remaining 90% are another. The 10% separations are beneficial to restrict the regions with extreme conductivities in relatively small areas. In addition, this also helps by reducing the limit parameters at the geometric deformations, resulting in better approximate mapping. The following considerations are to select the appropriate positive and negative layers for creating the 10% and 90% parts of the functional regions. For the 10% part in each type, p A(B) and p C(D) should be 1, and the values of m should be 22.4, 15, 12.8, and 11.8 for the triangle, square, pentagon, and hexagon schemes, respectively. In the 10% part of Type I, the relative conductivities for the triangle and square schemes are far larger than natural positive mediums, and the relative conductivities for the pentagon and hexagon are 291.52 and 199.75 W·m −1 ·K −1 with nearly 100% fractions. In order to satisfy the requirements of the conductivities for the 10% parts in Type I, the pure copper with a conductivity of 398 W·m −1 ·K −1 is selected for the triangle, square, and pentagon schemes. Aluminum alloy (6063) with a conductivity of 201 W·m −1 ·K −1 is employed in the hexagon scheme. Furthermore, the adjacent 10% parts of the Types II are all filled with polydimethylsiloxane (PDMS) to create essential anisotropies coupled with those former 10% parts at the geometric deformations (near the vertices).
For the 90% part in each type, p A(B) is 2, and the values of m and p C(D) should also be 22.4, 15, 12.8 and 11.8 for the triangle, square, pentagon, and hexagon schemes, respectively. Taking the above newly defined points into Equations (5), (6) and (8), (9), the corresponding conductivity components and related fractions of the positive and negative layers can be achieved by Equations (14)- (16). Here, the negative and positive layers are respectively defined as layers A and B. Owing to the symmetries of the distribution of the functional regions, the configurations of the conductivities in the positive and negative layers are uniform. The desired conductivities of positive layers B and related fractions for the 90% part of Types I and II are determined as a function of the conductivities of the negative layers A, as shown in Figure 2. For the 90% part in each type, pA(B) is 2, and the values of m and pC(D) should also be 22.4, 15, 12.8 and 11.8 for the triangle, square, pentagon, and hexagon schemes, respectively. Taking the above newly defined points into Equations (5), (6) and (8), (9), the corresponding conductivity components and related fractions of the positive and negative layers can be achieved by Equations (14)  It can be seen that the conductivities of layers B increased with those of the negative layers. Meanwhile, the composited fractions of layers B simultaneously reduced in order to satisfy the demands calculated by Equations (15) and (16) based on the effective medium theory [26]. In order to obtain the appropriate conductivities in the negative and positive layers for realizing such source illusions, the polydimethylsiloxane (PDMS) with a conductivity of 0.15 W·m −1 ·K −1 is employed to completely fill in the layers A, owing to its low conductivity and preference of industry. For the conductivities of layers B, the values for the proposed triangle, square, pentagon, and hexagon schemes should be 44.015, 32.172, 28.233, and 26.232 W·m −1 ·K −1 , respectively. Considering the previous studies [29,32,34], the conductivities inside layers B could be achieved by alternately arranging negative and positive mediums, owing to the approximate values calculated by Equation (16)   It can be seen that the conductivities of layers B increased with those of the negative layers. Meanwhile, the composited fractions of layers B simultaneously reduced in order to satisfy the demands calculated by Equations (15) and (16) based on the effective medium theory [26]. In order to obtain the appropriate conductivities in the negative and positive layers for realizing such source illusions, the polydimethylsiloxane (PDMS) with a conductivity of 0.15 W·m −1 ·K −1 is employed to completely fill in the layers A, owing to its low conductivity and preference of industry. For the conductivities of layers B, the values for the proposed triangle, square, pentagon, and hexagon schemes should be 44.015, 32.172, 28.233, and 26.232 W·m −1 ·K −1 , respectively. Considering the previous studies [29,32,34], the conductivities inside layers B could be achieved by alternately arranging negative and positive mediums, owing to the approximate values calculated by Equation (16) and κ P,l = α P, m κ P,m + α N,m κ N,m (where α P,m and α N,m denote the area fractions of the positive and negative media composited in layers B, and κ P,m and κ N,m are the corresponding conductivities for the positive and negative media). In order to design the schemes with fewer kinds of media, the negative media composited in layers B are also PDMS. Hence, the total fractions of employed mediums inside this part (90% part) can be obtained.
In Equations (17) and (18), α P,m and α N,m are the fractions of the positive and PDMS composited in the each positive layer. β P,T and β N,T denote the total fractions of the positive and negative media (PDMS) employed in each 90% part of all types, respectively. Furthermore, the more balanced configurations of the positive media and PDMS, the optimal expected performances would occur, owing to the higher anisotropies inside the functional regions.
Considering the 90% part of the separated structures discussed above, the appropriate positive media, whose total composited fractions approach 40% in Types I and 50% in Types II, should be adopted for all the proposed schemes. Hence, die-casting aluminum A380 with a conductivity of 96.2 W·m −1 ·K −1 , ductile iron with a conductivity of 75 W·m −1 ·K −1 , and steel SAE 1010 with a conductivity of 59 W·m −1 ·K −1 are selected as the potential candidates for the proposed schemes. The total contents of the candidates in each scheme are illustrated in Figure 3. In order to employ fewer kinds of media in each scheme, only one positive medium, which approximately and simultaneously satisfied the fractional demands in Types I and II. Hence, die-casting aluminum A380 is selected for the triangle scheme with total fractions of 40.54% and 45.05% for Types I and II. Ductile iron is used for Types I and II of the square scheme with corresponding total fractions of 41.78% and 49.42%. Moreover, steel SAE is employed both in the pentagon and hexagon schemes, due to its most approximate fractions to the demands, i.e., 42.70% and 50.83% for Types I and II of the pentagon scheme, and 39.71% and 44.12% for Types I and II of the hexagon scheme. the candidates in each scheme are illustrated in Figure 3. In order to employ fewer kinds of media in each scheme, only one positive medium, which approximately and simultaneously satisfied the fractional demands in Types I and II. Hence, die-casting aluminum A380 is selected for the triangle scheme with total fractions of 40.54% and 45.05% for Types I and II. Ductile iron is used for Types I and II of the square scheme with corresponding total fractions of 41.78% and 49.42%. Moreover, steel SAE is employed both in the pentagon and hexagon schemes, due to its most approximate fractions to the demands, i.e., 42.70% and 50.83% for Types I and II of the pentagon scheme, and 39.71% and 44.12% for Types I and II of the hexagon scheme. Considering the thermal expanding effects of the functional regions (Types I and II), the wedge cell structures are employed for fabricating the positive and negative medium layers. Taking into account the selected media and related fractions, both in the 10% and 90% parts of Types I, one pure positive wedge layer (copper/aluminum alloy) is ordinarily employed in the 10% part for the proposed schemes. Five PDMS wedge layers and four positive medium (die-casting aluminum A380, Considering the thermal expanding effects of the functional regions (Types I and II), the wedge cell structures are employed for fabricating the positive and negative medium layers. Taking into account the selected media and related fractions, both in the 10% and 90% parts of Types I, one pure positive wedge layer (copper/aluminum alloy) is ordinarily employed in the 10% part for the proposed schemes. Five PDMS wedge layers and four positive medium (die-casting aluminum A380, ductile iron, and steel SAE 1010) wedge cells are embedded into the 90% part, as part of the total fractions. In analogy to Types II, one pure PDMS layer is embedded in the 10% part. Five positive medium layers and four PDMS wedge layers are employed in the rest regions with homologous fractions. The geometrical models for the proposed source illusion schemes are illustrated in Figure 4. Different from the previous contributions on electromagnetic meta-devices [39,40], the proposed linear maps methodology has been expanded into thermal fields with the considerations of local radial and angular components. Considering the general electromagnetic applications [41,42], it is believed that the designed illusive schemes based on the proposed methodology could open up an avenue to utilize the excited heat source efficiently in practical applications of the thermal field.  In order to validate and demonstrate the proposed schemes, numerical simulations based on the finite volume method are employed with ANSYS Fluent to solve the energy differential equation by adopting the second-order difference scheme. The boundaries of the calculated domains are thermally insulated, and the point sources are located at the center of the entire systems with constant temperatures Ts = 373 K. Furthermore, a scheme of pure nickel steel (50% Ni) plate with dimensions of 400 mm × 400 mm is established to make fair contrasts with the proposed schemes.

Demonstration of the Proposed Schemes and Discussions
Based on the transformation derivations and relevant thermal expanding effects, the temperature distributions of the contrasting pure nickel steel plate and the proposed source illusion schemes are obtained at t = 1000 s. Figure 5 illustrates the temperature distributions of the pure nickel In order to validate and demonstrate the proposed schemes, numerical simulations based on the finite volume method are employed with ANSYS Fluent to solve the energy differential equation by adopting the second-order difference scheme. The boundaries of the calculated domains are thermally insulated, and the point sources are located at the center of the entire systems with constant temperatures T s = 373 K. Furthermore, a scheme of pure nickel steel (50% Ni) plate with dimensions of 400 mm × 400 mm is established to make fair contrasts with the proposed schemes.

Demonstration of the Proposed Schemes and Discussions
Based on the transformation derivations and relevant thermal expanding effects, the temperature distributions of the contrasting pure nickel steel plate and the proposed source illusion schemes are obtained at t = 1000 s. Figure 5 illustrates the temperature distributions of the pure nickel steel plate. This indicates that the thermal energy generated by the point source diffused uniformly and forced the isothermal lines in the calculated domains, approaching circles as shown in Figure 5a, owing to the isotropies in the contrast scheme. In order to describe the characteristics of the temperature distributions and to make striking contrasts for between the following source illusion schemes, two characteristic lines, including the isothermal line of T = 302 K and the measured line of r = 0.1 m, are selected and illustrated in polar coordinates. As shown in Figure 5b, the r direction was defined as the radii of the locations of T = 302 K. It can be seen that the isothermal line formed in the circle profile, as the radii on the isothermal line of T = 302 K were 0.091 m, which were in accordance with that shown in Figure 5a. The temperature distributions on the line of r = 0.1 m were shown in Figure 5c. Obviously, the temperatures on the line of r = 0.1 m were all 300.6 K. That is, the thermal energy diffused uniformly along all the directions inside the system.

Temperature Distributions of the Proposed Thermal Source Illusions
The temperature distributions of the proposed thermal source illusions at t = 1000 s are shown in Figure 6. The temperature distributions outside the inscribed polygonal regions were shaped according to the pre-designed polygonal profiles. That is, the outer thermal fields were reshaped into designed polygons, and all the illusive performances were significant.

Temperature Distributions of the Proposed Thermal Source Illusions
The temperature distributions of the proposed thermal source illusions at t = 1000 s are shown in Figure 6. The temperature distributions outside the inscribed polygonal regions were shaped according to the pre-designed polygonal profiles. That is, the outer thermal fields were reshaped into designed polygons, and all the illusive performances were significant.
For the triangle scheme, the thermal energy diffused uniformly inside the inner circle region, i.e., the distributions of isothermal lines were shaped in circular profiles, owing to the inside isotropic material configurations and the point sources. Moreover, the distributions of the thermal energy were adjusted in the functional region, where the Types I and II were located, owing to the thermal expanding effects inside the functional regions. The isothermal lines were manipulated by the function elements of Types I and II, which forced non-uniform diffusions of the thermal energy, due to the alternating configurations of the positive and negative media.
With the outwards expansive diffusion, the thermal energies, which were located at the central regions corresponding to the outer boundaries, diffused into the outer domains of the triangle functional region faster than those located at the geometric deformations (near the vertices). That is, the outward diffusions of thermal energies were enhanced, given the reduced distances from the points on the boundaries to the origin O (0, 0). Furthermore, the refined configurations of the 10% fractional part at the geometric deformations promoted the thermal diffusions in the wedge elements approaching the vertices, owing to the high conductivities of the selected media. Hence, the profiles of the isothermal lines varied violently near the vertices of the outer triangle, which optimized distributions of the thermal fields outside the functional regions (outer circle) and reshaped profiles of the isothermal lines into triangular distributions. Thus, thermal source illusions occurred in the outer domains, which hid the original characteristics of the sources inside the central circle regions, and changed those into the pre-designed thermal distribution characteristics. Similar to the other inscribed polygon schemes, the distributions of the thermal fields were reshaped, by virtue of the functional wedge elements inside the inscribed polygons. Moreover, the profiles of the isothermal lines, which approached to the polygon boundaries, would better approximate the shapes of the inscribed polygons with the increasing side numbers of the polygons. That is, the more similar temperatures outside the outer circles were observed along the corresponding profiles of the inscribed polygons, owing to the reduced areas of the domains between the outer circles and the inscribed polygons.

Temperature Distributions of the Proposed Thermal Source Illusions
The temperature distributions of the proposed thermal source illusions at t = 1000 s are shown in Figure 6. The temperature distributions outside the inscribed polygonal regions were shaped according to the pre-designed polygonal profiles. That is, the outer thermal fields were reshaped into designed polygons, and all the illusive performances were significant. For the triangle scheme, the thermal energy diffused uniformly inside the inner circle region, i.e., the distributions of isothermal lines were shaped in circular profiles, owing to the inside isotropic material configurations and the point sources. Moreover, the distributions of the thermal energy were adjusted in the functional region, where the Types I and II were located, owing to the thermal expanding effects inside the functional regions. The isothermal lines were manipulated by the function elements of Types I and II, which forced non-uniform diffusions of the thermal energy, due to the alternating configurations of the positive and negative media.
With the outwards expansive diffusion, the thermal energies, which were located at the central regions corresponding to the outer boundaries, diffused into the outer domains of the triangle functional region faster than those located at the geometric deformations (near the vertices). That is, the outward diffusions of thermal energies were enhanced, given the reduced distances from the

Characteristics of Temperature Distribution
In order to further investigate the characteristics of the temperature distributions of the proposed schemes, the isothermal lines of T = 302 K and the measured line of r = 0.1 m were also selected. The characteristics of the temperature distributions of the isothermal lines of T = 302 K for the proposed schemes are illustrated in Figure 7. Obviously, the distribution characteristics of the thermal fields changed violently compared with that of the bare plate scheme shown in Figure 4b. As illustrated in Figure 6a, the radii of the locations of T = 302 K drastically fluctuated in the ranges of θ = 0 • -120 • , 120 • -240 • , and 240 • -360 • , and centered on θ = 60 • , 180 • , and 300 • , i.e., the ranges of the changing directions corresponded to the boundaries of the inscribed triangle. The radii of the locations rapidly decreased in the range of 0 • -60 • and then increased in the range of 60 • -120 • , due to the thermal flux manipulations of the functional elements inside the inscribed polygons. This indicates that the distributions of thermal fields were regulated following the profile of the pre-designed triangle, as the radii of the locations of T = 302 K varied gradually following the trends of the geometrical structure corresponding to the azimuths. In addition, the maximums of the locational radii appeared at the azimuths of θ = 0 • , 120 • , and 240 • , while the minimums were observed at θ = 60 • , 180 • , and 300 • . Thus, the profiles of the isothermal line of T = 302 K were approximately shaped in the triangle profile upon the above findings. However, the changing trends of the radii approached the azimuths of θ = 0 • , 120 • , and 240 • were larger than those of the boundaries of the inscribed triangle, which followed the function of r l = (r 1 cos π N )/ sin θ + π 2N , where r l denotes the radius on the boundary of the inscribed triangle. That is, the nature of thermal diffusions led to the local expansion of energy distributions with pure nickel steel inside the domains between the inscribed triangle and the outer circle. Similar to the other schemes, the maximums occurred, where the azimuths approached the vertices of the inscribed polygons, and the minimums were observed at the azimuths where the symmetry axes were located, i.e., θ = 45 • +nπ/4, 36 • + nπ/5, and 30 • + nπ/6 for the square, pentagon, and hexagon schemes (n ranges from 0 to N − 1). Furthermore, the inner regions inside the isothermal lines were enlarged with the increasing side counts of the inscribed polygons, which forced boundaries to approach the outer circles caused by the decreasing central angles. Meanwhile, the above phenomena also meant that the thermal energy distributions were expanded by the functional elements compared with the contrast scheme shown in Figure 5b, i.e., the minimum radii of the isothermal lines of T = 302 K were larger than those in Figure 5b.  Figure 8 provides the temperature distributions of the measured lines of r = 0.1 m (the outer circles) for the proposed schemes. It can be seen that the highest temperatures were observed at the positions of the geometric deformations of the inscribed polygons for all the proposed schemes, i.e., the vertices of the inscribed polygons located at the azimuths of θ = 2nπ/N. The lowest temperatures appeared at the azimuths of θ = (n + 1)·π/N in the proposed schemes. In addition, the lowest temperatures on the measured lines increased with the increasing side numbers. That is, the distributions of thermal energies were expanded outwards, which were in accordance with those illustrated in Figure 7. Furthermore, the temperatures first rapidly decreased and then increased at the azimuths ranging from 2nπ/N to 2(n + 1)·π/N, and centered on the azimuths of θ = (n + 1)·π/N, where the lowest temperatures were observed. Moreover, the changing trends of the temperatures increased with the increasing side counts, owing to the decreasing central angles of the inscribed polygons. In general, the temperatures near the vertices were significantly higher and the associated isothermal lines were forced to run parallel to the boundaries of the inscribed polygons, benefiting from the optimal conductivities in the 10% parts of the functional domains and the directed diffusions processes inside the 90% parts. Hence, the profiles of the isothermal lines outside the outer circle domains were approximately shaped in pre-designed profiles, and performed illusions for the actual sources through regulating the outside distributions of the thermal fields.
To further examine the illusion effectiveness of the proposed schemes, another seven measured lines of r = 0.04, 0.06, 0.08, 0.12, 0.14, 0.16, and 0.18 m were selected to calculate the temperature deformations [30] between the proposed schemes and the contrast bare plate, i.e., TD = |Tsc − Tpl|, where Tsc and Tpl are the temperatures on the measured lines of the proposed schemes and bare plate,  for the proposed schemes. It can be seen that the highest temperatures were observed at the positions of the geometric deformations of the inscribed polygons for all the proposed schemes, i.e., the vertices of the inscribed polygons located at the azimuths of θ = 2nπ/N. The lowest temperatures appeared at the azimuths of θ = (n + 1)·π/N in the proposed schemes. In addition, the lowest temperatures on the measured lines increased with the increasing side numbers. That is, the distributions of thermal energies were expanded outwards, which were in accordance with those illustrated in Figure 7. Furthermore, the temperatures first rapidly decreased and then increased at the azimuths ranging from 2nπ/N to 2(n + 1)·π/N, and centered on the azimuths of θ = (n + 1)·π/N, where the lowest temperatures were observed. Moreover, the changing trends of the temperatures increased with the increasing side counts, owing to the decreasing central angles of the inscribed polygons. In general, the temperatures near the vertices were significantly higher and the associated isothermal lines were forced to run parallel to the boundaries of the inscribed polygons, benefiting from the optimal conductivities in the 10% parts of the functional domains and the directed diffusions processes inside the 90% parts. Hence, the profiles of the isothermal lines outside the outer circle domains were approximately shaped in pre-designed profiles, and performed illusions for the actual sources through regulating the outside distributions of the thermal fields. In Equation (19), Nm is the numbers of the measured points on related measured lines. The standard deviations of the temperature deformations on the measured lines are shown in Figure 9. Note that the larger standard deviations indicate higher discreteness of the temperature deformations between observed values and their averages. Consequently, larger deformations would appear on corresponding azimuths of the proposed schemes, owing to the symmetries illustrated in Figure 5. Hence, more significant performances of illusions would be observed. As Figure 9 illustrates, the standard deviation on each measured line was 0 for the contrast bare plate, which contributed to shape the profiles of the temperature distributions as circles. For the proposed schemes, the standard deviations on each measured line increased with decreasing side numbers, which meant that the more significant illusions were achieved in the schemes with fewer diffusion elements inside the functional regions. For the independent schemes, the trends of the standard deviations as a function of the radii of the measured lines were similar, i.e., first increased in the functional regions and then decreased in the outer regions. In addition, the more approximate deformations would be observed in the schemes with more diffusion elements, i.e., more side counts. Hence, the more approximate temperature distributions were observed on the related azimuths with fewer side counts, due to the reduced central angles, which were in accordance with those shown in Figure 5 and validated in the electromagnetic field [14]. For the thermal diffusions in the outer spaces, the standard deviations decreased as the radii increased, due to the natural characteristics of the diffusions inside the isotropic medium of the outer spaces with thermal losses. Furthermore, the thermal fluxes were able to transfer efficiently, which could be orthogonal to the pre-designed boundaries of the inscribed polygons in the outer domains, once the media inside the outer spaces were compiled following certain conditions. To further examine the illusion effectiveness of the proposed schemes, another seven measured lines of r = 0.04, 0.06, 0.08, 0.12, 0.14, 0.16, and 0.18 m were selected to calculate the temperature deformations [30] between the proposed schemes and the contrast bare plate, i.e., T D = |T sc − T pl |, where T sc and T pl are the temperatures on the measured lines of the proposed schemes and bare plate, respectively. The standard deviations of the temperature deformations on the measured lines, used to express the effectiveness, can be written as: In Equation (19), N m is the numbers of the measured points on related measured lines. The standard deviations of the temperature deformations on the measured lines are shown in Figure 9. Note that the larger standard deviations indicate higher discreteness of the temperature deformations between observed values and their averages. Consequently, larger deformations would appear on corresponding azimuths of the proposed schemes, owing to the symmetries illustrated in Figure 5. Hence, more significant performances of illusions would be observed. As Figure 9 illustrates, the standard deviation on each measured line was 0 for the contrast bare plate, which contributed to shape the profiles of the temperature distributions as circles. For the proposed schemes, the standard deviations on each measured line increased with decreasing side numbers, which meant that the more significant illusions were achieved in the schemes with fewer diffusion elements inside the functional regions. For the independent schemes, the trends of the standard deviations as a function of the radii of the measured lines were similar, i.e., first increased in the functional regions and then decreased in the outer regions. In addition, the more approximate deformations would be observed in the schemes with more diffusion elements, i.e., more side counts. Hence, the more approximate temperature distributions were observed on the related azimuths with fewer side counts, due to the reduced central angles, which were in accordance with those shown in Figure 5 and validated in the electromagnetic field [14]. For the thermal diffusions in the outer spaces, the standard deviations decreased as the radii increased, due to the natural characteristics of the diffusions inside the isotropic medium of the outer spaces with thermal losses. Furthermore, the thermal fluxes were able to transfer efficiently, which could be orthogonal to the pre-designed boundaries of the inscribed polygons in the outer domains, once the media inside the outer spaces were compiled following certain conditions.

Conclusions
In general, a methodology for generally achieving directed thermal diffusions with thermal source illusions is proposed by setting several functional elements with the linear mapping function. The employment of space rotations can promote the fabrications of arbitrary thermal illusions in global transformation domains considering the geometric deformations on any azimuth. Note that such a methodology could perform the best results, once the length of arc of is small enough compared with the radius in the original domain, due to the geometric approximations in the design process. The details of the appropriate natural material selections with relevant fractions in varied illusion schemes have been investigated based on the effective medium theory. Fewer kinds of natural homogenous media for each proposed scheme, including the triangular, square, pentagonal, and hexagonal, have been employed to avoid the extreme parameters and improve the feasibility. The expected illusive performances and noticeably reshaped profiles of the isothermal lines have been observed in transient states. That is, the directed thermal diffusions through homogenous medium devices were validated, owing to the thermal-expansive effects inside the functional types. Furthermore, the performances and effectiveness of the source illusions have been demonstrated by characteristics of the measured lines and the standard deviations of the temperature deformations. The contrasts indicated that more approximate profiles to inscribed polygons of the isothermal lines, i.e., the better directed thermal diffusions, would be observed with more polygonal boundaries. As a theoretical and analytical work, the details of the novel design methodology have been provided and demonstrated by numerical models, some further worthy experimental works considering the practical conditions and environments in applications can be implemented to improve the efficiency of the heat transfer and thermal source utilizations in excited thermal techniques. The findings obtained in this paper can be used to design novel potential TO devices, including lens, source illusion, and diffusion in thermal field. Furthermore, the linear mapping function considering the changes of azimuths employed in this paper can be used to force the fluxes in the Laplace fields to efficiently transfer along certain directions in order to realize the high-gain "N-beam" diffusions [14].

Conclusions
In general, a methodology for generally achieving directed thermal diffusions with thermal source illusions is proposed by setting several functional elements with the linear mapping function. The employment of space rotations can promote the fabrications of arbitrary thermal illusions in global transformation domains considering the geometric deformations on any azimuth. Note that such a methodology could perform the best results, once the length of arc of is small enough compared with the radius in the original domain, due to the geometric approximations in the design process. The details of the appropriate natural material selections with relevant fractions in varied illusion schemes have been investigated based on the effective medium theory. Fewer kinds of natural homogenous media for each proposed scheme, including the triangular, square, pentagonal, and hexagonal, have been employed to avoid the extreme parameters and improve the feasibility. The expected illusive performances and noticeably reshaped profiles of the isothermal lines have been observed in transient states. That is, the directed thermal diffusions through homogenous medium devices were validated, owing to the thermal-expansive effects inside the functional types. Furthermore, the performances and effectiveness of the source illusions have been demonstrated by characteristics of the measured lines and the standard deviations of the temperature deformations. The contrasts indicated that more approximate profiles to inscribed polygons of the isothermal lines, i.e., the better directed thermal diffusions, would be observed with more polygonal boundaries. As a theoretical and analytical work, the details of the novel design methodology have been provided and demonstrated by numerical models, some further worthy experimental works considering the practical conditions and environments in applications can be implemented to improve the efficiency of the heat transfer and thermal source utilizations in excited thermal techniques. The findings obtained in this paper can be used to design novel potential TO devices, including lens, source illusion, and diffusion in thermal field. Furthermore, the linear mapping function considering the changes of azimuths employed in this paper can be used to force the fluxes in the Laplace fields to efficiently transfer along certain directions in order to realize the high-gain "N-beam" diffusions [14].
Acknowledgments: This work is supported by the National Natural Science Foundation of China (grant Nos. 51776050 and 51536001).
Author Contributions: Guo-Qiang Xu and Hao-Chun Zhang conceived and designed the research; Guo-Qiang Xu and Liang Jin conducted the calculation and analyzed the data; Yan Jin contributed analysis tools; and Guo-Qiang Xu wrote the paper, with the assistance of Hao-Chun Zhang.

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

Appendix A
In this part, the details of the employed numerical models presented in the section of "Selections of Materials and Geometric Profiles of Thermal Source Illusions" are further summarized and shown in Table A1. In this part, we would like to make a comparison among calculated values of effective medium approximation (EMA) [9,23,43], the classically alternative layer configurations [19] without the reserved 10% area portions, and the proposed methodology, to further validate and support results. The isothermal lines of T = 302 K outside the devices of the above contrasts are selected to demonstrate the radial distribution characteristics of the related points with 302 K. The distributions of the isothermal lines of T = 302 K in each contrast are shown in Figure A1. As can be seen from Figure A1, the radii of the points on the isotherms observed from the proposed schemes were quite approximate to the ideal values calculated by assuming EMA inside the functional regions. Furthermore, the peak values (near geometric deformations of the functional regions) of the ideal schemes and proposed schemes were close, owing to the enhanced heat transfer processes by the reserved 10% area portions. However, the valley values (near the centers of the corresponding boundaries of the functional regions) of the proposed schemes were higher than those ideal ones, due to the practical heat diffusions inside alternative medium layers. Moreover, the related values of the contrast schemes without the reserved enhancing area portions, i.e., employing the classical medium configurations [26][27][28][29][30][31][32][33][34], were also presented above. As it indicates, the classical medium configuration scheme would also perform the approximate characteristics to the ideal schemes. However, such approximations were weaker than the proposed schemes, as the peak values could not approach the ideal ones, and the valley values were higher than those of the proposed schemes. That is, the thermal energy inside the central medium cells diffused faster than that in the cells near the geometric deformations, due to the smaller heat transfer distances of the central cells. Hence, the profiles of the isotherms T = 302 K of the proposed schemes would be more approximate to the ideal EMA ones, compared with the classical medium configuration schemes without the reserved thermal enhancing area portions.

Appendix B.2. Independence Analysis of the Presented Schemes
The independence analysis of the presented schemes based on the proposed methodology is operated in this part. In order to complete the independency analysis, four measured points, i.e., A (0.15, 0), B (−0.05, 0), C (0, 0.1), and D (0, −0.1), are selected in all the schemes. Furthermore, five types of meshes with different grid numbers were selected for each mentioned scheme and are shown in Table A2. As can be seen from Figure A1, the radii of the points on the isotherms observed from the proposed schemes were quite approximate to the ideal values calculated by assuming EMA inside the functional regions. Furthermore, the peak values (near geometric deformations of the functional regions) of the ideal schemes and proposed schemes were close, owing to the enhanced heat transfer processes by the reserved 10% area portions. However, the valley values (near the centers of the corresponding boundaries of the functional regions) of the proposed schemes were higher than those ideal ones, due to the practical heat diffusions inside alternative medium layers. Moreover, the related values of the contrast schemes without the reserved enhancing area portions, i.e., employing the classical medium configurations [26][27][28][29][30][31][32][33][34], were also presented above. As it indicates, the classical medium configuration scheme would also perform the approximate characteristics to the ideal schemes. However, such approximations were weaker than the proposed schemes, as the peak values could not approach the ideal ones, and the valley values were higher than those of the proposed schemes. That is, the thermal energy inside the central medium cells diffused faster than that in the cells near the geometric deformations, due to the smaller heat transfer distances of the central cells. Hence, the profiles of the isotherms T = 302 K of the proposed schemes would be more approximate to the ideal EMA ones, compared with the classical medium configuration schemes without the reserved thermal enhancing area portions.

Appendix B.2 Independence Analysis of the Presented Schemes
The independence analysis of the presented schemes based on the proposed methodology is operated in this part. In order to complete the independency analysis, four measured points, i.e., A (0.15, 0), B (−0.05, 0), C (0, 0.1), and D (0, −0.1), are selected in all the schemes. Furthermore, five types of meshes with different grid numbers were selected for each mentioned scheme and are shown in Table A2. The values of the five measured points in each scheme with different meshes are illustrated in Figure A2. It can be seen that the temperatures of each measured point had little fluctuations with the increasing grid numbers. That is, the observed results for the schemes were independent of the meshes, i.e., the results are accurate to support the novelties. The values of the five measured points in each scheme with different meshes are illustrated in Figure A2. It can be seen that the temperatures of each measured point had little fluctuations with the increasing grid numbers. That is, the observed results for the schemes were independent of the meshes, i.e., the results are accurate to support the novelties.