Integro-Differential Equation for the Non-Equilibrium Thermal Response of Glass-Forming Materials: Analytical Solutions

An integro-differential equation describes the non-equilibrium thermal response of glass-forming substances with a dynamic (time-dependent) heat capacity to fast thermal perturbations. We found that this heat transfer problem could be solved analytically for a heat source with an arbitrary time dependence and different geometries. The method can be used to analyze the response to local thermal perturbations in glass-forming materials, as well as temperature fluctuations during subcritical crystal nucleation and decay. The results obtained can be useful for applications and a better understanding of the thermal properties of glass-forming materials, polymers, and nanocomposites.


Introduction
The dynamic heat capacity c dyn (t) of glass-forming liquids and glasses, considered as a function of time t, has been intensively studied since the pioneering work of Birge and Nagel published in 1985 [1]. However, much earlier, experiments on the dispersion and absorption of ultrasonic waves in polyatomic gases and liquids were explained by the relaxation of the specific heat in these substances. The experiments were comprehensively reviewed by Herzfeld and Litovitz in 1959 [2]. The relaxation of the specific heat of polyatomic gases and liquids is caused by a slow energy exchange between the external (translational) and internal (vibrational and rotational) degrees of freedom. Thus, the energy exchange in polyatomic gases is characterized by a limited set of characteristic relaxation times τ i [2]. The spectrum of relaxation times of glass-forming liquids and glasses is extremely broad, and it can be considered as a continuous spectrum . Naturally, this broad spectrum is observed not only for the relaxation of the dynamic heat capacity [3][4][5][6][7][18][19][20] but also for dielectric susceptibility [8][9][10][11][12][13][14][15][16][17][18], light scattering [23], and viscosity [8,[24][25][26]. Since the specific heat of glass-forming substances depends on time, it follows that the thermal response of these materials to a thermal perturbation at time t depends on the temperature distribution T t, → r in the system in previous times. This effect is especially significant for fast local thermal perturbations [27,28].
The time dependence of the heat capacity of glass-forming liquids and glasses leads to a non-equilibrium (non-Fourier) thermal response of these materials to fast thermal perturbations. This non-equilibrium thermal response can be described by the integrodifferential heat equation considered in [27][28][29]. This equation can be solved similarly to the second-kind Volterra integral equations. Volterra integral equations have several applications in many branches of science, technology, and industry. Viscoelasticity is one Symmetry 2021, 13,256 2 of 17 of the fields of physics where the Volterra equations are often used [25]. The applications of Volterra equations in renewable energy is an example of their use in industry [30][31][32]. Volterra integro-differential equations are usually difficult to solve analytically; therefore, approximate solutions and numerical methods are often used [33]. An interesting method for the numerical solution to nth-order integro-differential equations has been developed based on the integral mean value theorem [34]. However, to better understand the nature of physical phenomena, it is crucial to establish qualitative relationships between physical parameters and the ongoing physical processes. Thus, we focused on finding analytical solutions of the integro-differential heat equation that describe the non-equilibrium thermal response of glass-forming materials.
The heat transfer equation considered in this article was solved analytically and may be of interest for various applications. In previous articles, we found solutions for the equation with rectangular pulsed heat sources [28,29]. In this study, we focused on analytical solutions of the equation for a heat source with arbitrary time dependence. These solutions can be useful for studying local thermal perturbations in glass-forming materials, as well as temperature fluctuations during subcritical crystal nucleation and decay; these temperature fluctuations can have a significant effect on the kinetics of crystal nucleation in glass-forming materials.
The heat transfer problem with memory was analyzed in general terms by Miller in 1978 [35]. For an external heat flux represented by a smooth function, he proved the existence, uniqueness, and continuous dependence on parameters of the solution for heat transfer with memory. In this study, we considered a special case of integro-differential heat equations with kernels corresponding to glass-forming materials with dynamic heat capacity; these equations can be solved analytically. We focused on an external heat flux Φ t, → r acting on finite intervals in space and time. In this study, we restricted ourselves to considering the multiplicatively separable heat flux Φ → r F(t). We considered the heat flux as a continuous and piecewise smooth function. Thus, we assumed that both F(t) and Φ → r were continuous piecewise smooth functions. This is a sufficient condition for the absolute and uniform convergence of the Fourier series of the functions F(t) and Φ → r [36]. In addition, we focused on the dynamic behavior of substances with a dynamic heat capacity and restricted ourselves to considering only homogeneous boundary value problems.
The rest of the paper is organized as follows. Sections 2 and 3 formulate the heat transfer problem to be solved. Sections 4 and 5 discuss solutions for rectangular-pulsed, sinusoidal, and arbitrary heat sources that were considered for planar and spherical geometries. In Section 6, the effect of the relaxation time distribution on the thermal responses is discussed for different temperatures and materials. It is shown that the effect of the time dispersion of the dynamic heat capacity on the thermal response T(t, r) during local heating was significant. Examples of solutions T(t, r) for real glass-forming materials were considered. The temperature dependence of the relaxation-time distribution was taken into account.

Applicability of the Heat Equation with Dynamic Heat Capacity
In this study, we focused on the heat equation for non-metallic glass-forming materials at temperatures T outside the low-temperature range. Nonlocal effects of heat conduction [37] were not taken into account for (∂ ln(T)/∂x) −1 l ph , where l ph is the phonon mean-free-path. Furthermore, for small temperature changes, we considered the thermal parameters of materials independent of T. However, the temperature dependence of the relaxation time spectrum of the dynamic heat capacity c dyn (t) is discussed in Section 6. This relaxation time spectrum is very broad in glass-forming materials. Therefore, the time dispersion of the dynamic heat capacity is significant over a wide range of time scales. Conversely, the thermal conductivity λ of non-metallic glass-forming materials can be considered an equilibrium (time-independent) thermal parameter, at least for t > 1 ns and at Symmetry 2021, 13, 256 3 of 17 temperatures above the low-temperature range. In fact, thermal excitations associated with thermal conductivity come to local equilibrium much faster than in 1 ns in glass-forming liquids and glasses [38][39][40][41], and the phonon mean free path is about 1 nm or less [41][42][43][44][45][46][47]. Therefore, we focused on the diffusion-type Fourier heat conduction. However, we took into account the fact that the local heat absorption at a given moment of time t is determined not only by a change in the local temperature T t, → r at this moment but also at each previous moment. The temporal dispersion of the dynamic heat capacity c dyn (t) of glass-forming materials can be described within the framework of the linear response theory [6,48], similarly to the temporal dispersion of the dielectric constant [49].

Heat Equation with Dynamic Heat Capacity
Heat transfer in glass-forming materials with a dynamic heat capacity c dyn (t) can be described using the following integro-differential heat equation [27,28]: where T t, → r is the external heat flux. Equation (1) can be used for glass-forming substances at least for t > 1 ns and length scales greater than 1 nm.
Consider the problem with zero initial conditions Φ t, → r = 0 and T t, → r = 0 for t ≤ 0. Thus, from Equation (1) we obtain: Equation (2) can be solved if the dynamic heat capacity c dyn (t) is a given function. Usually, relaxation phenomena in glass-forming materials are described by the stretched exponential Kohlrausch relaxation function exp −(t/τ K ) β for β ∈ (0, 1], where β and the Kohlrausch relaxation time τ K characterize the relaxation time spectrum (for more details, see [50,51]). As a completely monotonic function, the stretched exponent can be represented as a continuous sum of exponentials (see Bernstein's theorem [52]). It should be noted that typical relaxation functions, such as exponential, stretched exponential, and power-law relaxation functions, are completely monotonic [53]. Thus, we assumed that the dynamic heat capacity c dyn (t) is a completely monotonic function of time; therefore, c dyn (t) can be represented as a continuous sum of exponentials, as seen in Equation (3): where c 0 and c in are the equilibrium and initial heat capacities, that is, c dyn (t) → c in as t → 0 and c dyn (t) → c 0 as t → ∞ . The distribution function H(τ 0 ) is normalized as follows: ∞ 0 H(τ 0 )dτ 0 = 1. However, for practical use, c dyn (t) can be given on a finite interval [τ min , τ max ] if this interval is wide enough (see [28] for details). The distribution function H(τ 0 ) can be found by using broadband heat capacity spectroscopy [19], as was done in [27][28][29] (see Example 3).
As a first step, consider the solution to Equation (2) for the auxiliary problem, which corresponds to the dynamic heat capacity c dyn (t) obeying the Debye relaxation law: where ε 0 = (c 0 − c in )/c 0 . Subsequently, the final solution for any given c dyn (t) can be obtained using the solution for an arbitrary positive τ 0 if H(τ 0 ), c 0 , and c in are specified. This final solution can be represented as a continuous sum of solutions depending on τ 0 and distributed according to the normalized distribution function H(τ 0 ). Using Equations (2) and (4), we obtain the following integro-differential equation with a difference kernel: where D 0 = λ/ρc 0 is the equilibrium thermal diffusivity of the glass-forming substance, τ 0 is a positive parameter, and 0 < ε 0 < 1. Usually, in glass-forming materials, ε 0 is in the range 0.2-0.3 [19,54] and sometimes even more [55].
Next, consider a one-dimensional example for a sample with a flat geometry. After that, the spherically symmetric problem in three dimensions can be reduced to the onedimensional problem mentioned above by using T(t, r) = U(t, r)/r.

Heat Equation with Dynamic Heat Capacity: Plane Geometry
As a basic example, consider a one-dimensional problem with homogeneous boundary and initial conditions. Consider an infinite plate of thickness d. Let the x-axis be directed along the normal to the plate surface. In many practically important cases, the heat flux Φ(t, x) can be considered as a multiplicatively separable function Φ(x)F(t), for example, in laser [56,57] and Joule heating [58], as well as local heating due to crystal nucleation [59][60][61]. Thus, we consider the heat flux in the form Φ(x)F(t) with continuous piecewise smooth F(t) and Φ(x) functions. Let the heat flux Φ(x)F(t) be distributed on the domain [0, d] and act during the time interval 0, τ p , i.e., F(t) = 0 for t ≤ 0 and τ p ≤ t. Suppose the dynamic heat capacity is described by Equation (4). Then, from Equation (5), we obtain: We focused on the dynamic behavior of the thermal response T(t, x) to the external heat flux Φ(x)F(t) in materials with a dynamic heat capacity. Suppose the temperature increases from the initial (thermostat) temperature. Thus, consider a homogeneous boundary value problem for the boundary conditions T(t, 0) = 0, T(t, d) = 0, and the zero initial condition T(t, x) = 0 for t ≤ 0. The solution T(t, x) can be represented as a series: where the functions ψ n (t) must satisfy Equation (8): LetT(t, x) and ψ n (t) denote the solutions of the conventional Equations (6) and (8) with ε 0 = 0, respectively. The solution ψ n (t) can be represented by Equation (9) [62]: For example, consider the solution to Equation where θ(t) is the Heaviside unit step function (with the condition θ(0) = 0). Denote this solution by χ n (t). Then, we have χ n (t) = Φ n ρc 0 τ n (1 − exp(−t/τ n )). Next, for a sinusoidal heat source F m (t) of duration τ p , where F m (t) = sin πmt/τ p for t ∈ 0, τ p and F m (t) = 0 Symmetry 2021, 13, 256 5 of 17 outside 0, τ p , we obtain the solution ψ n,m (t) of Equation (8) with ε 0 = 0, as seen in Equations (10) and (11): Thus, for the conventional Equation (6) with ε 0 = 0 and the sinusoidal heat source F m (t), the solution is:T To find the solutions T(t, x) and ψ n (t) of Equations (6) and (8) with nonzero positive ε 0 and τ 0 , we first solve an auxiliary problem with F(t) = θ(t), where θ(t) is the Heaviside unit step function, as seen in Equation (13). The solution to this auxiliary problem is denoted by ϕ n (t): Equation (13) can be solved using the Laplace transform method (see Appendix A for details). In fact, the solution is obtained similarly to the Volterra integral equations with a difference kernel [63]. Thus, we have: where −γ n and −µ n are the roots of the polynomial . The parameters γ n and µ n are real-valued and positive. Moreover, γ n − µ n > 0 for any positive τ n and τ 0 , and 0 < ε 0 < 1. It can also be shown that ϕ n (t) continuously tends to the solution χ n (t) = Φ n ρc 0 τ n (1 − exp(−t/τ n )) of Equation (8) with Let us consider the problem with nonzero ε 0 and τ 0 for the sinusoidal heat source F m (t). The solution to this problem can be obtained using the Duhamel integral [62,63] (see Appendix A for details). Thus, the solution to Equation (8) for the sinusoidal heat source F m (t) can be represented by Equation (16): From Equation (16) we obtain: Symmetry 2021, 13, 256 It can be shown that ψ n,m (t) continuously tends to the solution ψ n,m (t) as ε 0 → 0 or τ 0 → 0 . Thus, the solution to Equation (6) with nonzero ε 0 and τ 0 for the sinusoidal heat source F m (t) is: Finally, in the case of nonzero ε 0 and τ 0 , we obtain the solution to the problem for an arbitrary F(t). Indeed, if F(t) is a continuous piecewise smooth function on the interval 0, τ p , then the Fourier series in Equation (20) converges absolutely and uniformly to F(t) [36]. Thus, we represent F(t) on the interval 0, τ p using the Fourier series, as seen in Equation (20): where C m are the Fourier coefficients C m = 2 (6) is linear with respect to T(t, x), it follows that the solution T(t, x) of Equation (6) can be represented as a linear combination of solutions T m (t, x) for sinusoidal heat sources F m (t) with corresponding coefficients C m . Therefore, we obtain a solution to Equation (6) for an arbitrary continuous piecewise-smooth heat flux Φ(x)F(t) distributed on the domain [0, d] and acting during the time interval 0, τ p , as seen in Equation (21). The series shown in Equation (21) converges uniformly and absolutely [36]: Example 1. Let the heat flux Φ(x) be uniformly distributed on the domain d−x 0 2 , d+x 0 2 with density Φ 0 . Then: Denote byT p (t, x) the solution to the conventional Equation (6) with ε 0 = 0 for the rectangular heating pulse of duration τ p , i.e., where χ n (t) = Φ n ρc 0 τ n (1 − exp(−t/τ n )) and Φ n is determined using Equation (22). Note that τ n ∼ 1 n 2 . Thus, the series in Equation (23)  The solution of Equation (6) with ε 0 = 0 for the sinusoidal heat source F m (t) with m = 1 is equal to:T where ψ n,1 (t) is determined using Equations (10) and (11) for m = 1. Here again, Φ n is given by Equation (22).
Similarly, solutions to Equation (6) with nonzero ε 0 and τ 0 for the rectangular and sinusoidal heating sources are represented by Equations (25) and (26): where ϕ n (t) and ψ n,1 (t) are determined using Equation (14) and Equations (17) and (18), respectively, with Φ n being the same as in Equations (23) and (24). For example, let ε 0 = 1/3, τ 0 = 30 ns, and τ p = 10 ns. Let d = 100 nm, x 0 = 50 nm, and the thermal parameters of the substance are typical for those of glass-forming polymers. Thus, we take ρ = 1 g/cm 3 , ρc 0 = 2 × 10 6 J/m 3 K, λ = 0.3 W/mK, D 0 = 1.5 × 10 −7 m 2 /s, and Φ 0 = ρh 0 , where h 0 = 200 J/g is the heat release during crystallization. We focused on the nanometer and nanosecond scales since they are close to real processes during subcritical crystal nucleation [59][60][61]. First, we verified that Equation (21) gave the correct solution T(t, x) for a rectangular pulsed heat source. Indeed, the solutionsT p (t, x) and T p (t, x) represented by Equations (23) and (25) coincided with the resultsT(t, x) and T(t, x) calculated using the Fourier series (see Equation (21)) for m up to 45 (see Figure 1b). In addition, Figure 1b shows that the effect of the temporal dispersion of the dynamic heat capacity on the temperature was significant. In addition, Figure 1b shows that the effect of the temporal dispersion of the dynamic heat capacity on the temperature was significant. As an example of an arbitrary function ( ), consider the continuous piecewise smooth function ( ) shown in Figure 2. Let ( ) denote the Fourier series approximating ( ), as seen in Equation (27): For example, let us take = 30 and = 30. Figure 2a shows that ( ) was well approximated by ( ) at = 30. The solution ( , ) for the heating source with ( ) is represented by Equation (28), where , ( ) are determined by Equations (17) and (18) with corresponding coefficients (see Equation (27)): ( , ) = , ( ) ( ) . As an example of an arbitrary function F(t), consider the continuous piecewise smooth function F A (t) shown in Figure 2. Let F A (t) denote the Fourier series approximating F A (t), as seen in Equation (27): For example, let us take N = 30 and M = 30. Figure 2a shows that F A (t) was well approximated by F A (t) at M = 30. The solution T A (t, x) for the heating source with F A (t) is represented by Equation (28), where ψ n,m (t) are determined by Equations (17) and (18) with corresponding coefficients C A m (see Equation (27)):  The solutions T A (t, x) and T 1 (t, x) can be compared withT A (t, x) andT 1 (t, x), respectively. Let δT(t, x) = (T(t, x) −T(t, x))/T(t, x) denote the relative contribution to the solution associated with the temporal dispersion of the dynamic heat capacity. This contribution reached about 50% (see Figure 2b). Thus, the effect of the temporal dispersion of the dynamic heat capacity on T(t, r) was significant, especially at the beginning of the heating process. The position of the peak of the time dependenceT 1 t, d 2 was shifted relative to the peak of the heating pulse F 1 (t), as well asT A t, d 2 relative to F A (t) (see Figure 2b). Indeed, the thermal response usually lagged behind the heat source. Thus, the peak of the dependencesT 1 t, d 2 andT A t, d 2 appeared at around 8 ns and 7 ns, respectively. However, the maxima of F 1 (t) and F A (t) were at 5 ns and about 4 ns, respectively. Interestingly, this shift decreased due to the dynamic heat capacity. Thus, the peaks of the dependences T 1 t, d 2 and T A t, d 2 appeared around 7.5 ns and 6.5 ns, respectively. In fact, the dynamic heat capacity was less than the equilibrium heat capacity, especially at the beginning of the heating process.
The cylindrically symmetric problem can be solved in the same way as in [29]. Below, we consider only an example of a spherically symmetric problem.

Heat Equation with Dynamic Heat Capacity: Spherical Geometry
Consider a spherically symmetric problem with a spherical heat source Φ(r)F(t) of radius r 0 that is concentrically located in a spherical sample of radius R. The boundary condition is T(t, R) = 0 and the initial condition is T(t, r) = 0 for t ≤ 0. From Equation (5), we have: By replacing rT(t, r) with U(t, r) in Equation (29), we obtain a one-dimensional problem, as seen in Equation (30); this problem is similar to that discussed earlier. The boundary and initial conditions are U(t, 0) = 0, U(t, R) = 0, and U(t, r) = 0 at t ≤ 0. where the prime means the derivative with respect to the time variable. Thus, the solution T(t, x) is: where the functions ψ n (t) must satisfy Equation (32): where Φ n = 2 R R 0 rΦ(r)sin(πnr/R)dr and τ −1 n = D 0 (πn/R) 2 . The functions ψ n (t) for different F(t) can be obtained in the same way as before (see Example 1).

Example 2.
Let the heat flux Φ(r) be uniformly distributed on the domain [0, r 0 ] with density Φ 0 . Then: For heating pulses F p (t) and F A (t), similar to those in Figures 1 and 2, consider the solutionsT p (t, r) and T p (t, r), as well asT A (t, r) and T A (t, r). In addition, consider the solutionsT 2 (t, r) and T 2 (t, r) for the sinusoidal heat source F m (t) with m = 2. Let R = 300 nm, r 0 = 30 nm, τ p = 2 ns, ε 0 = 1/3, and τ 0 = 5 ns, along with the same thermal parameters as in the above example. The series in Equation (31) converges as fast as the series S N = ∑ N n=1 1/n 2 , which converges to π 2 /6 as N → ∞ [64]. The remainder π 2 /6 − S N is less than 1% of π 2 /6 for N = 65. Thus, to obtain an accuracy better than 1%, it is sufficient to take the series in Equation (31) up to N = 100. All the examples below were calculated with N = 100 and, as before, M = 30. Since R r 0 , it does not matter how large the parameter R is, as long as we consider a sufficiently short time t. This was verified by direct calculation ofT A and T A at R = 300 nm and 1000 nm (see Figure 3). Moreover, the solutions for different heating pulses practically did not change with the distance at r > 80 nm (see Figure 3). As seen in Figure 3, the effect of the temporal dispersion of the dynamic heat capacity was significant. Interestingly, the effect of the temporal dispersion of the dynamic heat capacity could even lead to a change in the sign of the solution T 2 (t, r). Indeed, T 2 (t, r) had a sign that was opposite toT 2 (t, r) at r < 21 nm and t 0 = 1.5 ns (see Figure 3b). A similar effect can occur when the temperature changes during the nucleation and decay of subcritical crystals.
The effect of the temporal dispersion of the dynamic heat capacity was most significant in the case of fast and local heating; this effect was most pronounced at the beginning of the heating process. This effect increased with increasing τ 0 and reached saturation at τ 0 in the order of r 0 2 /D 0 , that is, at about 10 ns at r 0 = 30 nm and D 0 in the order of 10 −7 m 2 /s (see Figure 4). Note, in the case of glass-forming substances, the relaxation times τ 0 have a broad distribution. Next, we considered the influence of this distribution on the solution T(t, r) for the spherically symmetric problem. capacity was significant. Interestingly, the effect of the temporal dispersion of the dynamic heat capacity could even lead to a change in the sign of the solution ( , ). Indeed, ( , ) had a sign that was opposite to ( , ) at < 21 nm and = 1.5 ns (see Figure 3b).
A similar effect can occur when the temperature changes during the nucleation and decay of subcritical crystals. The effect of the temporal dispersion of the dynamic heat capacity was most significant in the case of fast and local heating; this effect was most pronounced at the beginning of the heating process. This effect increased with increasing and reached saturation at in the order of / , that is, at about 10 ns at = 30 nm and in the order of 10 −7 m 2 /s (see Figure 4). Note, in the case of glass-forming substances, the relaxation times have a broad distribution. Next, we considered the influence of this distribution on the solution ( , ) for the spherically symmetric problem.  The effect of the temporal dispersion of the dynamic heat capacity was mo cant in the case of fast and local heating; this effect was most pronounced at the of the heating process. This effect increased with increasing and reached sat in the order of / , that is, at about 10 ns at = 30 nm and in the ord m 2 /s (see Figure 4). Note, in the case of glass-forming substances, the relaxatio have a broad distribution. Next, we considered the influence of this distributi solution ( , ) for the spherically symmetric problem.

Dependence of the Solution T(t, r) on the Distribution of Relaxation Times τ 0
As before, consider a spherically symmetric problem with a spherical heat source Φ(r)F(t) of radius r 0 that is concentrically positioned in a spherical sample of radius R. The boundary and initial conditions are the same as in Section 5. Suppose the dynamic heat capacity is represented by Equation (3) or by the same equation in a finite interval [τ min , τ max ]. Then, using the distribution function H(τ 0 ), we can obtain the desired solution as a linear combination of solutions for different τ 0 . The distribution function H(τ 0 ) can be found using broadband heat capacity spectroscopy [19].

Example 3.
Consider the Kohlrausch relaxation law (stretched exponential), which is often used for glass-forming substances. Then, the dynamic heat capacity can be represented using Equation (34): For example, let β = 0.5. Then, the normalized distribution function is 4πτ K τ 0 [50,51], where the Kohlrausch relaxation time τ K determines the distribution width, and the average relaxation time is τ 0 AV = 2τ K [64]. In fact, the shape of the distribution function H(τ 0 ) is not very significant since the effect of the temporal dispersion of the dynamic heat capacity is saturated with increasing τ 0 . For example, consider the uniform distribution H u (τ 0 ) = 1 4τ K on the interval [0, 4τ K ] and H u (τ 0 ) = 0 outside of this interval. In this case, the average relaxation time is also 2τ K . Comparing the results calculated using the Kohlrausch distribution H(τ 0 ) = exp(−τ 0 /4τ K ) √ 4πτ K τ 0 and the uniform distribution H u (τ 0 ), it can be seen that the results are very close to each other (see Figures 5 and 6). The value of the average relaxation time is a more significant parameter.
The dynamic heat capacity c dyn (t) can be set on a finite interval [τ min , τ max ] if this interval is wide enough. It is sufficient to take [τ min , τ max ] = [0.001τ K , 15τ K ] for an accuracy of about 1%. Thus: The distribution function in Equation (35) depends on the Kohlrausch relaxation time τ K ; this relaxation time can be found using broadband heat capacity spectroscopy [19], as was done in [28]. Heat capacity spectroscopy allows one to obtain the frequency dependence of the dynamic heat capacity. Let ω max denote the frequency of the maximum of the imaginary part of the dynamic heat capacity. In fact, τ K = 0.74/ω max [28]. ω max can be obtained for different temperatures T from the Vogel-Fulcher-Tammann-Hesse relation log(ω max ) = A − B/(T − T 0 ), where A, B, and T 0 are measured in the experiment [19]. Thus, one can obtain the Kohlrausch relaxation time τ K (T) for different temperatures. Consider the solution T(t, r, τ K (T in )) for τ K (T in ), where T in is the initial temperature at t ≤ 0. Compare this solution with T(t, r, τ K (T max )), where T max is the maximum temperature reached by the sample during the pulsed heating. Since τ K (T) changes with heating between τ K (T in ) and τ K (T max ), it follows that the correct solution T(t, r) is between the solutions obtained for τ K (T in ) and τ K (T max ), i.e., T(t, r, τ K (T max )) < T(t, r) < T(t, r, τ K (T in )).
For example, let T in = 500 K. Consider the solution T 1 (t, r) for the sinusoidal heat source F m (t) with m = 1, where F m (t) = sin πmt/τ p for t ∈ 0, τ p and F m (t) = 0 outside 0, τ p . Let us take the same thermal and size parameters as in the above example. Thus, we obtained τ K (T in ) = 13.1 ns for T in = 500 K and τ K (T max ) = 3.2 ns for T max = 553 K, using the parameters A = 10.2, B = 388 K, and T 0 = 341.5 K for polystyrene (PS) [19]. The difference between the solutions obtained for τ K (T max ) and τ K (T in ) was less than 5% of the temperature change near the maximum of the curves shown in Figure 5. The influence of the temporal dispersion of the dynamic heat capacity was much greater than this difference (see Figure 5). Moreover, the effect of changing τ K (T) upon heating was insignificant at lower T in since τ K (T) increased with decreasing T and the influence of the temporal dispersion of the dynamic heat capacity saturated with increasing τ K (see Figure 6). Therefore, we obtained τ K (T in ) = 200 µs for T in = 400 K and τ K (T max ) = 120 ns for T max = 455 K using the parameters A, B, and T 0 for PS [19]. Thus, we obtained the solutions T(t, r, τ K (T in )) and T(t, r, τ K (T max )), which were practically the same (see Figure 6). In fact, the difference between T(t, r, τ K (T in )) and T(t, r, τ K (T max )) was negligible in a wide temperature range from the glass transition temperature to 400 K for PS. Similar calculations for poly(methyl methacrylate) (PMMA) provided an insignificant difference between the solutions T(t, r, τ K (T in )) and T(t, r, τ K (T max )) calculated in a wide temperature range from the glass transition temperature even up to 700 K (see Figure 7). For example, we obtained τ K (T in ) = 127 ns for T in = 700 K and τ K (T max ) = 107 ns for T max = 755 K using the parameters A = 7.3, B = 185 K, and T 0 = 354.3 K for PMMA [19]. Thus, the calculation of the influence of the temporal dispersion of the dynamic heat capacity for constant τ K (T in ) was quite accurate, even if the temperature change was about 50 K (see Figures 5-7).
Summing up, we concluded that the calculations can be carried out for a constant Kohlrausch relaxation time, determined at T in , and for sufficiently large (about 50 K) temperature changes during heating. This holds for a wide temperature range for T in . Summing up, we concluded that the calculations can be carried out f Kohlrausch relaxation time, determined at , and for sufficiently large (abo perature changes during heating. This holds for a wide temperature range fo

Conclusions
An integro-differential equation describes the non-equilibrium therma glass-forming substances with a dynamic (time-dependent) heat capacity to perturbations. We found that the corresponding heat transfer problem cou analytically for a heat source with an arbitrary time dependence and differen The solutions provide analytical expressions for fast thermal processes on and longer timescales. It was shown that the effect of the time dispersion of heat capacity on the thermal response ( , ) upon local heating was signifi forming materials. This effect was enhanced when the relaxation time of the d capacity increased. However, this effect reached saturation at tens of nano cause of this saturation, in many practical cases, the effect of the time disp dynamic heat capacity on the thermal response ( , ) can be calculated usin tion time distribution fixed at a constant initial temperature. This effect can highly localized and short processes, such as crystal nucleation. The method calculations described in this work can be applied to an arbitrary time dep local thermal perturbation, for example, the formation of a crystal nucleus i ing materials. The results obtained can be useful for applications and a better ing of the thermal properties of glass-forming substances, polymers, and nan   (a) Time dependences ofT 1 (t, 0), T 1 (t, 0, τ K (T in )), and T 1 (t, 0, τ K (T max )) (shown using squares, circles, and triangles, respectively) and (b) the corresponding spatial temperature distributions at t 0 = 1.5 ns for PMMA at T in = 700 K and a sinusoidal heat source F 1 (t). PMMA: poly(methyl methacrylate).

Conclusions
An integro-differential equation describes the non-equilibrium thermal response of glass-forming substances with a dynamic (time-dependent) heat capacity to fast thermal perturbations. We found that the corresponding heat transfer problem could be solved analytically for a heat source with an arbitrary time dependence and different geometries. The solutions provide analytical expressions for fast thermal processes on nanosecond and longer timescales. It was shown that the effect of the time dispersion of the dynamic heat capacity on the thermal response T(t, r) upon local heating was significant in glassforming materials. This effect was enhanced when the relaxation time of the dynamic heat capacity increased. However, this effect reached saturation at tens of nanoseconds. Because of this saturation, in many practical cases, the effect of the time dispersion of the dynamic heat capacity on the thermal response T(t, r) can be calculated using the relaxation time distribution fixed at a constant initial temperature. This effect can be crucial for highly localized and short processes, such as crystal nucleation. The method of analytical calculations described in this work can be applied to an arbitrary time dependence of a local thermal perturbation, for example, the formation of a crystal nucleus in glass-forming materials. The results obtained can be useful for applications and a better understanding of the thermal properties of glass-forming substances, polymers, and nanocomposites.