Received: 03 Oct 2019
Revised: 13 Dec 2019
Accepted: 16 Dec 2019
Published online: 17 Dec 2019

Anomalous Transient Heat Conduction in Fractal Metamaterials

Wuxi Lin, Shengpeng Huang and Jie Ren*
Center for Phononics and Thermal Energy Science, China-EU Joint Center for Nanophononics, Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology, School of Physics Sciences and Engineering, Tongji University, Shanghai 200092, China 

*Corresponding address: 



Transient dynamics of heat conduction in isotropic fractal metamaterials is investigated. By using the Laplacian operator in non-integer dimension, we analytically and numerically study the effect of fractal dimensionality on the evolution of the temperature profile, heat flux and excess energy under certain initial and boundary conditions. Particularly, with randomly distributed absorbing heat sinks in the fractal metamaterials, we obtain an anomalous non-exponential decay behavior of the heat pulse diffusion. and an optimal dimension for efficient heat absorption as a function of sink concentrations. Our results may have potential applications in controlling transient heat conduction in fractal media, which will be ubiquitous as porous, composite, networked, hierarchical meta-materials.


Table of Content

We developed a theory to study transient heat conduction in fractal media, which are ubiquitous as porous, composite, hierarchical-networked metamaterials.





Keywords: Anomalous heat diffusion; Fractal metamaterials;  Non-integer dimensional; Fractional  Laplacian operator 

1. Introduction

Heat problem has attracted lots of attention in meta-structures in the past decade.1-4 Particularly, heat propagation as well as mass or excitation diffusion in fractal media are of great importance in our everyday life.5 This is because fractal meta-media describe the porous, composite, networked, hierarchical metamaterials, inwhich a part of the structure resembles larger entities or the whole structure. Such self-affine structural patterns are ubiquitous in many fields such as physics, material science and life science. For example, in branching artery network,6 in photosynthesis,7 and in bones,8 the fractal media of statistical self-similarity is quite involved. Therefore, fractal metamaterials have attracted much attention across diverse research fields, ranging from mechanics9,10 and elastics,11 to acoustics12-14 and optics.15-17 Different from those research that focused on effects of fractal structure on wave dynamics, heat conduction and transfer in fractal metamaterials is in general a diffusion process. 

Effective thermal conductivity is often used to characterize the heat diffusion in fractal media of fractal dimensions. Various methods have been developed to investigate the effective thermal conductivity of porous media,18-25 that may be even used to build a thermal diode.26 Pitchumani27 applied fractal theory in the research of the effective thermal conductivity of unidirectional fibrous composites. Yu and Cheng28 developeda fractal model to calculate the effective thermal conductivity of mono- and bi-dispersed porous media. Using thermal-electrical analogy, Yu29 and Kou30 presented fractal models and fractal analysis of effective thermal conductivityof composites with embedded fractal-like tree networks and saturated fractal porous media, respectively. 

To study fractal media of fractal dimensions, two methods are often applied: one is fractional calculus31,32 and the other one is calculus in fractional dimension space. In fractional calculus, factorial is replaced by gamma function to expand the application scope of previous calculuses,33 where the differential and integral calculus of time are usually involved. While calculus in fractional dimension space pays more attention to the geometric property of the non-integer space. An axiomatic system is established first34 and then applied for excitons in fractional dimensional space,35 and sequentially generalized.36,37 Stillinger constructed the axiomatic basis for spaces with non-integer dimension, and gave the form of Laplacian operator in fractal space.34 Tarasov suggested a generalization of vector calculus for the case of non-integer dimensional space,37 and gave a solution of heat propagation in fractal pipe and rodunder cylindrical coordinate.38

Relevant experiments have also been carried to study heat conduction in fractal media. Rozanova-Pierrat investigated how the shape of a prefractal radiator may enhance global heat transfer at short time.39 Cervantes-Alvarez reported the thermal characterization of plate-like composite samples made of polyester resin and magnetite inclusions.40 


