Received: 03 Oct 2019

Revised: 13 Dec 2019

Accepted: 16 Dec 2019

Published online: 17 Dec 2019

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: Xonics@tongji.edu.cn

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.

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

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 mechanics^{9,10} and elastics,^{11} to acoustics^{12-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} Pitchumani^{27} applied fractal theory in the research of the effective thermal conductivity of unidirectional fibrous composites. Yu and Cheng^{28 }developeda fractal model to calculate the effective thermal conductivity of mono- and bi-dispersed porous media. Using thermal-electrical analogy, Yu^{29} and Kou^{30} 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 calculus^{31,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 first^{34} 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.

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:

(1)

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:

(2)

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

* * (3)

where the non-integer d_s≥1 denotes the fractional dimension. This is a very important result cited from F. H. Stilinger’s original work^{34} [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:

* * (4)

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

(5)

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.

Assume the fractal medium boundary is surrounded by heat sinks at the distance r=R with the fixed boundary condition T(r=R,t)=T

(6)

where ξ_{n} is the n-th zero point of the d_{s}/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):

(7)

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

(8)

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 theorem^{43,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 ⟨x^{2}⟩=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)-T_{0} 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(d_{s},t). Dimension d_{s} increases from 1 to 3 with increment of 0.4 in the direction of the arrow. Parameters are 1/D=8.1×10^{3} s/m^{2}, R=1 m, κ=518.52 W/(m^{ds}-2)⋅K), c=4.2×10^{3} J/(kg⋅K), ρ=10^{3} kg/m^{ds}), μ(r)=300 K, T_{0}=100 K, ψ(r,t)=0. For conduction driven by heat source: (d) Plots of the temperature evolution T(r,t)-T_{0} 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(d_{s},t) Dimension d_{s} increases from 1 to 3 with increment of 0.4 in the direction of the arrow, where μ(r)=T_{0}=100 K, ψ(r,t)=10^{3} 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 d_{s} 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 d_{s} 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)=T_{0}. From Fig. 2(d), we see that although with a lower saturated temperature, the temperature saturates faster in media of larger dimensions d_{s}. This indicates that the larger d_{s} can promote the faster heat conduction. As d_{s} 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

* * (10)

The amplitude C_{n} is given by

* * (11)

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

* * (12)

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 d_{s}-dimensional hyper-sphere of radius R (so as volume *V**d**s* 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)

where 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 d_{s} 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 d_{s} 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 d_{s} increases from 1 to 3 with increment of 0.4 indicated by the arrow direction; (c) Plots of the mean decay time ¯τ(d_{s}) 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 d_{s} with minimal decay time for efficient heat absorption. T_{0}=100 K, ψ(r,t)=0, T_{p}=10^3 K, a=0.01 m. Other parameters are the same as before. Note, we limit d_{s} 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)-T_{0} 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(d_{s},t) is faster as d_{s} 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 ¯τ(d_{s}) in Fig. 3, by applying sharp Gaussian pulse with T_{p}=10^{3} 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/m^{ds} and C=100/m^{ds}, 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 d_{s} increases,as the dash line in Fig. 3(c). When C is relatively small, ¯τ decreases as d_{s} 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 d_{s}, where the heat absorption is more efficient, see the solid line in Fig. 3(c) with enlarged view in the inset. The optimal d_{s} 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 (d_{s}→3) the better, while at higher sink concentrations, the lower dimension (d_{s}→1) the better.

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:

* * (18)

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:

(19)

; (20)

* *

* * (21)

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 d_{s}. Once the explicit dimension dependence is known, we can replace those constants as functions of d_{s}, 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.

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 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)=T

The T

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 ξ

(28)

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:

(29)

What the equation v_{1} satisfies is very similar to what the equation T_{1} satisfies, so we can easily write down its general solution. Then we make an integral, and can write down the the general solution of T_{1}:

* * * *

* *(30)

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:

* * (31)

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

Suppose that

* * (32)

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

* * (33)

In the present circumstances, z=1, then the integral’s approximate expression is . The R_{0} 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:

* * (34)

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.

**Acknowledgements**

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.

https://doi.org/10.1063/1.2951600

2. D. K. Ma, X. Wan and N. Yang, Phys. Rev. B, 2018, 98, 245420.

https://doi.org/10.1103/PhysRevB.98.245420

3. D. Ma, A. Arora, S. Deng, G. Xie, J. Shiomi and N. Yang, Mater. Today Phys., 2019, 8, 56-61.

https://doi.org/10.1016/j.mtphys.2019.01.002

