Local heating variations and transient effects in the coupling of thermal radiation and non-Fourier heat transport

Thermalization is the process by which a system reaches a state of thermal equilibrium with its surroundings. Conduction and radiation are two fundamental modes of heat transfer that play a crucial role in a wide range of applications ranging from thermal therapies Andreozzi et al. (2019) to heat management at the nanoscale Stewart et al. (2022); Smoyer and Norris (2019). In some cases, these two modes can interact, leading to a phenomenon known as conduction-radiation coupling. In the near-field, thermalization can occur on a much faster timescale than in the far-field due to the proximity of the two bodies Pendry (1999). Radiative heat transfer (RHT) between two bodies separated by a gap at different temperatures has become a topic of interest in energy management at the micro and nanoscale Lucchesi et al. (2021). For large separations, this phenomenon is well described by the Planck theory of black body radiation, which only considers propagating electromagnetic waves Cuevas and García-Vidal (2018). However, at subwavelength separations between the bodies, known as the near-field regime, strongly localized evanescent waves become the main contributors to RHT Pascale et al. (2023). Unlike the far-field case, the heat flux depends on the separation between the bodies. Close to room temperature of $300\;$K, thermal radiation in the far-field is limited to $\approx 6\;\text{Wm}^{-2}\text{K}^{-1}$ by Planck’s blackbody radiation law. When tens of nanometers separate surfaces, radiative thermal conductance can increase several orders of magnitude above the black body prediction Polder and Van Hove (1971); Kittel et al. (2005); De Wilde et al. (2006). In the original theory describing radiative heat exchanges between two closely spaced bodies introduced by Polder and van Hove Polder and Van Hove (1971), there is no coupling between the heat carriers inside the materials and the electromagnetic energy transferred across the vacuum gap.

This interplay between heat transport and radiative heat transfer in thermal relaxation has been explored recently Reina et al. (2020, 2021); Messina et al. (2016). In these works, the Boltzmann transport equation achieves a general coupling of heat radiation and conduction for bodies of arbitrary size. This leads to study systems exhibiting diffusive or semi-ballistic behavior Hoogeboom-Pot et al. (2015). These authors focus mainly on planar structures comprised of polar materials. It was concluded that the radiative heat flux exchanged between parallel slabs at nanometric distances is reduced due to the change in the temperature profiles within each body.

This paper addresses the case when heat is exchanged between non-Fourier materials. Non-Fourier heat transport refers to heat transfer in which the temperature gradient is not directly proportional to the heat flux. The heat flux can exhibit non-local, non-linear, or memory-dependent behavior, and the temperature field can propagate at a finite speed de la Rosa and Esquivel-Sirvent (2022).

To mend this unphysical consequence, Cattaneo and Vernotte Cattaneo (1958); Vernotte (1958); Spigler (2020); Ostoja-Starzewski (2009), independently proposed inserting a delay time $\tau$ between applying the temperature gradient and the heat flux Maillet (2019). One example of a material that exhibits non-Fourier behavior is a thin film of metal, such as gold or silver, which can have a much higher thermal conductivity than the surrounding air or substrate Sobolev and Dai (2022). When a laser pulse is focused on such a material Jaunich et al. (2008), the temperature of the metal can rapidly increase, leading to non-equilibrium behavior and the onset of thermalization. Further applications of non-Fourier heat conduction can be found in describing the temperature profile dynamics in biological tissue Camacho De La Rosa et al. (2021); Gupta and Srivastava (2019).

In this work, we will discuss the coupling between radiation and conduction between two closely separated solids, provided the conduction properties do not follow Fourier’s law of heat conduction. For the separation range considered, radiative energy exchange competes with respect to the conductive heat transport within the bodies Messina et al. (2016). We do not consider the transition from radiation to conduction at sub-nanometer gaps Chiloyan et al. (2015).

