# Thermal conductivities of amorphous SiO2 films by approach-to-equilibrium molecular dynamics

Louis Petit, Cheick O Diarra, Achille Lambrecht, Ibrahim Bel-Hadj, Mauro

Boero, Évelyne Martin

<figure><img src="image_1.png" /></figure>

► To cite this version:

Louis Petit, Cheick O Diarra, Achille Lambrecht, Ibrahim Bel-Hadj, Mauro Boero, et al.. Thermal conductivities of amorphous SiO2 films by approach-to-equilibrium molecular dynamics. Solid State Sciences, 2026, 175, 10.1016/j.solidstatesciences.2026.108262. hal-05514396

HAL Id: hal-05514396

https://hal.science/hal-05514396v1

Submitted on 17 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Thermal conductivities of amorphous SiO₂ films by approach-to-equilibrium molecular dynamics☆

Louis Petit<sup>ID</sup>, Cheick O. Diarra<sup>ID</sup>, Achille Lambrecht<sup>ID</sup>, Ibrahim Bel-Hadj<sup>ID</sup>, Mauro Boero<sup>ID</sup>, Evelyne Martin<sup>ID*</sup>

University of Strasbourg, CNRS, ICube Laboratory UMR 7357, 300 Boulevard Sébastien Brant, Illkirch, F-67412, France

## ARTICLE INFO

**Keywords:**

Thermal conductivity

Film

Amorphous SiO₂

## ABSTRACT

The thermal conductivity of stand-alone amorphous SiO₂ films is calculated both in and cross plane. In both cases, the approach-to-equilibrium molecular dynamics (AEMD) method is employed, showing how the method can be extended to calculate the cross-plane thermal conductivity. The study is conducted within the first-principles molecular dynamics framework. We ascertain that the thermal conductivity in the plane remains unaffected by the presence of free surfaces. Remarkably, even for a thickness of 2 nm, the thermal conductivity is found to be equal to (1.3 ± 0.1) W K⁻¹ m⁻¹, which is in accordance with the measured macroscopic value. With regard to cross-plane thermal conductivity, the value for a film with thickness *d* is equivalent to that obtained by AEMD in the case of a bulk with a periodicity of 2*d*. The cross-plane thermal conductivity is shown to decrease for thicknesses below 3 nm.

## 1. Introduction

A reduction in thermal conductivity is generally observed when the scale of materials is reduced from macroscopic dimensions to smaller sizes, as occurs in thin films. This effect has the consequence of impacting the dissipation of heat in nanotechnological devices. A premier material in this type of devices is SiO₂ in its vitreous state. Amorphous SiO₂ (aSiO₂) films are formed either intentionally, for example in SOI (silicon on insulator) technologies, or not, as is the case of native layers. The bulk thermal conductivity of bulk aSiO₂ is low. Yet, with the progressive size reduction of devices, it is crucial to ascertain if there exists a layer thickness beneath which the thermal conductivity might be further reduced. The available experimental data [1] fail to demonstrate a reduction for thicknesses down to 60 nm, other than those associated with structural changes, such as an increased porosity. On the theoretical forefront, it was found [2] that the cross-plane thermal conductivity of an amorphous silica film approximately 4 nm thick was not reduced. However, the interatomic potential employed in these simulations is restricted in its capacity to delineate the structure of the material, particularly the free surfaces of films.

In our previous work [3], we studied size effects in the thermal conductivity of aSiO₂ by using approach-to-equilibrium molecular dynamics (AEMD). The AEMD method has been developed to calculate the thermal conductivity of periodic bulk systems subjected to

a one-dimensional (1D) heat flux. Within these assumptions, thermal conductivity obtained is contingent on the period, and increasing this period leads to an increase in thermal conductivity until a plateau at the macroscopic value is reached. For aSiO₂, saturation occurs at 6 nm, at a value of (1.2 ± 0.1) W K⁻¹ m⁻¹, in quantitative agreement with measurements.

In the present study, we extend the application of AEMD to aSiO₂ films. First-principles molecular dynamics (FPMD) is the computational tool used to preserve the predictive power shown for bulk aSiO₂ at the density functional theory (DFT) level of accuracy. Specifically, the restriction represented by periodic boundary conditions (PBC) is lifted in the direction perpendicular to the film plane. In the first stage, the surfaces are meticulously equilibrated. Subsequently, the in-plane thermal conductivity is calculated to ascertain whether the free surfaces can impact heat transfer. The calculation is performed for a 2 nm-thick film. In a subsequent step, we focus on the cross-plane thermal conductivity. This calculation requires an extension of the AEMD to non-periodic systems which, in our case, translates into ultra-thin films with thicknesses of 1 and 2 nm.