Fig. 1 Schematic illustration of the transient heat conduction in fractal meta-media, which can be porous, composite, networked materials, showing a self-affine pattern that a part of the structure resembles larger entities or the whole structure. The heat diffusivity is described by the Laplacian operator in fractional dimension.

In this paper, we use calculus in fractional dimension to describe transient heat conduction in isotropic fractal media under specified boundary and initial conditions. For simplicity, we consider spherically symmetric media space. The initial temperaturedistribution in the spherical media is arbitrary and the boundary temperature keeps constant. We thus analytically and numerically study the evolution of the temperature profile, heat flux and excess energy. In particular, we analytically obtain the non-exponential decay of the heat pulse diffusion in the fractal media with randomly distributed absorbing heat sinks. In this case an optimal dimension exists for faster heat absorption, which depends on the heat sink concentrations. 

However, we have to mention that the mathematical operators in fractional dimension space may not fully describe the real fractal media, because the mathematics here does not maintain the strict self-similarity at all scales. Nevertheless, it captures themain characteristic of fractal media in large scale in the sense of effective media. For this reason, calculus in fractional space can offer us important insights about the diffusion behavior in the fractal media. Moreover, when varying dimensions, we consider the case of constant mass density following the same spirit as previous. 37,40 Our results and formulas can be readily extended to the case of varying mass density due to different arrangements of composite units at different dimensions. 

2. Theory of the Fractal Heat Conduction

In this section, we will provide the basic theory of the heat conduction in a fractal material. Although Fourier’s law may break down at nanoscale,41 throughout this work the fractal structures of interest are beyond the mesoscopic level and the unit size of the fractal structure is above the micrometer, so that Fourier’s law is valid. We consider the heat conduction in homogenous isotropic media, which is described by the Fourier’s law: J ⃗=-κ∇T, with J ⃗ the heat flux vector, κ the thermal conductivity of the medium, T the temperature as a function of location. The continuous equation of energy conservation requires the in-out flux balance: 


where f(r ⃗,t) denotes the intensity of the internal heat source, c is the heat capacity, ρ is the medium’s mass density, S and V denotes the surface and volume of the medium, respectively. As such, ∇⋅q ⃗+ρc ∂T/∂t=f(r,t). Considering both the divergence theorem and Fourier’s law, we arrive at the differential form: 


where ψ(r ⃗,t)=f(r ⃗,t)/κ is the source term and D=κ/(cρ) has the physical meaning of heat diffusion coefficient. The Laplacian operator in the fractal dimension can be written as 


where the non-integer d_s≥1 denotes the fractional dimension. This is a very important result cited from F. H. Stilinger’s original work34 [see also, the Laplace-Beltrami operator]. Stilinger’s work is a very important achievement in the area of fractal geometry. He establishes an axiomatic system under the condition that the space’s dimension is a fraction and on the basement he derives expressions of many basic geometry quantities such as the length and the volume element. On the basement he provides the expression of the Laplacian operator in the fractional dimension space. His results are widely cited in the papers about the researches of the physical characteristics in fractal media. For example, his theory to solve the problem about the light harvesting in a fluctuating antenna,7 Tarasov has referred to his theory to solve the problem about heat transfer in fractal media,38 and Wei-Ping Zhong et al. have used his theory to solve the problem about spatiotemporal accessible solitons in fractional dimensions.42 
In view of the fact that the homogenous isotropic thermal medium is spherically symmetric, all physical quantities can be represented as scalar functions of the radius distance, such as T(r ⃗,t)=T(r,t),ψ(r ⃗,t)=ψ(r,t). In such fractal media, the Laplacian operator can be simplified to: 


Therefore, we can express the general equation for the heat condition in isotropic fractal media as: 



3. Results and Simulations 

In what follows, we are going to present analytical formulas and numerical results for two typical transient diffusion cases in fractal dimensions: A) Heat diffusion driven by fixed temperature bias and heat sources; and B) Transient pulsed heat diffusionwith random absorbing sinks. 

3.1. Fixed Temperature Boundaries and Heat Sources
Assume the fractal medium boundary is surrounded by heat sinks at the distance r=R with the fixed boundary condition T(r=R,t)=T0 and the initial temperature profile T(r,t=0)=μ(r), we can obtain the solution: 