We consider two planar bodies of width $L$ separated by a vacuum gap $d$ at initial temperatures $T_{01}$ and $T_{02}$ respectively, as shown in Fig.1. Body 2 is also connected to a thermal bath at temperature $T_{02}$. The bodies exchange energy only via RHT. Energy transferred via electromagnetic waves is then absorbed into the body and distributed according to a hyperbolic thermal transport equation, unlike Fourier materials that follow a parabolic-type equation. Within a slab, changes of the internal energy density $u(\vec{x},t)$ are due to the conductive and radiative energy fluxes, $\Phi_{c}$ and $\Phi_{r}$, thus

A constitutive equation can describe the conductive energy flux in both Fourier and non-Fourier materials. When the medium follows the Fourier equation, we will have $-\kappa\nabla T=\vec{q}$, where $\kappa$ is the thermal conductivity [W m^-1 K^-1], we assume it is a scalar valid for an isotropic medium. On the other hand, when the medium is described by the Cattaneo-Vernotte equation, we have

where $\tau$ is a delay time. The second term on the right-hand side can be seen as a first-order approximation to the description of the time lag between the heat flux $\vec{q}$ and changes in the temperature gradient $\nabla T$. For the limiting case of instantaneous response $\tau\rightarrow 0$, the diffusive Fourier’s law is recovered.

To obtain an equation for the temperature, we first notice that slab 2 is in contact with a thermal bath at constant temperature $T_{20}$, it will be useful to define the auxiliary function

where $k=1$, 2 indicates the slab. The function $F_{k}$ describes deviations from the equilibrium temperature of the system and is therefore equal to zero when the slabs reach thermal equilibrium. To obtain a partial differential equation (PDE) for the temperature fields difference $F_{k}(\vec{x},t)$, we can use the energy conservation equation Eq.(1), the CV heat conduction equation Eq.(2), and the constitutive equation for the internal energy $u=c_{v}T$, where $c_{v}$ is the specific heat per unit volume [J kg^-1 K^-1]. We get

where $\alpha=\kappa/c_{v}$ is the thermal diffusivity [m² s^-1], $v^{2}=\alpha/\tau$ is known as the heat propagation speed, and where $\Phi_{r}$ has taken the role of an external energy source. Notice that while Eq.(4) is the general 3-d expression, due to the symmetry of the problem, we only need to consider the spatial variable $z$, the axis perpendicular to the slabs. Eq.(4) describes the time-dependent temperature profile within each slab. With no external energy source, Eq.(4) describes a damped wave. This starkly contrasts to the Fourier case, where the temperature profile decays exponentially with time without an external energy source.

To obtain an expression of the radiative energy flux $\Phi_{r}$, we further assume that the conductive channel is the dominant mode of energy transfer within each slab and therefore neglect RHT between portions of the same slab. Furthermore, we use the heat-sink approximation, valid for polar dielectrics, which considers RHT to occur only at the slab interface Reina et al. (2021). Application of this approximation allows us to neglect the source term $\Phi_{r}$ within each slab, and the PDE becomes homogeneous with initial and boundary conditions describing the RHT contribution of energy. We will, therefore, only require the description of $\left.\Phi_{r}\right|_{z=0}$ and $\left.\Phi_{r}\right|_{z=d}$ which will be imposed in the initial and boundary conditions.

The radiative energy flux $\Phi_{r}$ evaluated at the interfaces can be related to the radiative thermal conductance defined as

where $\Phi^{1\rightarrow 2}$ is the net flux per unit area received by body 2 from 1. Explicitly, the radiative thermal conductance can be written as