The paper is organized as follows. In Section 2, the AEMD approach is developed within the framework of cross-plane heat transport. A comparative analysis is then conducted with its derivation in the context of periodic material or in-plane heat transport. The FPMD framework and the atomistic models of the aSiO₂ films are described in

<footnote>☆ This article is part of a Special issue entitled: 'Computational modelling of materials' published in Solid State Sciences.</footnote>

<footnote>* Corresponding author.
E-mail address: evelyne.martin@cnrs.fr (E. Martin).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Configurations used to extract from the AEMD approach the in-plane (left) and cross-plane (right) thermal conductivity of a film of thickness $L_z$. The red and blue areas indicate the hot and cold blocks ($T_1 > T_2$), respectively.</figcaption></figure>

Section 3. Then, the AEMD is applied to in-plane heat transport, and the results are presented in Section 4. The issue of cross-plane thermal conductivity of aSiO$_2$ films is addressed in Section 5. Finally, a number of findings is provided in relation to both aSiO$_2$ and the utilization of AEMD for the estimation of cross-plane thermal conductivities.

## 2. General assessment of the in-plane and cross plane thermal conductivities

As a template model, we consider a homogeneous system on which periodic boundary conditions are applied along the x and y directions, namely a film of finite size in the z direction (Fig. 1). To extract the thermal conductivity, a thermal transient regime is initiated in one direction.

In-plane. For the in-plane case, the thermal transport is assumed to occur along the y direction (left part of Fig. 1). The heat conduction in the y direction is governed by the heat equation:

$$ \frac{\partial T(y, t)}{\partial t} = \frac{\kappa_{in}}{C \rho} \frac{\partial^2 T(y, t)}{\partial y^2} \qquad (1) $$

where $\kappa_{in}$ is the in-plane thermal conductivity in W K$^{-1}$ m$^{-1}$, C the heat capacity in J K$^{-1}$ and $\rho$ the number density in m$^{-3}$. The equation is then solved as in the case of a bulk material [4], using the continuity of the temperature field and heat flux $\vec{J} = -\kappa_{in} \nabla T$ across the periodic boundaries at $y = 0$ and $y = L_y$. The analytical solution is the Fourier series:

$$ T(y, t) = \frac{T_1 + T_2}{2} + \sum_{n=1}^{\infty} \frac{2(T_1 - T_2)}{(2n+1)\pi} \sin\left(\frac{2\pi(2n+1)y}{L_y}\right) e^{-(2n+1)^2 t/\tau_{in}} \qquad (2) $$

with

$$ \tau_{in} = \frac{L_y^2 C \rho}{4\pi^2 \kappa_{in}} \qquad (3) $$

The dominant contribution, a sinusoidal temperature profile of period $L_y$, has been observed in several former AEMD applications (see for instance [5–7]). The corresponding difference $\Delta T(t)$ between the average temperature of the hot and cold blocks is equal to:

$$ \begin{aligned} \Delta T(t) &= \frac{2}{L_y} \int_{L_y/2}^{L_y} T(y,t) dy - \frac{2}{L_y} \int_{0}^{L_y/2} T(y,t) dy \\ &= \sum_{n=1}^{\infty} \frac{8(T_1 - T_2)}{(2n+1)^2 \pi^2} e^{-(2n+1)^2 t/\tau_{in}}. \end{aligned} \qquad (4) $$

Cross-plane. The heat conduction in the z direction is governed by a heat equation very similar to the former one, namely:

$$ \frac{\partial T(z, t)}{\partial t} = \frac{\kappa_{cr}}{C \rho} \frac{\partial^2 T(z, t)}{\partial z^2} \qquad (5) $$

being $\kappa_{cr}$ the cross-plane thermal conductivity. In this case, the heat flux vanishes at the top and bottom surfaces because of the vacuum, resulting in the Dirichlet-like boundary conditions:

$$ \left. \frac{\partial T(z, t)}{\partial t} \right|_{z=0} = 0 \quad \text{and} \quad \left. \frac{\partial T(z, t)}{\partial t} \right|_{z=L_z} = 0. \qquad (6) $$