4. Y. Xiao, Q. Chen, D. Ma, N. Yang and Q. Hao, ES Mater. Manuf., 2019, 5, 2-18.

http://doi.org.10.30919/esmm5f237

5. B. Yu, Appl. Mech. Rev., 2008, 61(5), 050801.

https://doi.org/10.1115/1.2955849

6. V. V. Gafiychuk, I. A. Lubashevsky and B. Y. Datsko, Phys. Rev. E, 2005, 72, 051920.

https://doi.org/10.1103/PhysRevE.72.051920

7. J. Chmeliov, G. Trinkunas, H. van Amerongen and L. Valkunas, J. Am. Chem. Soc., 2014, 136, 8963-8972.

https://doi.org/10.1021/ja5027858

8. N. Reznikov, M. Bilton, L. Lari, M. M. Stevens and R. Kröger, Science , 2018, 360, eaao2189.

https://doi.org/10.1126/science.aao2189

9. D. Rayneau-Kirkhope, Y. Mao and R. Farr, Phys. Rev. Lett., 2012, 109, 204301. https://doi.org/10.1103/PhysRevLett.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.

https://doi.org/10.1073/pnas.1509120112

11. K. Billon, I. Zampetakis, F. Scarpa, M. Ouisse, E. Sadoulet-Reboul, M. Collet, A. Perriman and A. Hetherington, Compos. Struct., 2017, 160, 1042.

http://doi.org/10.1016/j.compstruct.2016.10.121

12. G. Y. Song, Q. Cheng, B. Huang, H. Y. Dong and T. J. Cui, Appl. Phys. Lett., 2016, 109, 131901 (2016).

http://doi.org/10.1063/1.4963347

13. M. Fellah, Z. E. A. Fellah, A. Berbiche, E. Ogam, F. G. Mitri and C. Depollier, Wave Motion., 2017, 72, 276.

http://doi.org/ 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. http://doi.org/ 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.

http://doi.org/ 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.

http://doi.org/10.1002/adom.201300023

17. D. Garoli, E. Calandrini, A. Bozzola, A. Toma, S. Cattarin, M. Ortolani and F. De Angelis, ACS Photonics, 2018, 5, 3408. http://doi.org/10.1021/acsphotonics.8b00676

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. http://doi.org/10.1016/j.applthermaleng.2007.01.031

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. http://doi.org/10.1115/1.2911205

28. B. Yu and P. Cheng, J. Thermophys. Heat Transfer, 2002, 16, 22. http://doi.org/10.2514/2.6669

29. B. Yu, Phys. Rev. E, 2006, 73, 066302. http://doi.org/10.1103/PhysRevE.73.066302

30. J. Kou, F. Wu, H. Lu, Y. Xu and F. Song, Phys. Lett. A, 2009, 374, 62.

http://doi.org/10.1016/j.physleta.2009.10.015

31. Y. Povstenko and J. Klekot, J. Appl. Math. Comput. Mech., 2017, 16, 101.

32. L. Filashtinskii, T. V. Mukomel and T. A. Kirichok, J. Math. Sci., 2011, 178, 557.

http://doi.org/10.1007/s10958-011-0569-2

33. V. V. Uchaikin, Fractional Derivatives for Physicists and Engineers, Springer, Berlin (2013). http://doi.org/10.1007/978-3-642-33911-0

34. F. H. Stillinger, J. Math. Phys., 1977, 18, 1224. http://doi.org/10.1063/1.523395

35. X. F. He, Phys. Rev. B, 1991, 43, 2063. http://doi.org/10.1103/PhysRevB.43.2063

36. C. Palmer and P. N. Stavrinou, J. Phys. A: Math. Gen., 2004, 37, 6987.

http://doi.org/10.1088/0305-4470/37/27/009

37. V. E. Tarasov, Phys. Lett. A., 2005, 341, 467; Commun. Nonlinear Sci. Numer. Simul., 2015, 20, 360. http://doi.org/10.1016/j.physleta.2005.01.024 http://doi.org/10.1016/j.cnsns.2014.05.025

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. http://doi.org/10.1103/physrevlett.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. http://doi.org/10.1103/PhysRevE.94.012216

43. R. Piessens, Math. Comp., 1984, 42, 195. http://doi.org/10.2307/2007570

44. Á. Elbert and P. D. Siafarikas, Canad. Math. Bull., 1999, 42, 56.

http://doi.org/10.4153/CMB-1999-007-4

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. http://doi.org/10.1038/nrc1451

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).