where ξn is the n-th zero point of the ds/2-1 fractional order Bessel function
. The first part in the integration describes the multi-time-scale relaxation from the initial temperature profile to the fixed boundary temperature due to the heat diffusion. The second part in the integration denotes the temperature raising due to the heat source flux. 

Therefore, we can calculate the heat flux as a function of location by making the gradient J(r,t)=-κ∇T(r,t): 

As the same, the first flux results from the heat diffusion due to the temperature difference between the interior profile and the exterior fixed boundary condition, and the second one results from the inject heat flux due to the source term, respectively. 
      Accordingly, using the hyper-sphere volume and integrating, we can also calculate the medium’s excessenergy, , as


Clearly, the thermal relaxation process is described by the series of time decay factors . As such, in the long time limit, the characteristic time-scale τ_c of the heat conduction relaxation is governed by the first term with the smallest exponent ξ_1^2, which can be written as 
                                     .    (9)
Clearly, increasing d_s will decrease the characteristic time, which means larger dimension can promote the heat diffusion. Note ξ_1 (d_s), as the first zero of the Bessel function J_(d_s/2-1), is a function of the fractal dimension d_s and has the theorem43,44:

We can compare it with the Brownian random motion of a single free particle immersed in a thermalized environment in multi-dimension ds. The mean square displacement (MSD) increases linearly with time as ⟨x2⟩=2ds Dt. When replacing theMSD ⟨x^2⟩ by the square radius R^2, it is clear to see that the upper bound τ_c in Eq. (9) is exactly the time for the MSD of the Brownian particle’s random walk reaching the square radius. In other words, the diffusiondynamics of heat conduction under fixed temperature boundaries is faster than the random diffusion of Brownian particles. This is understandable, because the former one is essentially a nonequilibrium process driven by external fields, including both theheat sources and temperature boundary conditions with inherent thermal bias, while the latter process is essentially an equilibrium under a uniform thermalized environment with no temperature bias. Therefore, the latter Brownian diffusion is slower with a larger time scale as the upper bound of the time scale of former heat diffusion. 


Fig. 2 Transient heat conduction driven by temperature difference (a)(b)(c) or heat source (d)(e)(f). For conduction driven by temperature difference: (a) Plots of the temperature evolution T(r,t)-T0 at r=0.6 m; (b) Plots of the density evolution of the heat flow q(r,t) at r=1 m; (c) Plots of the excess energy decay ΔE(ds,t). Dimension ds increases from 1 to 3 with increment of 0.4 in the direction of the arrow. Parameters are 1/D=8.1×103 s/m2, R=1 m, κ=518.52 W/(mds-2)⋅K), c=4.2×103 J/(kg⋅K), ρ=103 kg/mds), μ(r)=300 K, T0=100 K, ψ(r,t)=0. For conduction driven by heat source: (d) Plots of the temperature evolution T(r,t)-T0 with r=0.6m; (e) Plots of the density evolution of the heat flow J(r,t) with r=1m; (f) Plots of the incremental energy ΔE(ds,t) Dimension ds increases from 1 to 3 with increment of 0.4 in the direction of the arrow, where μ(r)=T0=100 K, ψ(r,t)=103 K/m(ds-1). Other parameters remain the same. 

In the following, we will present the associated numerical calculations. Figs. 2(a)(b)(c) plot the transient heat conduction only driven by temperature difference in the absence of heat source ψ(r)=0. The media of larger dimension ds can promote the heat propagation, so as to dissipate energy faster [see Fig. 2(a)]. This leads to smaller temperature gradient so as smaller heat flux density on the boundary of the media, as shown in Fig. 2(b). As shown in Fig. 2(c), ΔE also decreases faster as ds increases. This is intuitively because higher dimension means each physical site has more neighbors so that the propagation has more paths to spread out, more efficient and easier. 