The solution is then expressed by the series:

$$ T(z, t) = \frac{A_0}{2} + \sum_{n=1}^{\infty} A_n \cos\left(\frac{n\pi z}{L_z}\right) e^{-n^2 t/\tau_{cr}} \qquad (7) $$

with

$$ \tau_{cr} = \frac{L_z^2 C \rho}{\pi^2 \kappa_{cr}}. \qquad (8) $$

The Fourier coefficients

$$ A_n = \frac{2}{L_z} \int_{z=0}^{L_z} T(z, 0) \cos\left(\frac{n\pi z}{L_z}\right) dz \qquad (9) $$

can be calculated by applying the initial conditions given in the right part of Fig. 1, i. e. $T(z, 0) = T_1$ for $0 < z < L_z/2$ and $T(z, 0) = T_2$ for $L_z/2 < z < L_z$. We arrive then at the analytical solution:

$$ T(z, t) = \frac{T_1 + T_2}{2} + \sum_{n=1}^{\infty} \frac{4(T_1 - T_2)}{n\pi} \cos\left(\frac{n\pi z}{L_z}\right) e^{-n^2 t/\tau_{cr}}. \qquad (10) $$

The difference $\Delta T(t)$ between the average temperature of the hot and cold blocks is equal to:

$$ \Delta T(t) = \sum_{n=1}^{\infty} \frac{8(T_1 - T_2)}{n^2 \pi^2} e^{-n^2 t/\tau_{cr}}. \qquad (11) $$

Summary. A feature common to both cases is the exponential trend of the approach-to-equilibrium. More precisely, this equilibration occurs through a series of decreasing exponentials as explicitly given above. In the in-plane case, the decay time of order 2 is nine times smaller than that of order 1, whereas in the cross-plane case, the ratio between the two first decay times is equal to four, clearly smaller with respect to the in-plane case. On these grounds, we can infer that it should be easier to observe the second order in the cross-plane case. As for the temperature profile, it is sinusoidal with the period determined by the PBC ($L_y$) in the in-plane case and cosinusoidal with a period of 2 times the film thickness ($2L_z$) in the cross-plane case, to cope with the zero flux conditions at the surfaces. The relationships for obtaining the thermal conductivities in both cases are:

$$ \kappa_{in} = \frac{L_y^2 C \rho}{4\pi^2 \tau_{in}} \quad \text{and} \quad \kappa_{cr} = \frac{L_z^2 C \rho}{\pi^2 \tau_{cr}}. \qquad (12) $$

## 3. Models of aSiO$_2$ films

To apply Eq. (12) to the targeted system, FPMD is employed to generate heat flow in an amorphous SiO$_2$ film of thickness $d$. The selection of this material is driven by considerations outlined in the introduction, as well as the prior application of the AEMD method to calculate the thermal conductivity of bulk aSiO$_2$. The results obtained for bulk aSiO$_2$ as a function of the size of the systems, which is also equal to the period of the temperature profile, are shown in Fig. 5. The obtained values increase with the period, until a saturation is reached at 6 nm, with a value that is in agreement with available experimental data (highlighted as a pink area in the figure). In the case of the in-plane thermal conductivity of the film, our study targets a system with a period equal to 8 nm. This ensures that if a lower thermal conductivity
<table><caption>**Table 1**

Characteristics of the models of amorphous SiO₂ films used in this work.</caption><thead><tr><th>Name</th><th>N<sub>at</sub></th><th>Dimensions of the supercell (Å<sup>3</sup>)</th><th>Film thickness d (Å)</th></tr></thead><tbody><tr><td>Film S1</td><td>279</td><td>20.32 × 20.32 × 21.87</td><td>10.2 ± 1.0</td></tr><tr><td>Film S2</td><td>284</td><td>10.16 × 20.32 × 32.00</td><td>20.3 ± 1.0</td></tr><tr><td>Film L2</td><td>1114</td><td>10.16 × 81.28 × 32.00</td><td>20.3 ± 0.5</td></tr></tbody></table>

is obtained, it will be due to the presence of surfaces facing a vacuum, independently of the size effects that occur below 6 nm.