As previously mentioned, instead of appearing explicitly in Eq.(4), the initial conditions hold the contribution of $\Phi_{r}$ at the interfaces. Since we have a second-order equation for $t$ and $z$, we require 2 initial and 2 boundary conditions. These can be choose by first observing that at an initial time, $t=0$, each body has a constant temperature, $T_{01}$ and $T_{02}$, respectively so that $F_{1}(z,t=0)=\Delta T$ and $F_{2}(z,0)=0$. Notice that since the bodies are at constant temperatures at the initial time $t=0$, there is no heat flux in each body $\Phi_{c}=0$, so at the initial time $t=0$, the energy conservation equation gives $c_{v}\partial F_{1}(z,t)/\partial t\mid_{(t=0)}=\Phi_{r}=G_{R}\Delta T\delta(z)$ and $c_{v}\partial F_{2}(z,t)/\partial t\mid_{(t=0)}=\Phi_{r}=G_{R}\Delta T\delta(z-d)$. For body 1 we further require that the heat flux be zero at the position farthest from the gap and for slab 2 we impose $F_{2}(z,t)=0$, since it is in contact with a thermal bath. Finally, at the slab boundaries closest to the gap we require the conductive flux to be equal to the radiative flux. The initial and boundary conditions can be summarized as

Notice that compared to the Fourier solution, in the CV case there is an additional initial condition for each slab, corresponding to conditions 3) and 4). The imposition of the initial and boundary conditions leads to a summation and is shown in the supplemental material. Solutions of temperature profiles can be written as

which we then use to calculate these terms

Notice that in our system, the material of slabs 1 and 2, are the same so that the subindex $i$ in Eq.(10) can be dropped. The thermalization dynamics described by the Fourier equation can be found following an analogous procedure and is written here for the sake of completeness, a derivation of this expression can be found in the literature Reina et al. (2021),

We begin by considering a system composed of two slabs of $L=100$ $\mu$m composed of SiC separated by a vacuum gap of $d=1$ nm. As a first task, we show the deviations between the thermalization dynamics of the slabs when described by Fourier’s law and the hyperbolic heat equation. A time delay $\tau$ must be chosen for the SiC slabs to do this. The time constant for only a few materials has been reported K. Mitra and Moallemi (1995); Madhukar et al. (2019), and a general theory to calculate the CV relaxation time constant for a given material is still lacking. Past work which describes hyperbolic thermal transport in semiconductor materials proposed using a time constant extrapolated from the phonon relaxation time. Zul Karnain et al. (2020). While this approximation is reasonable, it is worth considering that the time constant $\tau$ in the CV describes the system heat flux time response to a thermal source, which is not necessarily the same as the phonon lifetime of the system. With this in mind, as a first approximation, we take the delay time in the upper limit but within the range of known phonon lifetimes for SiC Reina et al. (2020).

In Fig.2, we show the time-dependence of the spatial average temperature $\left<T_{1}(z,t)\right>=\frac{1}{L}\int_{-L}^{0}T_{1}(z,t)dz$ of slab 1 when described by the Fourier and CV equation with a delay time $\tau=10^{-8}$ s. Fig.2(a) shows a nearly equal dynamic between the Fourier (dashed line) and hyperbolic (solid line) equations for such a small value of $\tau$, with the slabs reaching thermal equilibrium within the first few milliseconds in correspondence with past work Reina et al. (2021). In Fig.2(b), we plot the temperature difference between the Fourier and CV case as a function of space and time within body 1, $\mathcal{T}_{1}(z,t)=T_{1}(z,t)-T_{1}^{F}(z,t)$. We show a zoom-in for $t<0.003$ ms $\approx 100\;\tau$ in the upper right side, where we expect the largest wave-like behavior. It can be seen that the temperature difference is small, with the difference tending to zero with increasing time, indicating the slab has reached thermal equilibrium. The inset of Fig.2(b) shows a zoom of the first half millisecond, where we observe that $\mathcal{T}_{1}$ exhibits an oscillatory behavior that eventually dies out, indicative of a transient behavior in the thermalization process. As expected, the difference between diffusive and hyperbolic behavior for small delay time values is negligible. However, as the delay time grows, the difference between both conduction models is more noticeable. To further understand the origin of this dynamics difference, we point out that the spatial part of the solution is the same for Fourier and hyperbolic equations (Eq.(9) and Eq.(13)). The difference in the dynamics stems from the temporal part of the solution.