Figs. 2(d)(e)(f) plot the transient heat conduction only driven by heat sources ψ(r,t)≠0. In other words, this case assumes that at the initial moment, the temperature in the media is the same as that in the exterior boundary, which means the uniform heat source exists to persist μ(r)=T0. From Fig. 2(d), we see that although with a lower saturated temperature, the temperature saturates faster in media of larger dimensions ds. This indicates that the larger ds can promote the faster heat conduction. As ds increases, since the temperature saturates faster to a lower value, the temperature gradient on the boundary is smaller so that the heat flux density is smaller on the boundary of the media [see Fig. 2(e)]. The behaviour of ΔE is more complicated. From Eq. (8), we see that at a short time, the change rate of ΔE is proportional to . While in the long time limit, the saturated ΔE is proportional to . We used in Ref. [45].) Correspondingly, as shown in Fig. 2(f), the change rate of ΔE increases as d_s increases, although the saturated ΔE decreases. 


3.2. Pulsed Heat Diffusion with Random Absorbing Sinks

It is worth noting that many subjects, such as life science and material science, have came across a similar problem: how to get the temperature profile and energy evolution when an initial heat pulse is excited in a fractal medium. That is to say, at the initial moment, the temperature focused at an spot is much higher than the other part of the media, and may be described by a Dirac delta function. For example, when a cancer tissue has been hit by the focused gamma ray beams, a spot of the cancer tissue can be at a very high temperature. If we can get the temperature profile and energy evolution in the cancer tissue, it will help to adjust the beam intensity or other parameters to achieve the goal of killing the cancer tissue. (Here is a document aboutusing gamma ray knife to treat cancer.46) Heat pulse excitation is also applied to materials to detect their intrinsic thermal diffusion and conductivity properties. 

This kind of problem is a specific form of the general problem we have risen above but if we take the Dirac delta function directly into the formula Eq. (6), we may not get the final result simply. So we use different scenario and mathematical methods to deal with this problem: apply a heat pulse to the center of a fractal sphere of uniform temperature, which may dissipate if introducing the distributed absorbing heat sinks.

Assume the heat pulse takes the form of Gaussian distribution with Tp a high temperature. When the length scale a is small enough, the heat pulse may be rewritten as , as a Dirac delta function. Therefore, denoting T0 the initial ambient temperature and the fixed temperature of heat sinks surrounding the fractal sphere boundary, the initial condition and boundary condition are given by

The temperature profile evolution under a heat pulse can be solved in an analytical form by setting ψ(r,t)=0 in Eq. (6), expressed as 


The amplitude Cn is given by 


with an omitted factor for small a≪R. The heat flux profile evolution is obtained accordingly: 


By integrating the distribution of temperature profile over the whole volume of the fractal sphere, we obtain the excess energy dwelling in the fractal medium at time t: 

                   .    (13)

As expected, this solution represents the kinetics of the excess energy injected by the heat pulse that decays from the initial one , with a dwelling fraction f_E (t)=ΔE(t)/ΔE(0) at time t. 

Therefore, the characteristic decay time of the heat pulse’s energy can be calculated to a simple expression: 

                            .    (14)

Surprisingly, this decay time of the heat pulse’s energy conducting in the fractal medium is just equal to the diffusion time needed for the Brownian random walk in this medium to reach the distance R. 

We note this observation does not conflict with the previous one [see Eq. (9) and associate discussions] where the nonequilibrium conduction is faster than the equilibrium diffusion. Here, the decay time is for the heat pulse that was excitedapproximately as a Dirac delta function, which makes the nonequilibrium process well described by the linear response theory, see Ref. [47]. And as is well known, in the linear response the nonequilibrium properties have connections to the equilibrium properties, such as the equivalence between heat conduction and diffusion,47 the connection between nonequilibrium transport coefficients and equilibrium flux autocorrelations in Green-Kubo formula.48 Thus, it is reasonablethat the decay time of the heat pulse is equal to the diffusion time of a Brownian motion in the fractal medium. 