The amorphous silica model systems used in [3] were obtained with a well assessed melt/quench procedure. The resulting structure, particularly the pair correlation functions, closely resembles what reported in the literature. In this work, the aforementioned structures are used as a starting model system to expand the simulation cell in one direction, with the scope of introducing surfaces. The empty space added to the system (10 Å of vacuum) ensures that electronic wave functions and densities protruding from the surfaces vanish (<10⁻⁵ a.u.) well before the next periodic image. This provides good separation within the applied periodic boundary conditions scheme. The film thus generated is then relaxed to its equilibrium configuration prior to being brought to room temperature. The FPMD specific implementation used in this work is the one proposed by Car and Parrinello [8], as implemented in the OpenCPMD code [9], and relying on the DFT formulation for the potential energy of the system. The exchange-correlation energy of the Kohn-Sham functional is composed by the exchange formula proposed by Becke [10] and the correlation one of Lee and coworkers [11] (BLYP). The valence-core interaction is described by numerical norm-conserving pseudopotentials of Trouiller and Martins [12]. The wave functions of the valence electrons are represented by a plane wave basis set with a relatively high cutoff of 80 Ry (1088.4 eV) to accurately depict the vanishing wave functions at the free surfaces facing the vacuum. The sampling of the Brillouin zone is restricted to the Γ point. The mass of the fictitious electronic degrees of freedom is set at 600 a.u., and the time step at 5 a.u. (0.12 fs) to ensure optimal conservation of the constants of motion. The ionic temperature is controlled by a Nosé-Hoover thermostat [13–15], while the fictitious electronic kinetic energy is controlled with Blöchl-Parrinello thermostat [16]. The masses of the thermostats (in units of corresponding frequencies) are equal to 580 cm⁻¹ for ions and 800 cm⁻¹ for fictitious electronic degrees of freedom.

The electron density is modified by the presence of the two surfaces, which cause surface atoms to displace and reform bonds that were broken when a vacuum was introduced into the simulation cell. Some atoms, particularly O atoms, detach from the surface and, therefore, were identified and removed from the system. The corresponding stoichiometric variation, less than 2 %, is considered negligible in the following. After a further relaxation of the *cleaned* configuration, the system was brought to room temperature. We remark that upon removal of the electronic thermostat, an increase of the fictitious kinetic energy occurs, accompanied by a decrease in ionic temperature. This is problematic since in Phase 2 of the AEMD, the system must be entirely free from any thermostat, and the dynamics must take place in the microcanonical (N, V, E) ensemble. To investigate the origin of this loss of adiabaticity, we analyzed the electronic density of states, identifying those states within the band gap that could affect the adiabaticity. According to the protocol assessed in Ref. Lambrecht et al. [5], the inverse participation ratio of the states was computed to identify dangling bonds (DBs), being a general source of instability in first-principles calculations. To get rid of the DBs, as routinely done in these cases, we passivated then with H atoms. As a consequence, for the case of the system of length equal to 8 nm, the number of atoms in the system varied from N<sub>at</sub> = 1116 to 1114 to include the 5 additional H atoms. The characteristics of this film labeled L2 are summarized in Table 1.

The same model film L2 was used to compute the cross-plane thermal conductivity, along with two other systems listed in Table 1.

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Atomistic representation of the three models systems of SiO₂ films.</figcaption></figure>

The systems discussed here are labeled as S1 and S2; both are constructed from our bulk supercell containing 279 atoms and extensively documented in Ref. Martin et al. [3]. The only difference is a rotation applied to system S1, resulting in a film with a thickness d = 1 nm, whilst S2 and L2 model films have a thickness d = 2 nm. In the case S2, five DBs were passivated by H atoms, thus bringing the total number of atoms to 284. This additional step was not needed for the surfaces of system S1. The systems S1, S2 and L2 are illustrated in Fig. 2.

## 4. In-plane thermal conductivity of aSiO₂ films

As mentioned, the thermal conductivity of a film with a thickness of d = 2 nm is studied using the model L2. The system is initially placed in an out-of-equilibrium state through the simultaneous application of two thermostats at temperatures of T = 200 and 400 K applied to two blocks delimited along the y-direction (0 < y < L<sub>y</sub>/2 and L<sub>y</sub>/2 < y < L<sub>y</sub>). As illustrated in Fig. 3.a, the temperatures of these blocks are stabilized during an initial stage refereed to as Phase 1, lasting generally from 2 to 5 ps.