The roots that define the exponential behavior of the CV equation temporal solution are given by $r^{\pm}_{n}$, which are written in Eq.(10). Unlike the Fourier case, with just one real root, when using the CV equation these roots can have an imaginary part depending on the system parameters. The dependence of the parameter $r^{\pm}_{n}$ as a function of the CV delay time $\tau$ is particularly interesting. The inverse of the real part of coefficient $r^{\pm}_{n}$ defines a sort of lifetime of the mode $t_{n}=1/r_{n}$, where $r_{n}=\text{min}(Re[r^{+}_{n}],Re[r^{-}_{n}])$, which can be compared to the analogous quantity obtained for the Fourier case, written as $1/r^{F}_{n}=L^{2}/(x_{n}^{2}\alpha)$. In Fig.3 we show the lifetimes of the first three modes $n=$ 1 (solid), 2 (dashed), and 3 (dotted) as a function of the CV response time $\tau$.

The series of the temperature solution (Eq.(9)) converges because as $n$ grows, the associated mode decreases its contribution. In Fig.3(a) for each $\tau$ we observe that $t_{1}(\tau)\geq t_{2}(\tau)\geq t_{3}(\tau)$. As we increase $\tau$, the imaginary part of $r_{n}$ takes relevance, indicating the threshold where is relevant the wave-like behavior.

In contrast, the number of modes required to accurately describe the temperature dynamics near the temporal origin will depend on the value of $\tau$. As expected, the temporal window where a non-Fourier description is required depends on the value of $\tau$. For a value of $\tau=10^{-8}$ s, modes with $n\geq 2$ die off after less than a tenth of millisecond. In Fig.3, we can see that, unlike the Fourier case, the temporal coefficient $r_{n}$ can have an imaginary part for a given $\tau$. It is interesting to note that for the case of $\tau=10^{-8}$s, there is no imaginary part to the first mode. An imaginary part of $r_{n}$, and therefore the oscillatory nature of the solution, will be presented only for large values of $n$ with lifetimes of the order of $\tau$.

To further illustrate the wave-like behavior of the solution, we point out that the physical temperature is obtained by taking the real part of Eq.(9). With this in mind, we can recast each term of the sum Eq.(9) in the form $T(z,t)=T_{02}+\sum_{n=1}^{\infty}\text{Re}[T_{1n}+T_{2n}]$, where each term $T_{in}$ satisfies that

where we write $r_{in}=r^{\prime}_{in}+ir^{\prime\prime}_{in}$, and the coefficient $A_{n}=A_{n}^{\prime}+iA_{n}^{\prime\prime}$, it saves the information of $a_{n}$ and $b_{n}$. This shows that the sum is just a superposition of damped waves. Furthermore, it allows us to define the velocity of each mode as $v_{n}=\frac{r^{\prime\prime}_{i,n}L}{x_{n}}$. The velocity cannot be defined when no imaginary part of $r_{i}$ exists. In Fig.3, we show the velocity of the first three modes as a function of the parameter $\tau$. We can see that the mode velocity increases as $n$ increases, and even the slowest mode ($n=1$) travels 100 $\mu$m in about 0.2 milliseconds.

While deviations from Fourier transport can be small when $\tau$ is of the order of the phonon lifetime, it is interesting to consider larger response times. Physically this can be obtained by inserting defects into the material so that the thermal response time no longer matches the material phonon lifetime. To this end, we consider larger values of $\tau$, which may lead to oscillating modes with larger lifetimes $t_{n}$. In Fig.4, we show the difference in thermalization dynamics for two different values of $\tau$. For these cases, the temperature differences are much larger than the corresponding small values of $\tau$, being of the order of tens of degrees. This larger deviation can be traced back to the coefficients $a_{n}(\tau)$ and $b_{n}(\tau)$ of Eq.(9), which depend on $\tau$ through parameter $r^{\pm}_{n}$. As in the previous case, the largest deviations are found after that onset the thermalization. The alternating differences in temperature shown in Fig.4 indicate the wave-like behavior of the CV solution.