So far, we have discussed the free heat diffusion in an unperturbed way in a fractal ds-dimensional hyper-sphere of radius R (so as volume Vds with no additional absorbing boundaries inside the medium. In reality, the medium interior could have thermal radiation spots or heat sinks randomly distributed as absorbing boundaries of thermal energy. Therefore, we assume that the heat sink distribution follows Poisson statistics, so that the probability to obtain no absorbing sinks inside (but on the boundary of) volume  is given by ,, where C is the concentration of absorbing heat sinks. By taking into account all possible routes of the heat pulse diffusion, we calculate the mean excess energy of the heat pulseremaining in the medium by averaging the dynamics above [see Eq. (13)] over different realizations of , written as: 
                                    .    (15)
In the long time limit t→∞, the decay dynamics of heat pulse’s energy is governed by the first term of the series (contains ξ1 only), which with the slowest decay dominates the asymptotic behavior. Thus, we can apply the Laplace’s method, orsay, the saddle-point approximation, to obtain the asymptotic 
                               ,    (16)
and  are coefficients given by A, . This asymptotic form clearly shows a non-exponential decay behavior. More importantly, the time-dependent decay behavior depends mainly on two undetermined parameters: the dimension of the fractal medium d_s and the concentration of absorbing heat sinks C, with thermal diffusivity D just rescaling the time, which all can be fitted out from the time-dependent experimental measures. 

Similarly, taking into account all possible routes of the heat pulse diffusion, we calculate the mean decay time by averaging over the probability distribution of , as: 
                     .    (17)
The same, by measuring the average decay time for different heat sink concentration C, we can also fit out the fractal dimension ds and the thermal diffusivity D. It is worth noting that ¯τ has a nontrivial dependence on dimension for difference C. For example, only at small C, increasing ds will decrease the decay time, which means larger dimension can promote the heat diffusion; while at large C, things are reversed. 
Fig. 3 Heat pulse diffusion with random absorbing sinks in fractal media. (a)(b) Plots of ¯ΔE(t) from Eq. (16) for C=30/m(ds) and C=100/m(ds), respectively, show complicated dimension dependences. Dimension ds increases from 1 to 3 with increment of 0.4 indicated by the arrow direction; (c) Plots of the mean decay time ¯τ(ds) from Eq. (17) with C=2/m(ds) in upper dot black line, 3/m(ds) in middle solid red line, 4/m(ds) in lower dash blue blue; (d) Optimal dimension ds with minimal decay time for efficient heat absorption. T0=100 K, ψ(r,t)=0, Tp=10^3 K, a=0.01 m. Other parameters are the same as before. Note, we limit ds between 1 and 3, consideringreality.


Also, we present here the associated numerical calculations. For the case of pulsed heat diffusion, the temperature evolution T(r,t)-T0 and the density evolution of the heat flow J(r,t) show similar behaviors as those in Figs. 2(a)(b). The excess energy decay ΔE(ds,t) is faster as ds increases. Therefore, we do not show them here repeatedly. For the heat pulse with random absorbing sinks, we plot the numerical results of ¯ΔE(t) and ¯τ(ds) in Fig. 3, by applying sharp Gaussian pulse with Tp=103 K and a=0.01 m. Figs. 3(a)(b) show clearly the non-exponential behaviors of the mean excess energy decay ¯ΔE(t) for heat sink concentration C=30/mds and C=100/mds, respectively. They have very complicated dimension dependences, as we can see from Eq. (16). 

The mean decay time ¯τ, [see Eq. (17)], also has a nontrivial dependence on d_s, which is affected by the concentration C of absorbing heat sinks. When C is relatively large, ¯τ increases as ds increases,as the dash line in Fig. 3(c). When C is relatively small, ¯τ decreases as ds increases, as the dot line in Fig. 3(c). If C is set to be an intermediate value in the middle, then ¯τ is allowed to obtained a minimal value with an optimal dimension ds, where the heat absorption is more efficient, see the solid line in Fig. 3(c) with enlarged view in the inset. The optimal ds is depicted in Fig. 3(d), which indicates that at intermediate C, the optimal dimension for efficient heat absorption is closed to the solid line. At lower sink concentrations, the larger dimension (ds→3) the better, while at higher sink concentrations, the lower dimension (ds→1) the better. 


4. Discussions

As we have noted in the beginning, the fractal structures studied at present are beyond the mesoscopic level and the unit size of the fractal structure is above the micrometer instead of at nanoscale, so that Fourier’s law is valid. To include the non-Fourier cases, one can generalize the master equation 2 (without loss of generality, we omit the source term): with constant heat diffusion coefficient D=κ/(cρ) to a non-Markov diffusion equation with memory: 
The retarded function K(t) is called the memory kernel, which plays the important role in diffusion dynamics so that the future will not only depend on the present state but also on the history. Some special kernels will lead to familiar equations as follows: 

                      ;    (20)

A completely no memory kernel K(t)∼δ(t) leads to the normal diffusion process, while a full memory of the history K(t)∼Θ(t) gives the ballistic wave equation. In between, the decaying memory kernel leads to the telegraph equation that combines the diffusion and ballistic wave equation. By applying the fractional dimension Laplacian operator ,, we can generalize the present discussions to the general diffusion process with non-Markovian memory kernels. 

So far, we restricted ourself to the constant diffusion coefficient D=κ/(cρ). We tried to fix the diffusion constant to merely see the pure dimension effect from varying d_s. This kind of scheme also follows the treatment in Refs. [38,40], where they also consider the mass density ρ, thermal conductivity κ, and so as diffusion coefficient D as a constant, when dealing with the problem about heat conduction in fractal media. Nevertheless, we need to note that in general the properties κ,c,ρ will have rich dependences on the dimension ds. Once the explicit dimension dependence is known, we can replace those constants as functions of ds, to see more rich and complicated heat conduction behavior in fractional dimension. 

Moreover, given spatial dependent system parameters, the fractional dimensional diffusion equation with different memory kernels, can be straightforwardly applied to the transformation thermodynamics,1 to design the transformed thermal cloaking, camouflage and so on, with fractional dimensions and anomalous non-Fourier thermal behaviors.


5. Conclusion

In this paper, we have described the transient heat conduction in fractal media in the framework of calculus in fractional dimension space. We have studied the influence of dimension to the evolution of the temperature distribution, the density of the heat flux and the excess energy, and several examples are analyzed to illustrate the results. We have found that in general larger dimension can promote heat propagation, but may with a complicated dependence on the system parameters. A special case for the heat pulse has been considered. With randomly distributed heat sinks in the media, we have obtained a non-exponential decay behavior of the excess energy, and the time-dependent decay behavior depends mainly on two undetermined parameters: the dimension of the fractal medium ds, and the concentration of absorbing heat sinks C. At lower sink concentrations, the large dimension promotes the heat absorption, while at higher sink concentrations, the lower dimension the better. An optimal dimension forefficient heat absorption emerges for intermediate C. By experimentally measuring the time-dependent kinetics, one may fit out the dimension of the fractal medium ds, the concentration of absorbing heat sinks C, and the thermal diffusivity D. 

Our results may have implications in material science, life science and medical science to describe transient heat conduction in fractal media like porous media, living tissue and composite. We hope they can be used to guide the design of thermal device, controlling transient heat conduction in ubiquitous fractal media, such as porous, composite, networked materials. 

5.1. The Simple Derivation of the General Solution

The heat diffusion in the fractal medium is expressed as: 
Suppose that this problem’s boundary conditions and initial conditions are T(r,t=0)=μ(r),T(r=R,t)=T0. Separate this temperature in its mathematical form T=T1+T2+T0, for the T1 and T2, the following equations are satisfied respectively: 
The T1,T2 satisfy the following initial conditions and boundary conditions respectively:

We solve the first equation by using the variable separation , then we get these two eigenvalue equations: 
Because of the boundary condition, we can confirm that the eigenvalue-related parameter λ we introduce above, is 
where ξn is the nth zero point of the (ds/2-1) order Bessel function
. Solve all the equations, then we get the general solution of the first equation: 

For the second equation, we use impulse theorem to solve it. As people often do in solving this kind of typical mathematical physics problems, we suppose that: 
What the equation v1 satisfies is very similar to what the equation T1 satisfies, so we can easily write down its general solution. Then we make an integral, and can write down the the general solution of T1


Considering the initial conditions, we can finally write down the special solution of , with 



which is exactly the Eq. (6). Then we can calculate the heat flux by means of making a differential [see Eq. (7)], and calculate the energy by making an integral [see Eq. (8)]. 


5.2. The Saddle-Point Method to get the Asymptotic Result
In order to get the result of the ¯ΔE(t) [see Eq. (15)], we use a mathematical method called saddle-point method. Because the dominate factor that influences the asymptotic result of the integral is the first term of the series, so we can calculate only the first term, and use the result to approximate the accurate result. We can provide the concrete form of the integral formula: 

Then we suppose that , , then the integral (excluding the coefficient) can be written as . In this formula,


Suppose that
By using steepest descent method, we can know that if R is always a real number, the f(z) can be approximated by the beneath expression: 49 
In the present circumstances, z=1, then the integral’s approximate expression is . The R0 in the formulas is the zero point of the function  . In this problem it is . We can substitute this result in the approximate expression of the integral, then we can get: 
Then we substitute the coefficient and the concrete expression of α and β in the above formula, and we get the result Eq. (16) in the main text. 


Our work is supported by the National Natural Science Foundation of China with grant No.11935010, No. 11775159, the Fundamental Research Funds for the Central Universities, and the Opening Project of Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology. 

1.   C. Z. Fan, Y. Gao and J. P. Huang, Appl. Phys. Lett. , 2008, 92, 251907.
2.   D. K. Ma, X. Wan and N. Yang, Phys. Rev. B, 2018, 98, 245420.
3.   D. Ma, A. Arora, S. Deng, G. Xie, J. Shiomi and N. Yang, Mater. Today Phys., 2019, 8, 56-61.
4.   Y. Xiao, Q. Chen, D. Ma, N. Yang and Q. Hao, ES Mater. Manuf., 2019, 5, 2-18.
5.   B. Yu, Appl. Mech. Rev., 2008, 61(5), 050801. 
6.   V. V. Gafiychuk, I. A. Lubashevsky and B. Y. Datsko, Phys. Rev. E, 2005, 72, 051920.
7.   J. Chmeliov, G. Trinkunas, H. van Amerongen and L. Valkunas, J. Am. Chem. Soc., 2014, 136, 8963-8972.
8.   N. Reznikov, M. Bilton, L. Lari, M. M. Stevens and R. Kröger, Science , 2018, 360, eaao2189.
9.   D. Rayneau-Kirkhope, Y. Mao and R. Farr, Phys. Rev. Lett., 2012, 109, 204301. 
10.   L. R. Meza, A. J. Zelhofer, N. Clarke, A. J. Mateos, D. M. Kochmann and J. R. Greer, Pnas, 2015, 112, 11502.
11.   K. Billon, I. Zampetakis, F. Scarpa, M. Ouisse, E. Sadoulet-Reboul, M. Collet, A. Perriman and A. Hetherington, Compos. Struct., 2017, 160, 1042.
12.   G. Y. Song, Q. Cheng, B. Huang, H. Y. Dong and T. J. Cui, Appl. Phys. Lett., 2016, 109, 131901 (2016).
13.   M. Fellah, Z. E. A. Fellah, A. Berbiche, E. Ogam, F. G. Mitri and C. Depollier, Wave Motion., 2017, 72, 276. 10.1016/j.wavemoti.2017.04.006
14.   X. Zhao, G. Liu, C. Zhang, D. Xia and Z. Lu, Appl. Phys. Lett., 2018, 113, 074101. 10.1063/1.5038431
15.   F. Miyamaru, Y. Saito, M. W. Takeda, B. Hou, L. Liu, W. Wen and P. Sheng, Phys. Rev. B , 2018, 77, 045124. 10.1103/physrevb.77.045124
16.   H. X. Xu, G. M. Wang, M. Q. Qi, L. Li and T. J. Cui, Adv. Opt. Mater., 2013, 1, 495.
17.   D. Garoli, E. Calandrini, A. Bozzola, A. Toma, S. Cattarin, M. Ortolani and F. De Angelis, ACS Photonics, 2018, 5, 3408.
18.   K. Ramani and A. Vaidyanathan, J. Compos. Mater., 1995, 29, 1725. 
19.   M. R. Islam and A. Pramila, J. Compos. Mater., 1999, 33, 1699 (1999).
20.   K. Bakker, Int. J. Heat Mass Transfer, 1997, 40, 3503.
21.   C. R. Havis, G. P. Peterson and L. S. Fletcher, J. Thermophys., 1989, 3, 416. 
22.   J. F. Thovert, F. Wary and P. M. Adler, J. Appl. Phys., 1990, 68, 3872 (1990). 
23.   D. Veyret, S. Cioulachtjian, L. Tadrist and J. Pantaloni, ASME J. Heat Transfer., 1993, 115, 866.
24.   J. Y. Qian, Q. Li, K. Yu and Y. M. Xuan, Sci. China Ser. E, 2004, 47, 716.
25.   X. L. Huai, W. W. Wang and Z. G. Li, Appl. Therm. Eng., 2007, 27, 2815.
26.   W. Zhu, G. Wu, H. Chen H and J. Ren, Front. Energy Res., 2018, 6, 9.
27.   R. Pitchumani and S. C. Yao, J. Heat Transfer, 1991, 113, 788.
28.   B. Yu and P. Cheng, J. Thermophys. Heat Transfer, 2002, 16, 22.
29.   B. Yu, Phys. Rev. E, 2006, 73, 066302.
30.   J. Kou, F. Wu, H. Lu, Y. Xu and F. Song, Phys. Lett. A, 2009, 374, 62.
31.   Y. Povstenko and J. Klekot, J. Appl. Math. Comput. Mech., 2017, 16, 101. 
32.   L. Fila€shtinskii, T. V. Mukomel and T. A. Kirichok, J. Math. Sci., 2011, 178, 557.
33.   V. V. Uchaikin, Fractional Derivatives for Physicists and Engineers, Springer, Berlin (2013).
34.   F. H. Stillinger, J. Math. Phys., 1977, 18, 1224.
35.   X. F. He, Phys. Rev. B, 1991, 43, 2063.
36.  C. Palmer and P. N. Stavrinou, J. Phys. A: Math. Gen., 2004, 37, 6987.
37.   V. E. Tarasov, Phys. Lett. A., 2005, 341, 467; Commun. Nonlinear Sci. Numer. Simul., 2015, 20, 360.
38.   V. E. Tarasov, Int. J. Heat Mass Transfer, 2016, 93, 427. 
39.   A. Rozanova-Pierrat, D. S. Grebenkov and B. Sapoval, Phys. Rev. Lett., 2012, 108, 240602.
40.   F. Cervantes-Alvarez, J. J. Reyes-Salgado, V. Dossetti and J. L. Carrillo, J. Phys. D: Appl. Phys., 2014, 47, 235303. 
41.   N. B. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi and B. Li, Rev. Mod. Phys., 2012, 84, 1045. 
42.   W. Zhong, M. R. BeliÄ, B. A. Malomed, Y. Zhang and T. Huang, Phys. Rev. E, 2016, 94, 012216.
43.   R. Piessens, Math. Comp., 1984, 42, 195.
44.   Á. Elbert and P. D. Siafarikas, Canad. Math. Bull., 1999, 42, 56.
45.   G. N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge University Press, Page 502 (1995). 
46.   J. Bernier, E. J. Hall and A. Giaccia, Nat. Rev. Cancer, 2004, 4, 737.
47.   S. Liu, P. Hänggi, N. Li, J. Ren and B. Li, Phys. Rev. Lett., 2014, 112, 040601. 
48.   M. S. Green, J. Chem. Phys., 1954, 22, 398; R. Kubo, J. Phys. Soc. Jpn., 1957, 12, 570. 
49.   Z. X. Wang and D. R Guo, Special Functions, Peking University Press, 276 (2012).