In a second step, termed Phase 2, the thermostats are turned off to allow the system to approach equilibrium. Fig. 3.a shows that, for the system under consideration, it takes approximately 6 ps for the temperature difference between the hot and cold blocks to become negligible. Fig. 3.b reports this temperature difference on a semi-log scale, evidencing a clear exponential decrease, according to the main contribution of Eq. (4). The temperature profile, averaged over the total duration of Phase 2, is shown in Fig. 3.c and assumes a sinusoidal shape with a period equal to the size L<sub>y</sub> in the heat transport direction. Oscillations superimpose on the temperatures of the hot and cold blocks and are visible in Fig. 3.a. These oscillations are due to a vibrational mode in the (x, y) plane, as shown in Fig. 4 by comparing the positions of the atoms at the beginning of Phase 2 and at a time corresponding to the maximum amplitude of the oscillation. Since PBC are applied in these two directions, the thin film structure remains unaffected, and the effect is only visible in the temperatures of the blocks.

The temperature difference during the approach to equilibrium is adjusted to extract the decay time τ<sub>in</sub>, which turns out to be 4.0 ps in the example of Fig. 3. Eq. (12) linking τ<sub>in</sub> to κ<sub>in</sub> is then applied to extract the in-plane thermal conductivity. The heat capacity, C = dU/dT, entering into the equation, is determined from the values of the total energy U at different temperatures and results to be equal to 3N<sub>at</sub>k<sub>B</sub> within an error bar of 1%. In order to estimate the statistical variation of the thermal conductivity calculated according to these guidelines, two auxiliary calculations were done. These involved extending Phase 1 by a few ps and then switching to Phase 2. This leads to the determination of a second value of τ<sub>in</sub> and therefore of κ<sub>in</sub>. This procedure is then repeated to obtain a third value of κ<sub>in</sub>. These three values are shown in Fig. 5 by small green dots. Their average value and associated standard deviation are 1.35 ± 0.14 W K<sup>-1</sup> m<sup>-1</sup>, and are shown as a large dot and an associated error bar in the figure. The individual and average values are also reported for completeness in Table 2.

Compared to our former calculations [3] on systems periodic in the three directions, we see here for the film that the thermal conductivity remains basically unaltered. Hence, we can infer that heat transfer is not affected by the presence of surfaces, even though the layer
<figure><img src="image_6.png" /><figcaption>Fig. 3. AEMD for the in-plane configuration: (a) temperature of the hot and cold blocks during Phases 1 and 2; (b) difference of temperature between the two blocks during Phase 2 (black line) and exponential fit (green line) ($\Delta T_0 \exp(-t/\tau)$, with $\tau = 4.0$ ps); (c) temperature profile during Phase 2 (black dots) and sinusoidal fit (purple line). System L2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Schematic atomic representation of the SiO₂ film at the beginning ($t = 0$ ps) of Phase 2 of the AEMD, and after 2.6 ps. Si atoms are drawn in beige, O in red and H in white.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. In-plane film (green dots) and bulk thermal conductivity (grey dots) (reproduced from Ref. Martin et al. [3]) of amorphous SiO₂ versus the period of the temperature profile. The small dots represent values from different Phases 2, while the large dots and their error bars represent the mean value and its corresponding standard deviation. Case of film of thickness equal to $d = 2$ nm, and a cross section equal to $S = 2$ nm² (Film L2). The pink area represents the range [1.1–1.5] W K⁻¹ m⁻¹ of experimental values for the bulk reported in the literature [17–21].</figcaption></figure>

<caption>Table 2</caption>

In-plane thermal conductivities (W K⁻¹ m⁻¹) of the 2 nm-thick amorphous SiO₂ film.

<table><tr><td>1st determination</td><td>1.42</td></tr><tr><td>2nd determination</td><td>1.15</td></tr><tr><td>3rd determination</td><td>1.47</td></tr><tr><td>Average</td><td>1.35 ± 0.14</td></tr></table>

is only about 2 nm thick. This result differs from what has been observed, for instance, in a silicon film [22] where thermal conductivity is reduced by a factor of 4 to 5 when moving from the bulk material to an 8 nm-thick film. This difference can be attributed to the distinct thermal transport mechanisms observed in aSiO₂ and Si. In silicon, heat conduction occurs primarily through propagative modes (phonons). In contrast, in silica, heat is predominantly transmitted through non-propagative modes, by thermal diffusion.

## 5. Cross-plane thermal conductivity of aSiO₂ films