We have established the influence of the RHT in the boundaries $z=0$ and $z=d$, acting as energy sources. In the diffusive case, this energy perturbation is felt instantaneously in the whole slab, while for the hyperbolic model, the introduction of $\tau$ implies a delay in the transport of energy; thus the farthest position felt later the thermal perturbation than the first point, giving rise to a wave-like behavior. In Fig.4, we can see a minimum in $\mathcal{T}_{1}$ that reflects back and forth between the two interfaces of slab 1, it means that Cattaneo’s model predicts an oscillation around Fourier’s law for the first moments, indicative of a transitory effect in the thermalization. As time increases, the hyperbolic solution predicts a lower temperature than the diffusive’s case, until both temperatures coincide and $\mathcal{T}_{1}\rightarrow 0$. Both models reach the thermal equilibrium at similar times but with different rates of change. For the case of $\tau=10^{-4}$ s in Fig.4(b) a minimum in $\mathcal{T}_{1}$ is appreciable for nearly $\approx$ 1 ms = $10\;\tau$ and travels at a speed roughly equal to $v_{1}\approx 1000\;\mu$m/s, suggesting that this minimum is due mainly to the interference of mode n = 1 of the thermal waves (see Fig.3(b)).

In Fig.5 we show the temperature difference evaluated at the interface of body 1 with the gap. This is equivalent to a vertical cut of the temperature maps of Fig.4 in $z=0$. In Fig.5(a) we can see that the oscillation periods are not constant during the thermalization, this is consistent with the $\mathcal{T}_{1}=0$ lines (isoclines) shown in Fig.4. Besides, this agrees with the fact that at small times $t$, modes with larger $n$ and smaller periods contribute to the temperature profile. After half a millisecond the period becomes more regular since at this time, fewer modes have a perceptible contribution to the thermalization dynamic. In Fig.5(b), we show the averaged temperature within body 1, we observe as the diffusive and hyperbolic models reach the thermal equilibrium at similar times, but the rate of change of $\left<T_{1}(z,t)\right>$ is faster for $\tau=10^{-4}\;$s. It means that modifying the delay time affects the rate of change of the average temperature, decreasing faster for larger response times. While inside the body, the difference in the temperature profile can be up to around 40 degrees, the average temperature inside body 1 of Cattaneo’s model for $\tau=10^{-4}\;$s deviates only about 10 degrees with respect to the Fourier case.

We observe that the solutions of the diffusive and hyperbolic models reach the thermal equilibrium in similar times, the dissimilarity between both models depends on the delay time. The two bodies exchange energy via near-field thermal radiation, and when coupled to heat conduction, it appears as a modification in the dynamics of the temperature profile. On the other hand, the diffusive model predicts a slower transitory process than the hyperbolic case for large values of $\tau$, however as time increases the difference between both models is negligible.

In this work, we have studied the thermalization dynamics of two bodies that exchange energy via near-field thermal radiation using a hyperbolic heat conduction equation. For the case of SiC, we have concluded that if the thermal time response is of the order of the systems phonon lifetime, the thermalization dynamics described with the CV equation remains the same, to the order of a millidegree, when compared to the Fourier description. As the response time increases, the thermalization dynamics change, showing that the rate or speed of thermalization is slowed. Furthermore, at the onset of the thermalization process, the variation in temperature at the surface of the hotter body exhibits oscillations typical of transient behavior. These oscillations decrease with time, reaching a steady-state thermalization rate. Large deviations in the temperature profile of tens of degrees can be expected when material response time is a tenth of a millisecond. This could have important implications in applications such as scanning probe lithography, where precise control of the temperature profile is required, particularly when considering temperature profiles in highly inhomogeneous material. Future work could also consider the case of layered systems of dielectrics and biological-like tissue, in particular, to have thermal shields with the potential of drastically slowing down the thermalization process.

Finally, thermalization in the near-field is an important process in many physical systems and can be accurately described using the hyperbolic heat conduction equation. This equation considers the non-Fourier behavior of materials and allows for the rapid propagation of heat waves at velocities that exceed the speed of sound. Understanding thermalization in the near field is essential for designing microelectronic devices, optimizing laser-based manufacturing processes, and developing new materials with tailored thermal properties.