The three systems listed in Table 1 are used to investigate the cross-plane thermal conductivity. For each one, the AEMD Phase 1 is conducted through the delimitation of hot and cold blocks in the z direction. The temperature profile obtained during Phase 2 in the case of the system L2 is displayed in Fig. 6.a, along with a section of cosine with a period equal to 4 nm, amounting to twice the film thickness $d = L_z = 2$ nm ($\cos(2\pi z/(2L_z)) = \cos(\pi z/L_z)$). The strong correspondence between the FPMD outcome and the cosine fit validates the hypothesis that our results follow Eq. (10). During this Phase 2, the temperature difference between the hot and cold blocks decreases exponentially, as shown in Fig. 6.b.

The value of $\tau_{cr}$ can be extracted from this decrease, thus allowing for the calculation of $\kappa_{cr}$ via the relationship reported on the right-hand side of Eq. (12). The results obtained for the three systems are shown in Fig. 7. The data points represent the values obtained by repeatedly executing several AEMD, specifically three for the system L2 and five for the systems S1 and S2. The larger dots, along with their error bars, denote the mean value and standard deviation, respectively. An average
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** AEMD in the cross-plane configuration: (a) temperature profile during Phase 2 (black dots) and cosinusoidal fit (pink line).; (b) difference of temperature between the hot and cold blocks during Phase 2 (black line) and exponential fit (green line) ($\Delta T_0 \exp(-t/\tau)$, with $\tau = 1.3$ ps). System L2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Cross-plane thermal conductivity as a function of the aSiO₂ film thickness *d*, for system S1 (purple dots), S2 (blue dots) and L2 (red dots), and bulk thermal conductivity versus half-length of periodicity (grey dots) (reproduced from Ref. Martin et al. [3]). The small dots represent values from different Phases 2, while the large dots and their error bars are the average values and corresponding standard deviations.</figcaption></figure>

<caption>**Table 3**</caption>

Cross-plane thermal conductivities (W K⁻¹m⁻¹) of the 1 nm- and 2 nm-thick amorphous SiO₂ films.

<table><thead><tr><th rowspan="2"></th><th colspan="1">d = 1 nm</th><th colspan="2">d = 2 nm</th></tr><tr><th>S1</th><th>S2</th><th>L2</th></tr></thead><tbody><tr><td>1st determination</td><td>0.39</td><td>0.89</td><td>0.90</td></tr><tr><td>2nd determination</td><td>0.66</td><td>1.09</td><td>0.82</td></tr><tr><td>3rd determination</td><td>0.43</td><td>0.89</td><td>0.89</td></tr><tr><td>4th determination</td><td>0.51</td><td>1.02</td><td>-</td></tr><tr><td>5th determination</td><td>0.44</td><td>1.17</td><td>-</td></tr><tr><td>Average</td><td>0.49 ±0.10</td><td>1.01 ± 0.11</td><td>0.87 ± 0.03</td></tr></tbody></table>

value and an error bar are also provided for the film thickness $d = L_z$ to account for the roughness of the film surface of our models. The individual and average values are also reported for completeness in Table 3. As evidenced in Fig. 7, the cross-plane thermal conductivity of a 2 nm thick film is lower than that of bulk aSiO₂, and this decrease is even more pronounced for a 1 nm thick film. It is also evident that the cross-section does not exert a substantial influence on the result. However, an increased cross-section induces a reduction in the error bar because of the increased number of atoms per block, which improves statistics, thus rendering temperature fluctuations less pronounced.

For completeness, Fig. 7 includes also results reported in Ref. Martin et al. [3] for bulk aSiO₂ as a function of the half-period of the temperature profiles. This choice is motivated by the fact that the boundary conditions in the film require the temperature profile to be a cosine whose half-period is equal to the thickness of the film. We

stress the fact that this is a purely geometrical requirement independent on the chemical nature of the atoms composing the system. The result is a fairly good agreement of the cross-plane thermal conductivities of the films and of the bulk. On these grounds, we can infer that for aSiO₂ films of thickness larger than 3 nm, the cross-plane thermal conductivity will not decrease with respect to the macroscopic value.

A further noteworthy finding that emerges from these calculations is that the thermal conductivities $\kappa^{\text{Bulk}}(L)$ obtained by AEMD as a function of the period $L$ of systems periodic in the three directions are sufficient to characterize the reduction of cross-plane thermal conductivity of films $\kappa^{\text{Film}}(d)$, using the relation:

$$ \kappa_{cr}^{\text{Film}}(d) = \kappa^{\text{Bulk}}(2d). \quad (13) $$

A key benefit of this approach is that it eliminates the need for re-calculations, and the additional steps required to eliminate DBs and surface states that could possibly jeopardize the adiabaticity of the FPMD simulations.

## 6. Conclusion

The present study focused on the thermal conductivities of films, intended as an extension of the AEMD method we developed and used for bulk materials. The amorphous silica has been selected as a prototype material to assess this extension in view of both its importance in practical applications and the availability of reliable experimental data. This same material has been formerly targeted [3], thus providing a good template for the proposed extension. At variance with the fully periodic bulk previously studied, here we apply periodic boundary conditions in only two directions, thereby creating two free surfaces. The reconstruction of these surfaces was carefully analyzed and existing DBs passivated by hydrogen, as routinely done both in simulations and synthesis. The subsequent step involved the simulation of temperature transients. In the initial phase, heat flow occurred in the plane, thereby resulting in a sinusoidal temperature profile with a period equal to the size of the system. The thermal conductivity was determined from the time required by the system to revert back to equilibrium conditions. The outcome of our simulations indicate that no reduction occurs for the 2 nm-thick thin film. We could then deduce that this lack of sensitivity to the presence of surfaces is related to the diffusion mode characterizing conduction in amorphous silica.

In a subsequent step, the thermal transient has been generated in the thickness of the film to obtain the cross-plane thermal conductivity. A specifically suited equation has been formulated to relate the transient decay time to the cross-plane thermal conductivity, thereby accounting for the cancellation of the heat flux at the surfaces facing the vacuum. The application of this relation has enabled the determination of the cross-plane thermal conductivity of films with thicknesses of 2 and 1 nm. In comparison to the macroscopic thermal conductivity of glassy silica, these values exhibit a decrease. Moreover, the cross-plane thermal conductivity of a film of thickness $d$ has been shown to be equivalent to that of the bulk material of length $2d$ in the heat transport direction as characterized by AEMD. We can infer from this finding that the cross-plane thermal conductivity of an aSiO₂ film will only be reduced for thicknesses less than 3 nm.
CRediT authorship contribution statement

**Louis Petit:** Investigation. **Cheick O. Diarra:** Methodology. **Achille Lambrecht:** Methodology. **Ibrahim Bel-Hadj:** Methodology. **Mauro Boero:** Software, Investigation. **Evelyne Martin:** Writing – original draft, Supervision, Investigation, Conceptualization.

Code availability statement

The Open CPMD code is available under MIT License at https://github.com/OpenCPMD

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Evelyne Martin reports financial support was provided by French National Research Agency. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the French ANR via the TACOS project (ANR-24-CE-09-5163). The authors acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsac calcul/Big Data. Calculations on the larger systems were performed by using resources from GENCI (Grand Equipement National de Calcul Intensif) (grant no. 0910296 and 0905071).

Data availability

Data will be made available on request.

References

[1] M.C. Wingert, J. Zheng, S. Kwon, R. Chen, Thermal transport in amorphous materials: a review, Semicond. Sci. Technol. 31 (11) (2016) 113003, http://dx.doi.org/10.1088/0268-1242/31/11/113003.

[2] W. Zhu, G. Zheng, S. Cao, H. He, Thermal conductivity of amorphous SiO₂ thin film: A molecular dynamics study, Sci. Rep. 8 (2018) http://dx.doi.org/10.1038/s41598-018-28925-6.

[3] E. Martin, G. Ori, T.-Q. Duong, M. Boero, C. Massobrio, Thermal conductivity of amorphous SiO₂ by first-principles molecular dynamics, J. Non-Cryst. Solids 581 (2022) 121434, http://dx.doi.org/10.1016/j.jnoncrysol.2022.121434.

[4] E. Lampin, P.L. Palla, P.-A. Francioso, F. Cleri, Thermal conductivity from approach-to-equilibrium molecular dynamics, J. Appl. Phys. 114 (3) (2013) 033525, http://dx.doi.org/10.1063/1.4815945.

[5] A. Lambrecht, G. Ori, C. Massobrio, M. Boero, E. Martin, Assessing the thermal conductivity of amorphous SiN by approach-to-equilibrium molecular dynamics, J. Chem. Phys. 160 (9) (2024) 094505, http://dx.doi.org/10.1063/5.0193566.

[6] I. Bel-Hadj, M. Guerboub, A. Lambrecht, G. Ori, C. Massobrio, E. Martin, Thermal conductivity of crystalline $Ge_2Sb_2Te_5$; lattice contribution and size effects in the cubic phase quantified by approach-to-equilibrium molecular dynamics, J. Phys. D: Appl. Phys. 57 (23) (2024) 235303, http://dx.doi.org/10.1088/1361-6463/ad316b.

[7] C.O. Diarra, C. Massobrio, E. Martin, Heat transfer through covalent versus non-covalent bonding: a case study on crystalline π-conjugated P3HT polymer using approach-to-equilibrium molecular dynamics, Phys. Chem. Chem. Phys. 27 (2025) 18487-18495, http://dx.doi.org/10.1039/D5CP02084K.

[8] R. Car, M. Parrinello, Unified approach for molecular dynamics and density-functional theory, Phys. Rev. Lett. 55 (1985) 2471-2474, http://dx.doi.org/10.1103/PhysRevLett.55.2471.

[9] CPMD, Copyright IBM Corp. 1990-2022, copyright MPI für Festkörperforschung Stuttgart 1997-2001, 2025, Available on GitHub under MIT License. URL: https://github.com/OpenCPMD.

[10] A.D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38 (1988) 3098-3100, http://dx.doi.org/10.1103/PhysRevA.38.3098.

[11] C. Lee, W. Yang, R.G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 37 (1988) 785-789, http://dx.doi.org/10.1103/PhysRevB.37.785.

[12] N. Troullier, J.L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B 43 (1991) 1993-2006, http://dx.doi.org/10.1103/PhysRevB.43.1993.

[13] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1) (1984) 511-519, http://dx.doi.org/10.1063/1.447334.

[14] W.G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31 (1985) 1695-1697, http://dx.doi.org/10.1103/PhysRevA.31.1695.

[15] C. Massobrio, I.A. Essomba, M. Boero, C. Diarra, M. Guerboub, K. Ishisone, A. Lambrecht, E. Martin, I. Morrot-Woisard, G. Ori, C. Tugène, S.D. Wansi Wendji, On the actual difference between the Nosé and the Nosé-Hoover thermostats: A critical review of canonical temperature control by molecular dynamics, Phys. Status Solidi (b) 261 (1) (2024) 2300209, http://dx.doi.org/10.1002/pssb.202300209.

[16] P.E. Blöchl, M. Parrinello, Adiabaticity in first-principles molecular dynamics, Phys. Rev. B 45 (1992) 9413-9416, http://dx.doi.org/10.1103/PhysRevB.45.9413.

[17] R. Anufriev, S. Tachikawa, S. Gluchko, Y. Nakayama, T. Kawamura, L. Jalabert, M. Nomura, Cross-plane thermal conductivity in amorphous Si/SiO₂ superlattices, Appl. Phys. Lett. 117 (9) (2020) 093103, http://dx.doi.org/10.1063/5.0023782.

[18] J.L. Braun, C.H. Baker, A. Giri, M. Elahi, K. Artyushkova, T.E. Beechem, P.M. Norris, Z.C. Leseman, J.T. Gaskins, P.E. Hopkins, Size effects on the thermal conductivity of amorphous silicon thin films, Phys. Rev. B 93 (2016) 140201, http://dx.doi.org/10.1103/PhysRevB.93.140201.

[19] D.G. Cahill, M. Katiyar, J.R. Abelson, Thermal conductivity of a-Si:H thin films, Phys. Rev. B 50 (1994) 6077-6081, http://dx.doi.org/10.1103/PhysRevB.50.6077.

[20] M. Okuda, S. Ohkubo, A novel method for measuring the thermal conductivity of submicrometre thick dielectric films, Thin Solid Films 213 (2) (1992) 176-181, http://dx.doi.org/10.1016/0040-6090(92)90280-O.

[21] K. Goodson, M. Flik, L. Su, D. Antoniadis, Annealing-temperature dependence of the thermal conductivity of LPCVD silicon-dioxide layers, IEEE Electron Device Lett. 14 (10) (1993) 490-492, http://dx.doi.org/10.1109/55.244740.

[22] H. Zaoui, P.L. Palla, S. Giordano, F. Cleri, M. Verdier, D. Lacroix, J.-F. Robillard, K. Termentzidis, E. Martin, Thermal conductivity of deca-nanometric patterned Si membranes by multiscale simulations, Int. J. Heat Mass Transfer 126 (2018) 830-835, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.06.004.