MnLargeSymbols’164 MnLargeSymbols’171

# Thermal conductivity of glasses above the plateau: first-principles theory
and applications

Michele Simoncelli<EMAIL_ADDRESS>Theory of Condensed Matter Group of the
Cavendish Laboratory, University of Cambridge (UK)  Francesco Mauri
Dipartimento di Fisica, Università di Roma La Sapienza, Italy  Nicola Marzari
Theory and Simulation of Materials (THEOS) and National Centre for
Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

###### Abstract

Predicting the thermal conductivity of glasses from first principles has
hitherto been a prohibitively complex problem. In fact, past works have
highlighted challenges in achieving computational convergence with respect to
length and/or time scales using either the established Allen-Feldman or Green-
Kubo formulations, endorsing the concept that atomistic models containing
thousands of atoms — thus beyond the capabilities of first-principles
calculations — are needed to describe the thermal conductivity of glasses. In
addition, these established formulations either neglect anharmonicity (Allen-
Feldman) or miss the Bose-Einstein statistics of atomic vibrations (Green-
Kubo), thus leaving open the question on the relevance of these effects. Here,
we present a first-principles formulation to address the thermal conductivity
of glasses above the plateau, which can account comprehensively for the
effects of structural disorder, anharmonicity, and quantum Bose-Einstein
statistics. The protocol combines the Wigner formulation of thermal transport
with convergence-acceleration techniques, and is validated in vitreous silica
using both first-principles calculations and a quantum-accurate machine-
learned interatomic potential. We show that models of vitreous silica
containing less than 200 atoms can already reproduce the thermal conductivity
in the macroscopic limit and that anharmonicity negligibly affects heat
transport in vitreous silica. We discuss the microscopic quantities that
determine the trend of the conductivity at high temperature, highlighting the
agreement of the calculations with experiments in the temperature range above
the plateau where radiative effects remain negligible
($50{\lesssim}T{\lesssim}450$ K).

##  I Introduction

The thermal conductivity of glasses is a key property for many and diverse
technological applications, including, e.g., the microelectronics 1,
pharmaceutical 2, aerospace 3, 4, optic5, and construction 6 industries. To a
large degree of universality, for glasses there are three characteristic
intervals in the temperature dependence of the thermal conductivity
$\kappa(T)$ 7, 8, 9, 10: (i) the low temperature region ($T{\lesssim}1\;K$)
where $\kappa(T){\sim}T^{2}$; (ii) the intermediate temperature region
($5{\lesssim}T{\lesssim}25\;K$) where the conductivity displays a plateau
($\kappa(T){\sim}{\rm constant}$); (iii) the above-the-plateau region
($T{\gtrsim}30$ K) where the thermal conductivity increases with temperature.
The trend observed in region (i) is usually explained relying on the
phenomenological 7, 11, 12, 13 “tunneling two-level system” (TLS) model 14,
15. The plateau found in the region (ii) has motivated several studies, which
explored its possible connection with the boson peak 16, 17, 18, 19. The
above-the-plateau region (iii) — which is the most relevant one for the
aforementioned technological sectors — has up to now been investigated using
the Allen-Feldman (AF) formulation 8 or molecular dynamics (MD) methods (the
Green-Kubo formulation in combination with classical 20, 21, 22, 23, first-
principles 24, 25, 26, 27, or machine learning 28 MD simulations (GKMD), non-
equilibrium MD (NEMD) 29, 20, 30, or the approach-to-equilibrium MD (AEMD) 31,
32). Past works relying on these formulations have highlighted two major
challenges. First, achieving computational convergence with respect to
supercell sizes and simulations times is onerous and problematic 33, 20, 34,
leading to the conclusion that accurate calculations require models containing
thousands of atoms, and simulation times of the order of tens to hundreds of
picoseconds; this yields a computational cost that is prohibitively high for
direct first-principles calculations. Second, neither Allen-Feldman nor MD
formulations can immediately ensure a comprehensive description of thermal
transport, since the former neglects anharmonicity (thus it is in principle
accurate only in the low-temperature regime where anharmonicity phases out 9),
and the latter, while accounting for anharmonicity exactly, misses the quantum
Bose-Einstein statistics of atomic vibrations 35, relevant at low
temperatures.

The recently developed Wigner formulation 36, 37 has shown that the two
established microscopic heat-conduction mechanisms for crystals and glasses —
i.e. the particle-like propagation of phonon wavepackets described by the
Peierls-Boltzmann equation 38, 39, and the wave-like tunneling mechanisms
described by the AF equation8, 9, 111here the “diffuson” mechanisms described
by Allen-Feldman theory is denoted with “wave-like tunneling” because it bears
some analogies to the electronic Zener interband tunneling 59, more details
about the terminology can be found in Refs. 36, 37., respectively — both
emerge as limiting cases from a unified transport equation. Such a Wigner
formulation offers a comprehensive approach to describe heat transport in
solids, encompassing: (i) “simple crystals”, characterized by phonon interband
spacings much larger than the linewidths, where particle-like propagation
dominates and the Peierls-Boltzmann limit is accurate; (ii) glasses, where
wave-like tunneling is relevant and the Wigner formulation extends AF theory
accounting for anharmonicity; (iii) the intermediate regime of “complex
crystals”, where particle-like propagation and wave-like tunneling are
simultaneously present and neither Peierls-Boltzmann or Allen-Feldman are
sufficient. This intermediate regime is common and prevalent crystals
characterized by phonon interband spacings smaller than the linewidths and
featuring ultralow thermal conductivity (e.g. those used in thermoelectrics41
or thermal barrier coatings37). It is worth mentioning that the thermal
conductivity expression derived from the Wigner formulation has been obtained
also following a many-body Green’s function approach33, 42, 43. The Wigner
formulation for transport has paved the way, as we will see, to tackle the
aforementioned challenges of comprehensively describing transport in glasses
accounting for the interplay between Bose-Einstein statistics, anharmonicity,
and disorder. Still, the challenge of predicting the conductivity of glasses
from first principles remains, due to a structural complexity larger than
simple crystals 44, disordered crystals45, and complex crystals, for which
several advances have recently been made 36, 46, 41, 37, 47.

Here, we present a protocol that addresses both the aforementioned challenges,
enabling accurate first-principles predictions for the thermal conductivity of
glasses, and combining the Wigner transport equation36, 37 (WTE) with
convergence-acceleration techniques; it is first discussed in the limiting
case of vanishing anharmonicity, where the WTE conductivity reduces to the
harmonic AF8, 9 conductivity, and then extended to account for the effects of
anharmonicity and evaluate from finite-size models of glasses the bulk limit
of the anharmonic WTE conductivity. We showcase the protocol in vitreous
silica ($v$-SiO2) — a paradigmatic glass that is employed in many and diverse
technologies 1, 5, 2, 4, 3, 6 — extending the simulations to very large cell
sizes thanks to the use of machine learned potentials (the recently developed
Gaussian Approximation Potential48 (GAP) for silica polymorphs 49). The GAP
potential we employ has been generated from first-principles calculations
performed using the SCAN functional50, and yields quantum-accurate predictions
for the vibrational properties of various silica polymorphs49. Here, we employ
GAP to describe the vibrational properties and evaluate the thermal
conductivity — both in the harmonic Allen-Feldman or in the anharmonic Wigner
framework — of models of $v$-SiO2 having very different sizes, showing how our
protocol allows to accurately evaluate the bulk limit of the harmonic or
anharmonic conductivity using models containing less than 200 atoms, that are
in very good agreement with the macroscopic cell limit of 5184 atom 49.

After having validated the protocol, we employ it to study the conductivity of
$v$-SiO2 fully from first principles using models of $v$-SiO2 containing
10851, 14452, 53, and 192 atoms54. We discuss how the AF or WTE conductivities
change if the widely used semi-empirical van Beest-Kramer-van Santen55, 56
(BKS) potential, or the state-of-the-art GAP potential for silica polymorphs
49, are used in place of first-principles calculations to compute the
vibrational properties of $v$-SiO2. We elucidate analogies and differences
between the harmonic and the anharmonic conductivity, relying on the latter to
show how the high-temperature trend of $\kappa(T)$ is determined by the
variation of the velocity-operator elements, whose value regulates the amount
of heat transferred by wave-like tunneling between vibrational eigenstates,
when the energy difference between coupled eigenstates increases. Finally, we
rationalize these results at the microscopic level, relying on Allen-Feldman’s
harmonic diffusivity and extending such quantity to account for anharmonicity,
but also showing that anharmonicity has negligible effects on the conductivity
and diffusivity of $v$-SiO2.

##  II Wigner formulation of thermal transport

As anticipated, the WTE formalism is general and can be used to describe
solids with variable degrees of disorder, ranging from ordered crystals with
small primitive cells, to disordered glasses with diverging primitive cell (in
this latter case, for sufficiently large primitive cells periodic boundary
conditions become irrelevant and the Brillouin zone (BZ) reduces to the point
$\bm{q}{=}\bm{0}$ only). In practice, in numerical simulations non-periodic
glasses can be approximatively described as crystals with large but finite
primitive cells, thus having a small but finite BZ that includes wavevectors
different from $\bm{q}{=}\bm{0}$. We will discuss in the next section the
conditions under which non-periodic glasses can be simulated in periodic-
boundary conditions, after having introduced here the salient features of the
WTE formulation in the general case where vibrations can have a wavevector
different from $\bm{q}{=}\bm{0}$.

For systems with low conductivity, which are the focus of this work, the WTE
can be solved accurately within the so-called single-mode relaxation-time
approximation (SMA, see Methods for details), which allows to write the
conductivity in the following compact form:

$\begin{split}\kappa{=}\frac{1}{\mathcal{V}{N_{\rm
c}}}{\sum_{\bm{q},s,s^{\prime}}}&\\!\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}}{4}\\!\left(\frac{C(\bm{q})_{s}}{\omega(\bm{q})_{s}}{+}\frac{C(\bm{q})_{s^{\prime}\\!}}{\omega(\bm{q})_{s^{\prime}\\!}}\right)\\!\frac{\rVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\lVert^{2}}{3}\\\
\times&\pi\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}]}(\omega(\bm{q})_{s}-\omega(\bm{q})_{s^{\prime}})\;,\end{split}$
(1)

where the wavevector $\bm{q}$ and the mode index $s$ label a vibrational
eigenstate having energy $\hbar\omega(\bm{q})_{s}$, anharmonic linewidth
$\hbar\Gamma(\bm{q})_{s}$, and specific heat

$C(\bm{q})_{s}{=}C[\omega(\bm{q})_{s}]{=}\frac{\hbar^{2}\omega^{2}(\bm{q})_{s}}{k_{\rm
B}T^{2}}{\mathsfit{N}}(\bm{q})_{s}\big{(}{\mathsfit{N}}(\bm{q})_{s}{+}1\big{)}$
(2)

(${\mathsfit{N}}(\bm{q})_{s}{=}[\exp(\hbar\omega(\bm{q})_{s}/k_{\rm
B}T){-}1]^{-1}$ is the Bose-Einstein distribution at temperature $T$); the
quantity

$\vspace*{-0.5mm}\rVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\lVert^{2}{=}\sum_{\alpha=1}^{3}\mathsfit{v}^{\alpha}(\bm{q})_{s,s^{\prime}}\mathsfit{v}^{\alpha}(\bm{q})_{s^{\prime},s}\vspace*{-0.5mm}$
(3)

denotes the square modulus of the velocity operator 37 between eigenstates $s$
and $s^{\prime}$ at the same wavevector $\bm{q}$ ($\alpha$ in these
expressions denotes a Cartesian direction, and since vitreous solids are in
general isotropic, the scalar conductivity (1) is computed as the average
trace of the tensor $\kappa^{\alpha\beta}$ given by Eq. (47) of Ref. 37;
$N_{\rm c}$ is the number of $\bm{q}$-points entering in such a summation and
$\mathcal{V}$ is the primitive-cell volume). Finally, $\mathcal{F}$ is a
Lorentzian distribution having a full width at half maximum (FWHM) equal to
$\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}$:

$\begin{split}&\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}]}\big{(}\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}}\big{)}\\\
&=\frac{1}{\pi}\frac{\frac{1}{2}\big{(}\Gamma(\bm{q})_{s}+\Gamma(\bm{q})_{s^{\prime}}\big{)}}{\big{(}\omega(\bm{q})_{s}-\omega(\bm{q})_{s^{\prime}}\big{)}^{2}+\frac{1}{4}\big{(}\Gamma(\bm{q})_{s}+\Gamma(\bm{q})_{s^{\prime}}\big{)}^{2}}\;.\end{split}$
(4)

In a crystal the primitive cell and the BZ have a finite volume and can be
univocally chosen relying on crystallographic conditions 57; wavevectors are
good quantum numbers and phonon group velocities are well defined. Under these
circumstances, it is useful to rewrite the WTE conductivity (1) as sum of two
terms, $\kappa{=}\kappa_{P}{+}\kappa_{C}$; specifically, the term $\kappa_{P}$
(referred to as “populations conductivity”37, 36) is determined by the
diagonal ($s{=}s^{\prime}$) or perfectly degenerate ($s{\neq}s^{\prime}$ with
$\omega(\bm{q})_{s}{=}\omega(\bm{q})_{s^{\prime}}$) terms in the summation in
the conductivity (1); in crystals such a term describes the Peierls-Boltzmann
particle-like heat conduction (averaged over the spatial directions), since
the average trace of the Peierls-Boltzmann conductivity tensor can be written
as $\kappa_{P}{=}\tfrac{1}{3}\sum_{\alpha}\kappa^{\alpha\alpha}_{P}$ with
$\kappa_{P}^{\alpha\alpha}{=}\frac{1}{\mathcal{V}N_{C}}\sum_{\bm{q}s}C[\omega(\bm{q})_{s}]{\mathsfit{v}^{\alpha}(\bm{q})_{s,s}}\Lambda^{\alpha}(\bm{q})_{s}$,
i.e. as particle-like vibrations having absolute energy
$\hbar\omega(\bm{q})_{s}$ (thus specific heat $C[\omega(\bm{q})_{s}]$) and
propagating between collisions over a length222Here we have exploited the
possibility to diagonalize at least one Cartesian components ${\alpha}$ of the
velocity operator in the degenerate subspace 108, 37, thus
${\mathsfit{v}^{\alpha}(\bm{q})_{s,s}}$ is the propagation velocity in
direction $\alpha$ and ${\mathsfit{v}^{\alpha}(\bm{q})_{s,s^{\prime}}}{=}0$
for $s{\neq}s^{\prime}$ in the degenerate subspace
$\Lambda^{\alpha}(\bm{q})_{s}{=}{\mathsfit{v}^{\alpha}(\bm{q})_{s,s}}[\Gamma(\bm{q})_{s}]^{-1}$.
Conversely, non-degenerate off-diagonal elements (“coherences”37, 36) account
for a different “Wigner” conduction mechanisms: they do not have an absolute
energy akin to that of a particle-like excitation, but are characterized by an
energy difference $\hbar\omega(\bm{q})_{s}{-}\hbar\omega(\bm{q})_{s^{\prime}}$
and describe a wave-like tunneling conduction mechanisms akin to the
electronic Zener interband tunnelling 59. It has been shown in Refs. 36, 37,
60 that in simple crystals particle-like mechanisms dominate and thus
$\kappa_{{}_{\rm P}}{\gg}\kappa_{{}_{\rm C}}$, while in complex crystals both
these mechanisms are relevant, and $\kappa_{{}_{\rm P}}$ and $\kappa_{{}_{\rm
C}}$ are of the same order. Finally, Refs. 36, 37 have shown that in the
ordered limit describing a harmonic glass8, 9, Eq. (1) reduces to the AF
formula for the conductivity of glasses. Specifically, such an ordered limit
requires first to describe a (structurally stable333With “structurally stable
glass” we mean a glass in which atoms vibrate around equilibrium positions,
thus atomic diffusion is negligible. See e.g. refs. 135, 136, 137, 138, 139,
140, 141, 142 for a discussion of how network structure and chemical
composition affect the structural stability of glasses and related phenomena.)
glass as the limiting case of a disordered but periodic crystal with an
increasingly larger primitive cell (i.e. $\mathcal{V}{\to}\infty$ and thus
with the BZ reducing to the point $\bm{q}{=}\bm{0}$ only 8, 9), and then
letting each linewidth go to the same infinitesimal broadening8, 9
$\hbar\eta$, $\hbar\Gamma(\bm{q})_{s}{\to}\hbar\eta{\to}0$, $\forall\;s$ and
$\bm{q}{=}\bm{0}$. Under these ideal circumstances only $\bm{q}{=}\bm{0}$ is
considered in the sum in Eq. (1), and the Lorentzian distribution (4) becomes
a Dirac $\delta$,

$\begin{split}\lim_{\eta{\to}0}\\!\\!\left[\lim_{\mathcal{V}\to\infty}\mathcal{F}_{[2\eta]}(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})\right]{=}\delta\left({\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}}}\right),\end{split}$
(5)

implying that the WTE conductivity (1) reduces exactly to the AF formula for
the conductivity of glasses (Eq. (3) of Ref.8). In practice, this ideal bulk-
glass limit cannot be reached in numerical calculations, and anharmonic
linewidths strongly vary with temperature. In the next sections we first
discuss a protocol that allows to accurately describe glasses at the AF
harmonic level using finite-size models, and then we rely on the WTE to extend
such a protocol accounting for anharmonicity.

##  III Describing glasses with finite-size models

###  III.1 Protocol to evaluate the harmonic Allen-Feldman conductivity

From a microscopic viewpoint, AF conductivity (Eq. (1) and (5) with
$\bm{q}=\bm{0}$ only) describes heat transport in an ideal bulk glass as being
mediated by a transfer of energy between vibrations that are degenerate in
energy (see Eq. (5)) and not localized in the Anderson sense62 (this
requirement is needed to have non-zero velocity operator elements in Eq. (1)).
In actual calculations a glass is approximately described as a crystal having
a primitive cell containing a large but finite number of atoms $N_{\rm at}$.
Such an approximation has two important implications. First, the BZ
corresponding to the (large) finite-size model does not reduce to
$\bm{q}{=}\bm{0}$ only but has a (small) finite volume. Second, in a realistic
finite-size model of a glass perfectly degenerate vibrational modes are
supposed to be absent (this because exact degeneracies are due to point-group
symmetries63, which are not present in amorphous systems); thus, the
vibrational spectrum is characterized by an average spacing between
vibrational energy levels equal to

$\hbar\Delta\omega_{\rm avg}=\frac{\hbar\omega_{\rm max}}{{3N_{\rm at}}},$ (6)

where $\hbar\omega_{\rm max}$ is the maximum vibrational energy of the solid
(which strongly depends on the chemical composition and negligibly depends on
disorder, see Fig. 10 in Methods) and $3N_{\rm at}$ is the number of
vibrational modes at a fixed $\bm{q}$ point. Eq. (6) underlines how degenerate
vibrational modes determining the harmonic conductivity emerge as a
consequence of disorder, since increasing the accuracy in the description of
disorder (i.e. increasing $N_{\rm at}$) yields a decrease of the average
energy-level spacing, implying that in the ideal limit of a bulk glass
($N_{\rm at}{\to}\infty$) adjacent vibrational eigenstates become degenerate
(hereafter we will use the term “quasi-degenerate” to refer to these adjacent
vibrational eigenstates that become degenerate in the ideal limit of a bulk
glass).

These two properties offer important insights on how to evaluate the strength
of finite-size effects in glasses and consequently extrapolate from the
finite-size “reference cell” of the model the behavior of the ideal (infinite)
glass. Specifically, the presence of a BZ with a small but finite volume
implies that Fourier interpolations can be used to sample vibrations in a
$n{\times}n{\times}n$ supercell of the finite periodic model. This procedure
simulates a system where vibrations are commensurate to a system that is
$n{\times}n{\times}n$ times larger than the reference cell used, but where the
disorder length scale remains limited by the size of the reference cell.
Within this scheme one can obtain information about finite-size effects in
multiple ways: (i) for a given finite model, one can study the differences
between a calculation performed at $\bm{q}{=}\bm{0}$ only, and a calculation
on a $n{\times}n{\times}n$ $\bm{q}$-mesh; (ii) one can repeat the analysis at
the previous point employing models having larger and larger reference cells —
for sufficiently large models one expects to achieve computational
convergence, with the calculation at $\bm{q}{=}\bm{0}$ only and that on the
$\bm{q}$-mesh yielding indistinguishable results. Moreover, in order to retain
the key physical property that couplings between quasi-degenerate eigenstates
can occur and contribute to heat transport, the Dirac $\delta$ (5) needs to be
replaced with a smooth distribution having a broadening $\eta$ of the order of
the average energy-level spacing8, 9, 34 $\hbar\Delta\omega_{\rm avg}$;
otherwise, no couplings would take place in any system represented by a
finite-size periodic supercell. Within such a numerical scheme $\eta$ is just
a computational parameter, and as such one expects to find a range of values
for $\eta$ (the domain of convergence) for which the conductivity is
independent from $\eta$. The considerations above show that the calculation of
the harmonic AF conductivity in finite-size models bears some procedural
analogies to electronic-structure calculations in metals, where the BZ is
integrated using a discrete mesh and the Dirac $\delta$ identifying the Fermi
level is broadened with an aptly chosen smearing function 64, 65, 66, and one
needs to identify a range of “converged values” for the smearing and for the
size of the BZ-integration mesh for which quantities such as the total energy
are practically independent from the smearing.

Pioneering work8, 9, 10, 34 evaluated the AF conductivity broadening the Dirac
$\delta$ (5) with a fat-tailed444We employ the nomenclature used in
statistics, which classifies as “light-tailed” (“fat-tailed”) all the
distributions having tails smaller (larger) than the exponential distribution.
Lorentzian distribution having FWHM of the order of $\hbar\Delta\omega_{\rm
avg}$, limited calculations to $\bm{q}{=}\bm{0}$ to reduce the computational
cost, and relied on smoothness arguments 9 to extrapolate the bulk limit from
finite-size models. Nevertheless, recent work 20, 33, 68 have highlighted
challenges in achieving computational convergence following such computational
protocol. Therefore, here we take inspiration from the computational
techniques alluded to above for electronic-structure calculations, to tackle
the problem of achieving convergence in the calculation of the AF
conductivity, using vitreous silica as a paradigmatic test case. We first show
in Fig. 1 that state-of-the-art models of vitreous silica49 ($v$-SiO2)
containing 192 or 5184 atoms (both studied using a recently developed GAP
potential49, see Methods for details) yield bulk AF conductivities that are
perfectly compatible, provided these are computed using a light-tailed
Gaussian broadening having a FWHM larger than the average energy-level spacing
(more details on the convergence domain for the broadening parameter are
discussed later), and a Fourier interpolation is used to extrapolate the bulk
limit from the small 192-atom model (the Fourier interpolation has a
negligible effect on the large 5184-atom model, confirming the aforementioned
expectations).

Figure 1:  Convergence of the Allen-Feldman theory in finite-size $v$-SiO2
models. Solid lines are the AF conductivities computed at $\bm{q}{=}\bm{0}$
only, empty circles and squares (crosses) are the AF conductivities computed
using the Fourier interpolation on a $3{\times}3{\times}3$
($5{\times}5{\times}5$) $\bm{q}$ mesh; results for a 192-atom (5184-atom)
model are in blue (red). All calculations are performed using the same
Gaussian broadening $\eta{=}4\;\rm{cm}^{-1}$ for the Dirac $\delta$ appearing
in the AF conductivity expression (such a value corresponds to a FWHM larger
than the average energy-level spacing of these models, and for this broadening
computational convergence is achieved in both structures, see text and Fig.
2). The discrepancy between the solid blue line and all the other data shows
that the Fourier interpolation is necessary to achieve computational
convergence using a small 192-atom model; the good agreement between scatter
points and the red line shows that computational convergence can be achieved
with a calculation at $\bm{q}{=}\bm{0}$ only for the larger 5184-atom model.
Figure 2: Convergence of the AF conductivity for $v-$SiO2 with respect to the
broadening $\eta$ for the Dirac $\delta$, for a 192-atom model (left) and for
a 5184-atom model (right). Dashed red lines represent a fat-tailed Lorentzian
having FWHM $2\eta$, solid black lines represent a light-tailed Gaussian with
variance $\eta^{2}{\pi/2}$ (i.e. the same maximum $(\pi\eta)^{-1}$ of the
Lorentzian with FWHM $2\eta$). The top row shows that evaluating the AF
conductivity using the $\bm{q}$-interpolation technique (see text) and the
Gaussian broadening yields a range of value for $\eta$ for which the AF
conductivity is not sensitive to the value of $\eta$ (i.e. a “convergence
plateau”) for both the 192- and 584-atom models. The horizontal orange and
blue lines show the computationally converged value for the AF conductivity at
100 and 300 K, respectively; we highlight that the 192- and 5184-atom models
yield compatible (indistinguishable) conductivities. In general, but
especially in the 192-atom model, the Gaussian broadening yields a wider and
more clear convergence plateau compared to the Lorentzian broadening; thus,
from a numerical viewpoint, the former has some advantage over the latter. The
vertical lines show that when the $\bm{q}$-interpolation technique is adopted,
a broadening $\eta=3\Delta\omega_{\rm avg}$ yields a computationally converged
thermal conductivity at both 100 and 300 K. In the bottom line we show that
the evaluation of the AF conductivity at $\bm{q}{=}\bm{0}$ only yields in the
large 5184-atom model a convergence plateau that is smaller compared to that
obtained in the same structure using $\bm{q}$-interpolation; at
$\bm{q}{=}\bm{0}$ the small 192-atom model instead does not show clear
convergence with respect to the broadening. The green arrows show the value of
$\eta=4$ cm-1, which yields a computationally converged conductivity in the
192-atom model in the $\bm{q}$-mesh calculation, as well as in 5184-atom model
in both the $\bm{q}$-mesh and $\bm{q}{=}\bm{0}$ only calculations. The origin
of the opposite trend of the broadening-conductivity curve obtained using
$\bm{q}$ interpolation or $\bm{q}{=}\bm{0}$ only is discussed in Sec. III.2.

The details on the domain of convergence for the broadening parameter $\eta$
and mesh used for the BZ sampling are reported in Fig. 2, where we show that
computational convergence can be achieved for both the 192- and 5184-atom
models, i.e. there exists a range of values for $\eta$ over which the
conductivity is not sensitive to the value of $\eta$ (“convergence plateau”).
In the larger 5184-atom model the convergence plateau extends to smaller
values of $\eta$; this in agreement with the expectation, based on Eq. (6),
that larger models allow for more accurate approximations of the Dirac
$\delta$. Importantly, we also show that the Gaussian broadening yields a
wider convergence plateau compared to the Lorentzian broadening, especially at
low temperatures where anharmonicity phases out and thus the harmonic AF
theory is accurate. We note in passing that the improved computational
performance found here for the Gaussian broadening compared to the Lorentzian
broadening bears analogies to density-functional-theory calculations for
metals, where refined representations of the Dirac delta 69, 66, 64, 65, 66
are used in place of the Fermi-Dirac thermal broadening 70 to improve
convergence with respect to Brillouin-zone sampling. The analysis of Fig. 2
demonstrates that the evaluation of the AF conductivity at $\bm{q}{=}\bm{0}$
only does not show clear convergence with respect to the broadening for the
192-atom model; a convergence plateau emerges instead if the Fourier
interpolation is employed to extrapolate to the bulk limit. In contrast, in
the 5184-atom model computing the AF conductivity at $\bm{q}=\bm{0}$ only or
using an interpolation mesh does not produce significant differences, provided
a value of $\eta$ belonging to the convergence plateau is used. This shows
that the protocol of using the $\bm{q}$-mesh interpolation and (ideally)
Gaussian broadening to compute the AF conductivity allows to achieve
computational convergence in three small models (containing $\lesssim 200$
atoms) typically affordable in first-principles studies; in $v$-SiO2 the model
containing 5184 atoms already at $\bm{q}=\bm{0}$ is accurate, and allows to
determine a convergence plateau for the broadening $\eta$.

###  III.2 Extension of the protocol to evaluate the anharmonic Wigner
conductivity

Figure 3:  Vibrational frequencies and anharmonic linewidths of vitreous
silica. Panel a): the scatter points represents the linewidths as a function
of frequency and temperature for the 192-atom model of $v-$SiO2; these are
computed at $\bm{q}{=}\bm{0}$ only accounting for third-order anharmonicity71,
72, 73, 74, 75 and scattering due to isotopic mass disorder at natural
abundance76. The solid lines represent a coarse graining of the frequency-
linewidth distribution into a single-valued function77, 78
$\Gamma_{a}[\omega]$ (see Methods for details). The horizontal black (gray)
lines show the average energy-level spacing for the 192-atom (5184-atom)
model; for the 192-atom (5184-atom) model all the vibrations below the black
(gray) line are not accurately accounted for by the bare WTE (Eq. (1) and
(4)), and need to be regularized using the protocol described in the main
text. Panel b): the solid red line represents the vibrational density of state
(vDOS) computed using the 5184-atom model and considering $\bm{q}{=}\bm{0}$
only, the overlapping dashed orange line is the vDOS computed using the
192-atom model over a $3{\times}3{\times}3$ $\bm{q}$-mesh (both these
calculations use the same Gaussian broadening $\eta=4\rm{cm}^{-1}$, discussed
also in Figs. 1,2). The vertical lines show how the vibrations of $v-$SiO2 are
sampled using the 192-atom model at $\bm{q}{=}\bm{0}$ only (top, corresponding
to the abscissas of the scatter points in panel a)), using the 192-atom model
and relying on Fourier interpolation on a $3{\times}3{\times}3$ $\bm{q}$-mesh
(center), and using the 5184-atom model at $\bm{q}{=}\bm{0}$ only (bottom).

We start by recalling that the WTE also generalizes the AF model accounting
for anharmonicity; in particular, the Lorentzian distribution (4) appearing in
the WTE conductivity (1) has a FWHM determined by the anharmonic linewidths.
From a computational perspective, and recalling that the size of the model
determines the average energy-level spacing (6), one expects that evaluating
the WTE using a finite-size model should lead to results negligibly affected
by finite-size effects whenever all the anharmonic linewidths of the finite-
size model are larger than its average energy-level spacing. To proceed, we
first show in Fig. 3a) the frequency-linewidth distributions for the 192- and
5184-atom $v-$SiO2 models, and compare these with the average energy-level
spacing. Each cloud of points represents the frequency-linewidth distribution
at a given temperature for the 192-atom model, evaluated explicitly at
$\bm{q}=\bm{0}$ only and accounting both for third-order anharmonicity71, 72,
73, 74, 75 and for scattering due to isotopic mass disorder at natural
abundance76 (see Methods for details). Each solid line is a coarse-grained
interpolation of the frequency-linewidth distribution into a single-valued
function $\Gamma_{a}[\omega]$; such a coarse graining is inspired by past work
77, 78 and validated in Fig. 11 in the Methods, where we show that evaluating
the WTE conductivity using the exact frequency-linewidth distribution or
linewidths determined using the coarse grained function
$\Gamma(\bm{q})_{s}{=}\Gamma_{a}[\omega(\bm{q})_{s}]$ yields practically
indistinguishable results. Clearly, in the 192-atom model a significant
portion of the vibrational modes have linewidths below the average energy-
level spacing at room temperature and below (this happens both when vibrations
are sampled at $\bm{q}{=}\bm{0}$ only, as well as when they are more densely
sampled on a $3\\!{\times}\\!3\\!{\times}\\!3$ $\bm{q}$-mesh). For the
5184-atom model, instead, the average energy-level spacing is more than one
order of magnitude smaller than that of the 192-atom model, implying that
already at 50 K most of the linewidths are larger than the average energy-
level spacing.

To better understand the combined effect of anharmonicity and of the finite
size of the model on conductivity, we recall once again that when the WTE
conductivity (1) and (4) is evaluated repeating periodically a finite-size
reference cell, the conductivity can be decomposed in the sum of a populations
and a coherences contribution ($\kappa_{\rm P}$ and $\kappa_{\rm C}$,
respectively, see Sec. II). We stress that the populations (coherences)
conductivity can be rationalized in terms of a particle-like (wave-like)
conduction mechanisms exclusively in actual crystals, where the periodicity
and symmetries of the atomic structure imply that the BZ and the vibrational
spectrum can be defined unambiguously57; consequently, the group velocities
with which the particle-like phonon wavepackets propagate, and the interband
spacings that characterize the wave-like tunneling of phonons, are well
defined. In finite-size glass models, instead, the decomposition in
populations and coherences conductivities is useful to understand finite-size
effect (more on this later), but cannot be further interpreted in terms of
particle-like and wave-like conduction mechanisms. In fact, such an
interpretation is well defined only when the BZ and vibrational spectrum can
be unambiguously defined, while in finite-size models of glasses the BZ volume
decreases as the size of the model increases. Finally, we also note that the
distinction between populations and coherences conductivities in finite-size
models of glasses is different from the distinction between “propagons” and
“diffusons” conductivities discussed by Allen et al.10 for an ideal
(infinitely large) glass. The propagons and diffusons conductivities
differentiate the contributions to heat transport originating from vibrational
modes that propagate ballistically (akin to phonon wavepackets) from those
that do not in an ideal glass; future work will aim at decomposing the total
Wigner conductivity into those contributions.

In finite-size models of glasses the populations conductivity $\kappa_{\rm P}$
is given by the terms diagonal in the mode index ($s{=}s^{\prime}$ in the sum
appearing in Eq. (1)), since perfectly degenerate vibrational modes are absent
(or a negligible fraction of the total number of modes, since degeneracies are
due to point-group symmetries63 and these are supposed to be absent in a
realistic finite-size model of a glass). For an ideal, infinitely large model
of a glass, which can be described accurately by a calculation at
$\bm{q}{=}\bm{0}$ only, one would obtain $\kappa_{\rm P}{=}0$, since at
$\bm{q}{=}\bm{0}$ acoustic vibrations have zero specific heat and the time-
reversal symmetry of the dynamical matrix63 implies that the group velocities
for optical vibrations are zero; the only non-zero contribution to the total
conductivity would be the coherences conductivity $\kappa_{\rm C}$, provided
some non-zero off-diagonal velocity-operator elements exists (a condition that
is verified when vibrations are not Anderson-localized 9, 10).

Figure 4:  Convergence of the anharmonic Wigner conductivity in finite-size
$v$-SiO2 models. Top, bare WTE conductivity (Eq. (1) and (4)) obtained using
the 192-atom (blue) or 5184-atom (orange) model, and evaluated at
$\bm{q}{=}\bm{0}$ (solid lines) or on a $3{\times}3{\times}3$ $\bm{q}$ mesh to
extrapolate to the bulk limit (dashed lines). Importantly, the crystal-like
divergence of the WTE conductivity computed using the $\bm{q}$-mesh is a
finite-size effect, which occurs at lower temperatures for larger models. For
$T{>}250$ K, in the 5184-atom model the conductivity computed at
$\bm{q}{=}\bm{0}$ is practically equal to the fully $\bm{q}$-sampled
conductivity; for the 192-atom model, instead, the conductivity computed at
$\bm{q}{=}\bm{0}$ underestimates that limit. Importantly, for $T{>}250$ K the
fully-sampled conductivity obtained from the 192-atom model is compatible with
that obtained from the 5184-atom model. Bottom panel, regularized WTE
conductivity (rWTE, see text for a discussion of the regularization) for the
192- (solid black) or 5184-atom (empty diamonds) model, computed in both cases
using a $3{\times}3{\times}3$ $\bm{q}$-mesh, corresponding to the full
$\bm{q}$-sampling limit. The rWTE smoothly connects the correct low- and high-
temperature limits, i.e. the fully $\bm{q}$-sampled harmonic AF conductivity
(as in Fig. 1, the dashed-dotted red line is the 192-atom model, empty red
circles are the 5184-atom model) and the bare WTE (dashed blue, reported from
top panel), respectively. Inset: the bare WTE conductivity evaluated at
$\bm{q}{=}\bm{0}$ for the 5184-atom model is practically equal to the fully
$\bm{q}$-sampled rWTE conductivities of the 192- or 5184-atom models. This
demonstrates that the full $\bm{q}$-sampling limit corresponds to the bulk
limit.

As discussed in the previous section, a technique to extrapolate the bulk
limit from finite-size models consists in relying on Fourier interpolation to
sample vibrations in a $n{\times}n{\times}n$ Born von-Karman supercell of the
finite periodic reference cell, while remaining aware of the limitations
stemming from having a disorder length scale limited by the size of the
reference cell. Fig. 3b) shows that such Fourier interpolation allows to
greatly improve the accuracy of the thermodynamic predictions for the 192-atom
model; specifically, the vibrational density of states (vDOS) of the 192-atom
model computed using a $3{\times}3{\times}3$ $\bm{q}$-mesh is in very good
agreement with the vDOS of the 5184-atom model computed at $\bm{q}{=}\bm{0}$
only. For reference cells larger than a certain size, one expects the
conductivities computed at $\bm{q}{=}\bm{0}$ only or using a
$n{\times}n{\times}n$ Fourier interpolation to be practically
indistinguishable555To support this expectation, we note that the sum over the
modes appearing in Eq. (1) contains $(3N_{at})^{2}$ term, out of which only
$3N_{at}$ are diagonal ($s=s^{\prime}$) and thus vanishing at
$\bm{q}{=}\bm{0}$. This implies the difference between a conductivity
calculation at $\bm{q}{=}\bm{0}$ only and one on the $n{\times}n{\times}n$
mesh goes to zero with a speed roughly proportional to the ratio between the
number of diagonal elements and total number of elements, i.e.
$\frac{3N_{at}}{(3N_{at})^{2}}=\frac{1}{3N_{at}}$.. In practice, to test the
accuracy of the bulk-limit extrapolations performed using the $\bm{q}$
interpolation, and to verify the aforementioned expectations, one has to
perform calculations in models having different sizes (ideally in one large
model that already describes the bulk limit and for which the $\bm{q}$
interpolation is therefore not needed, and in one small model where using or
not the interpolation technique is expected to yield appreciable effects). In
Fig. 4 (upper panel) we compare the WTE conductivity (Eq. (1) and (4))
obtained using the 192-atom model (blue) or the 5184-atom model (orange)
calculated either using $\bm{q}{=}\bm{0}$ only (solid lines) or using a
$3{\times}3{\times}3$ Fourier interpolation (dashed lines). These two
calculations differ only in the low-temperature limit, where they yield
conductivities having opposite trends. The divergence of the WTE conductivity
computed using the Fourier interpolation is a finite-size effect arising from
the periodic boundary conditions (a reminiscence of the divergence at low
temperatures of the conductivity of bulk crystals that is cutoff in real
crystals by the scattering with grains’ or samples’ boundaries); such an
effect occurs at lower temperatures for larger modes; thus, it is expected to
vanish in the ideal glass limit, i.e. for $N_{\rm at}{\to}\infty$ $\kappa(T)$
in the above-the-plateau regime is expected to follow the same trend observed
in the calculation at $\bm{q}{=}\bm{0}$ and in experiments (more on this
later). We note, in passing, that this reasoning also explains the opposite
trend of the broadening AF conductivity curve computed at $\bm{q}{=}\bm{0}$
only or using the Fourier interpolation shown in Fig. 2. In fact, in the limit
of vanishing broadening $\eta{\to}0$ only the term $s{=}s^{\prime}$ determines
the value of the sum in Eq. (1), yielding
$\kappa{\approx}\frac{1}{\mathcal{V}N_{C}}\sum_{\bm{q}s}C[\omega(\bm{q})_{s}]\frac{|\\!|{\mathsfit{v}(\bm{q})_{s,s}}|\\!|^{2}}{\eta}$.
Such a limiting expression for the conductivity diverges when the Fourier
interpolation is adopted, since at $\bm{q}{\neq}\bm{0}$ the diagonal velocity-
operator elements are non-zero; in contrast, in a calculation at
$\bm{q}{=}\bm{0}$ only the time-reversal symmetry implies that the diagonal
velocity-operator elements are zero; thus, in this case the conductivity
approaches zero when the broadening approaches zero.

We note that for the 5184-atom model at temperatures higher than 100 K the WTE
conductivity obtained relying on the Fourier interpolation is practically
equal to the conductivity obtained evaluating the WTE conductivity at
$\bm{q}{=}\bm{0}$ only. In contrast, for the 192-atom model the conductivity
obtained from a calculation at $\bm{q}{=}\bm{0}$ only is always significantly
different from that obtained using a computationally converged
$3{\times}3{\times}3$ $\bm{q}$ mesh666Using a denser $5{\times}5{\times}5$
mesh produces practically equivalent results, implying that in the 192-atom
model a $3{\times}3{\times}3$ mesh is sufficient to obtain computational
convergence.. Most importantly, for temperatures higher than 250 K the WTE
conductivity obtained from the 192-atom model using the Fourier interpolation
is practically indistinguishable from that obtained from the 5184-atom model,
confirming that the interpolation allows to accelerate convergence in the
calculation of the bulk limit.

The considerations above show that in finite-size models the WTE conductivity
with Fourier interpolation is accurate up to a lower-bound temperature $T_{L}$
(roughly defined as the temperature at which the temperature-conductivity
curve changes concavity). Below $T_{L}$ finite-size effects lead to a crystal-
like divergence, which emerges as a consequence of having a significant number
of vibrational modes with anharmonic linewidths $\hbar\Gamma(\bm{q})_{s}$
smaller than $\hbar\Delta\omega_{\rm avg}$ (see Fig. 3); from a microscopic
viewpoint, this implies that couplings in the distribution (4) between quasi-
degenerate eigenstates — which are present in an ideal glass, see Sec. III.1 —
are not correctly accounted for. Thus, the finite-size model fails to
represent the harmonic conduction described by the distribution (5) and
accurate for an ideal glass at very low temperature777here “very low
temperature” must be interpreted keeping into account that this work is
focused on the above-the-plateau temperature regime, i.e. $T>30K$.. This
limitation can be overcome relying on the fact that in the low-temperature
limit anharmonicity progressively phases out; thus, the AF model becomes
increasingly more accurate and can be evaluated using the protocol discussed
in the previous section. Therefore, we introduce a regularization protocol for
the WTE that allows to determine its bulk limit using Fourier interpolation,
and accounts for the prescriptions needed to correctly evaluate the low-
temperature harmonic limit discussed before. Specifically, we choose for our
protocol a Voigt profile 82 — a two-parameter distribution
$\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}},\eta]}$
obtained as a convolution between a Lorentzian with FWHM
$\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}$ and a Gaussian with
variance $\eta^{2}\pi/2$ — in place of the one-parameter distribution
$\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}]}$ appearing
in Eq. (1). By doing so Eq. (1) reduces to the AF harmonic limit at low
temperatures (where $\Delta\omega_{\rm
avg}{\sim}\eta{\gg}\frac{1}{2}(\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}})$),
and to the anharmonic WTE at intermediate and high temperatures (where
$\Delta\omega_{\rm
avg}{\sim}\eta{\ll}\frac{1}{2}(\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}})$).
Hereafter the conductivity computed using the Fourier interpolation (to
determine the bulk limit) and the Voigt distribution (to correctly describe
the low-temperature harmonic limit) will be referred to as “regularized WTE”
(rWTE), to distinguish it from the “bare” WTE (Eq. (1) with the Lorentzian
distribution (4)). As shown in the bottom panel of Fig. 4, the rWTE
conductivity reduces to the AF harmonic conductivity for temperatures lower
than $T_{L}$, and to the bare anharmonic WTE conductivity for temperature
higher than $T_{L}$, smoothly connecting these two limits. The inset of Fig. 4
highlights that the rWTE conductivity evaluated on a $3{\times}3{\times}3$
$\bm{q}$-mesh for the 192-atom model is practically indistinguishable from the
rWTE conductivity of the 5184-atom model (we note in passing that for the
5184-atom model the WTE conductivity evaluated at $\bm{q}{=}\bm{0}$ only is
equal to the rWTE conductivity, showing that the 5184-atom model is large
enough to describe the bulk limit without relying on the $\bm{q}$
interpolation and on the regularization). Therefore, the regularization
protocol greatly accelerates the convergence of the conductivity calculation,
allowing to determine the bulk limit of the conductivity of $v-$SiO2 using
models containing less than 200 atoms.

##  IV First-principles calculations and comparison with experiments

Figure 5:  Thermal conductivity of $v$-SiO2: theory vs experiments. Panel a),
solid (dotted) lines are the rWTE (AF) conductivities for a 192-atom model of
$v$-SiO2 generated using first-principles techniques54 (denoted with “192(D)”
in the text), repeated periodically $3{\times}3{\times}3$ times and then
studied using directly first-principles DFT calculations (red), a GAP
potential (green), or the semi-empirical BKS potential (purple). The dashed
lines are the classical AF conductivities (represented using the same color
code as before), corresponding to atomic vibrations having specific heat
determined by equipartition and equal to $k_{B}$ for all modes 35. Cyan and
purple scatter points are predictions from first-principles and classical
(BKS-based) molecular dynamics simulations, respectively, from Ercole et
al.27. These MD calculations are governed by classical equipartition 35, they
generally agree with the corresponding classical AF conductivities
(increasingly larger differences between MD and classical AF observed lowering
temperature might be due to the decrease of the accuracy with which atomic
vibrations are sampled in MD simulations when temperature is lowered). We also
note that in the high-temperature limit, when the correct quantum Bose-
Einstein statistics of vibrations yields a specific heat approaching the
classical frequency-independent limit, the MD calculations are in good
agreement with the corresponding (first-principles- or BKS-based) rWTE
predictions. Panel b), gray scatter points are experiments from Wray et al.
83, Kanamori et al. 84, Touloukian et al. 85, Sergeev et al. 86, Cahill 87,
Lee and Cahill (190 nm sample) 88, Abdulagatov et al. 89, Yang et al. 90, and
Heraeus 91. The red line represents the (most accurate) rWTE conductivity
obtained from first principles, which is in agreement with experiments in the
temperature range where radiative effects are negligible ($50\lesssim
T\lesssim 450$ K). Above 500 K, our calculations depart from the increasing
trend observed in experiments, showing that additional theoretical work is
needed to analyze the regime where radiative effects may be relevant.

Hitherto, no theoretical work has managed to evaluate the thermal conductivity
of $v$-SiO2 (and more generally of any amorphous solid) from first principles
and accounting for the interplay between anharmonicity, disorder, and the
Bose-Einstein statistics of atomic vibrations. More precisely, past works
studied the thermal conductivity of $v$-SiO2 using a variety of approaches,
including: (i) the AF model in combination with empirical 92 or semi-empirical
20 potentials; (ii) classical molecular dynamics 29, 20, 30, 23; (iii) first-
principles molecular dynamics 27. We note that the first-principles GKMD
simulation by Ercole et. al.27 represents a significant step forward in the
field — it is the first study that evaluates the conductivity of $v$-SiO2 from
first principles — with Ref.27 focusing on the high-temperature regime, where
the difference between the actual quantum Bose-Einstein occupation numbers of
vibrations and the classical (equipartition-determined) occupation numbers
implicit in the GKMD simulations35 is minimal. We have seen the rWTE protocol
discussed in the previous section allows to accurately determine the bulk
thermal conductivity of glasses using models with a size ($\lesssim$ 200
atoms) that is within the reach of first-principles techniques; so we now
employ these to compute the conductivity of $v$-SiO2.

We show in Fig. 5a) the bulk rWTE conductivity (solid lines) of a “192(D)”
model of $v$-SiO2, which contains 192 atoms and was generated relying on
density-functional theory (DFT) 54, 51 (this model is different from the
192-atom model generated with GAP and discussed in Sec. III; to avoid
confusion we will henceforth denote the 192-atom model discussed in Sec. III
with “192(G)”). To assess the effects of anharmonicity on the conductivity,
also the harmonic AF conductivity is shown (dotted lines). All the parameters
entering in the rWTE or AF expressions have been evaluated either from first-
principles (red, see Methods for details), or using a state-of-the-art GAP
potential49 (green), or using the well known semi-empirical BKS potential55,
56 (purple)888In all these cases we used a broadening $\eta=4$ cm-1 for the
Voigt distribution, this value was determined from a convergence test
analogous to that reported in the upper-left panel of Fig. 2.. Orange and
purple scatter points are results from first-principles999Obtained using the
PBE exchange-correlation functional. GKMD and classical GKMD using the BKS
potential, respectively, by Ercole et al. 27; these are meant to be accurate
at high temperature where the quantum Bose-Einstein vibration occupation
numbers approach the classical (equipartition) occupations underlying the GKMD
simulations 35. Importantly, Fig. 5a) shows that the conductivity of $v$-SiO2
is negligibly affected by anharmonicity. In fact, the rWTE perturbatively
accounts for anharmonicity at the lowest cubic order, and yields a thermal
conductivity that is (i) practically indistinguishable from that obtained
employing the harmonic Allen-Feldman theory over the entire temperature range
analyzed (30-1300 K); (ii) in good agreement at high temperature
($T{\gtrsim}500$ K) with the conductivity obtained by Ercole and Baroni using
GKMD 27, which accounts for anharmonicity exactly (in the high-temperature
limit our rWTE calculations, either from first-principles or based on the BKS
potential, are in very good agreement with the corresponding GKMD
predictions). Fig. 5a) reinforces the notion that the quantum Bose-Einstein
statistics of vibrations plays a crucial role in determining the thermal
conductivity 35 at low temperatures, since the conductivities computed
accounting for the quantum statistics (and for the correct bulk limit)
increase up to saturation with temperature, while the GKMD conductivities are
governed by classical equipartition and are roughly independent from
temperature. In addition, we highlight how studying the same 192(D) structure
using first-principles techniques (red) or the GAP potential (green) yields
very similar results, further endorsing the notion that machine-learned GAP
potentials can have an accuracy comparable to that of first-principles
techniques95, 96. In contrast, using the semi-empirical BKS potential yields
results that are significantly different from those obtained from first-
principles.

Figure 6: rWTE conductivity of various $v$-SiO2 models. Solid red is the
first-principles rWTE conductivity of the 192-atom model generated relying on
first-principles calculations54, computed using a $3{\times}3{\times}3$
$\bm{q}$ mesh; the empty diamonds show that using a denser
$5{\times}5{\times}5$ $\bm{q}$ mesh does not yield significant changes. The
solid-gray and dashed-blue lines show the first-principles rWTE conductivity
computed using a $5{\times}5{\times}5$ $\bm{q}$ mesh for a 108-51 and
144-atom52 model, respectively; both models have been generated relying on
first-principles calculations. For completeness, the conductivities of the
192- (black line) and 5184-atom (dashed orange) models discussed in Fig. 4 are
reported (these are labeled with the suffix “(G)” to distinguish them from the
structures generated relying on first-principles techniques, which are labeled
with “(D)”); the horizontal dashed-dotted lines highlight the $9\%$ difference
between the rWTE conductivities of the 192(D) and 192(G) models, both
described with GAP. All these atomistic models are available on the
MaterialsCloud archive97.

Importantly, Fig. 5b) shows that the first-principles rWTE conductivity is in
good agreement with experiments in the temperature range from 50 to $\sim$450
K; at higher temperatures the rWTE predicts a saturating trend for the
conductivity (which mirrors the saturating trend of the specific heat, see
Fig. 12 in the Methods section), while experiments display an increasing
drift. Past experimental work (e.g. Ref. 98 and references therein) hinted
that radiative energy transfer causes an increase of the thermal conductivity
in vitreous silica at high temperature, and highlighted how accurately
distinguishing the conduction and radiation contributions to heat transfer is
particularly challenging90 (among the various sources of difficulty, Refs. 84,
99 mentioned how the interplay between conduction and radiation depends from
the size and shape of the sample used in the experiment). Thus, in Fig. 5 the
discrepancy between theory and experiments at high temperature might be due to
having non-negligible radiative contributions in the experiments 90, 84, and
to not accounting for these radiative contributions in our calculations.

In Fig. 5, the bulk limit has been computed using a $3{\times}3{\times}3$
Fourier-interpolation mesh (thus corresponding to vibrations in a system
containing $192{\cdot}3^{3}{=}5184$ atoms); Fig. 6 demonstrates that using
this sampling achieves computational convergence, since increasing the mesh to
$5{\times}5{\times}5$ yields practically indistinguishable results (empty red
diamonds refer to a $5{\times}5{\times}5$ mesh, solid red line refers to a
$3{\times}3{\times}3$ mesh). Fig. 6 also shows that studying both the 192(D)
and the 192(G) models with GAP gives high-temperature limits for their rWTE
conductivities differing by about 9$\%$. Such a small difference is
particularly reassuring, given the two different methods employed to produce
these models (bond switching in the first case54, 51 and melt-quench in the
second case49). Moreover, Fig. 6 shows that the rWTE conductivity computed
from first principles using the 192(D) model (solid red) is very similar to
the conductivity computed from first-principles density-functional theory
using models containing 108 51 (solid gray) or 144 atoms 52, 53 (dashed blue),
both generated relying on first-principles calculations (therefore, we
henceforth refer to these models to as “108(D)” and “144(D)”, respectively).

##  V Velocity operator, anharmonicity, and trend of $\kappa(T)$

Figure 7: Velocity operator of vitreous silica and conductivity saturation
with temperature. Panel a), average square modulus of the velocity-operator
elements $\big{<}|\mathsfit{v}^{\rm avg}_{\omega_{\rm a}\omega_{\rm
d}}|^{2}\big{>}$ for the 192(D) model of $v$-Si$O_{2}$, computed from first
principles and represented as a function of the energy differences
($\hbar\omega_{\rm
d}{=}\hbar(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})$) and averages
($\hbar\omega_{\rm
a}{=}\hbar\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}}{2}$) of the
two coupled eigenstates (having wavevector $\bm{q}$ and modes $s$,
$s^{\prime}$; see text for details). The one-dimensional projections in panel
b) show that the elements $\big{<}|\mathsfit{v}^{\rm avg}_{\omega_{\rm
a}\omega_{\rm d}}|^{2}\big{>}$ are almost unchanged at a given average for
increasingly large energy differences. For increasingly larger temperatures,
these almost-constant elements drive the saturation the rWTE conductivity (Eq.
(1) with the Voigt distribution discussed in Sec. III.2), yielding results
very close to the Allen-Feldman conductivity curve (Fig. 5).

Now we want to rely on the results obtained for $v$-SiO2 to gain general
insight on how anharmonicity affects the high-temperature trend for the
$\kappa(T)$ curve in glasses. We start by recalling that, according to the
expression for the rWTE conductivity (Eq. (1) with the Voigt distribution
discussed in Sec. III.2), at high temperatures heat conduction in glasses is
mainly determined by couplings between vibrational eigenstates, and the
strengths of these couplings are given by the average square modulus of the
velocity-operator elements
$\frac{1}{3}\lVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\rVert^{2}$. More
precisely, in the high-temperature limit the specific heat is close to the
classical limit and thus can be considered independent from temperature, and
the linewidths are generally larger than the average energy-level spacing
(Fig. 3); therefore, the rWTE conductivity reduces to a Lorentzian-weighted
average of velocity-operator elements. The weights in such average depend on
the energy difference between the eigenstates coupled, and the temperature-
dependent anharmonic linewidths are the scale parameters (FWHM) of the
Lorentzian distributions (Eq. (4)). Increasing temperature yields larger
linewidths (Fig. 3), resulting in a broader Lorentzian distribution (4) that,
when determining the conductivity, gives more weight to velocity-operator
elements with larger energy differences. Therefore, the variation of the
elements of the velocity operator with respect to the energy differences
$\hbar\omega_{\rm
d}{=}\hbar(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})$ determines how
the conductivity varies with increasing temperature101010In this work the
velocity operator is considered independent from anharmonicity and
temperature, since we employ the standard approximation of considering the
force constants to be independent from anharmonicity and temperature (for
vitreous silica this approximation is realistic, as discussed in the
Methods).: matrix elements increasing or decreasing with $\omega_{d}$ imply a
conductivity increasing or decreasing with temperature. From this reasoning
and from the saturating trend of the temperature-conductivity curve for
vitreous silica shown in Fig. 5, we expect the velocity-operator elements for
$v$-SiO2 to be approximately constant with respect to $\omega_{d}$. To verify
this prediction, we plot in Fig. 7a) the velocity-operator elements for
$v$-SiO2 (we show the velocity operator for the 192(D) model and computed from
first principles; the other models yield practically indistinguishable results
when analyzed from first principles or using GAP) as a function of the energy
difference $\hbar\omega_{\rm d}$ and of the energy average $\hbar\omega_{\rm
a}{=}\hbar({\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}})/{2}$:

$\begin{split}&\big{<}|\mathsfit{v}^{\rm avg}_{\omega_{\rm a}\omega_{\rm
d}}|^{2}\big{>}=[\mathcal{G}(\omega_{\rm a},\omega_{\rm
d})]^{-1}\frac{1}{\mathcal{V}N_{\rm
c}}{\sum_{\bm{q},s,s^{\prime}}}\frac{\rVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\rVert^{2}}{3}\\\
&\hskip 28.45274pt{\times}\delta\left(\omega_{\rm
d}{-}(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})\right)\delta\Big{(}\omega_{\rm
a}{-}\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}}{2}\Big{)}\;;\end{split}$
(7)

$\mathcal{G}(\omega_{\rm a},\omega_{\rm d})$ is a density of states that
serves as normalization (see Methods for details). The one-dimensional
projections in Fig. 7b) show that these velocity-operator elements are almost
constant at varying $\hbar\omega_{\rm d}$ for all values of $\hbar\omega_{\rm
a}$. These findings validate the reasoning above: the variation of the
velocity-operator elements
$\frac{1}{3}\lVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\rVert^{2}$ with
respect to the energy differences $\hbar\omega_{\rm d}$ determines how the
conductivity varies with increasing temperature, implying that a saturating
temperature-conductivity curve is obtained when most of the velocity-operator
elements do not vary with $\hbar\omega_{\rm d}$. The reasoning above also
explains the small or negligible effects of anharmonicity on the conductivity
of $v$-SiO2 discussed before, since having a velocity operator displaying a
negligible dependence on $\hbar\omega_{\rm d}$ implies that the rWTE
conductivity does not vary appreciably when the linewidths (broadening of the
Lorentzian distribution) vary, and therefore it does not significantly differ
from the AF conductivity (computed using a constant broadening $\eta$).

##  VI Thermal diffusivity

In order to gain further insight on the microscopic mechanisms underlying
conduction it is useful to resolve how each vibrational mode contributes to
transport; i.e., the quantity of heat that it carries and the rate at which it
diffuses. It is possible to extract the contribution of a single vibration to
thermal transport by factorizing the single-vibration specific heat
$C(\bm{q})_{s}$ in the regularized rWTE conductivity (1) recasting it as
$\kappa=\frac{1}{\mathcal{V}N_{\rm
c}}\sum_{\bm{q}s}C(\bm{q})_{s}D(\bm{q})_{s}$, with $D(\bm{q})_{s}$ being the
“anharmonic thermal diffusivity”. The expression for $D(\bm{q})_{s}$ is
determined by such factorization and by the requirement that in the
coherences’ coupling between two vibrations $(\bm{q})_{s}$ and
$(\bm{q})_{s^{\prime}}$ each contributes to the coupling with a weight equal
to the relative specific heat 37 (e.g. for vibration $(\bm{q})_{s}$ the weight
is $\tfrac{C(\bm{q})_{s}}{C(\bm{q})_{s}+C(\bm{q})_{s^{\prime}}}$, and
correspondingly for vibration $(\bm{q})_{s^{\prime}}$ the weight is
$\tfrac{C(\bm{q})_{s^{\prime}}}{C(\bm{q})_{s}+C(\bm{q})_{s^{\prime}}}$):

$\begin{split}D(\bm{q})_{s}=&{\sum_{s^{\prime}}}\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}\\!}}{{2[C(\bm{q})_{s}{+}C(\bm{q})_{s^{\prime}\\!}]}}\\!\left[\\!\frac{C(\bm{q})_{s}}{\omega(\bm{q})_{s}}{+}\frac{C(\bm{q})_{s^{\prime}}}{\omega(\bm{q})_{s^{\prime}}}\\!\right]\\!\\!\frac{\lVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\lVert^{2}\\!}{3}\\\
&\times\pi\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}},\eta]}(\omega(\bm{q})_{s}-\omega(\bm{q})_{s^{\prime}})\;.\end{split}$
(8)

The goal of this decomposition is to resolve the rate at which the heat
carried by a vibration with wavevector $\bm{q}$ and mode $s$ diffuses. We note
that Eq. (8) accounts for the effects of anharmonicity on vibrations’
diffusion by means of the linewidths, depends on temperature through both the
specific heat and the linewidths, and applies to both glasses and crystals —
in the former case the wavevector $\bm{q}$ is just a label without direct
physical meaning (we will discuss later that for glasses the diffusivity has
to be represented as a function of frequency to be well defined), while in the
latter case such an expression is accurate if and only if the SMA is accurate
and $\eta{=}0$ is used111111This last condition implies that the Voigt
distribution analytically reduces to the Lorentzian distribution (4). It is
worth mentioning that in the case of simple crystals, characterized by
$\kappa_{P}{\gg}\kappa_{C}$, the term $s^{\prime}{=}s$ in Eq. (8) yields the
well known expression obtained from Peierls’s theory, which interprets the
diffusivity (averaged over the three Cartesian directions) as
$D(\bm{q})_{s}=\frac{1}{3}|\\!|\bm{\mathsfit{v}}(\bm{q})_{s,s}|\\!|^{2}\tau(\bm{q})_{s}$,
where $\bm{\mathsfit{v}}(\bm{q})_{s,s}$ is the free propagation velocity of
the particle-like heat carrier with wavevector $\bm{q}$ and mode $s$, and
$\tau(\bm{q})_{s}$ is the inter-collision time. In finite-size models of
glasses it is most informative to represent the diffusivity as a function of
frequency, first because the specific heat of a vibration depends on its
frequency $\omega$ ($C(\bm{q})_{s}{=}C(\omega(\bm{q})_{s})$, see Eq. (2)), and
second because the vibrational frequencies determine measurable quantities
such as the vibrational density of states (in contrast, as mentioned before,
in finite-size models of glasses quantities such as the wavevectors $\bm{q}$
span a BZ that depends on the model and are used only as a mathematical tool
in the determination of the bulk limit). Thus, we represent the thermal
diffusivity as a function of frequency with
$D(\omega,T){=}[g(\omega){\mathcal{V}N_{\rm
c}}]^{-1}\sum_{\bm{q},s}D(\bm{q})_{s}\delta(\omega{-}\omega(\bm{q})_{s})$
(here $g(\omega){=}({\mathcal{V}N_{\rm
c}})^{-1}\sum_{\bm{q},s}\delta(\omega{-}\omega(\bm{q})_{s})$ is the
vibrational density of states (vDOS), which can be considered independent from
temperature, as also shown in Fig. 12 in the Methods; the Dirac $\delta$ is
broadened with a Gaussian distribution having a broadening determined from the
convergence test discussed in Sec. III.1). In the low-temperature and
infinite-reference-cell limit, Eq. (8) reduces to the temperature-independent
harmonic diffusivity introduced by Allen and Feldman 8 (this follows from the
properties of the Voigt distribution
$\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}},\eta]}$
discussed in Sec. III.2); in the following the dependence from temperature
will be shown explicitly for the sake of clarity. In this frequency-dependent
representation the conductivity reads

$\kappa(T)=\int_{0}^{\infty}g(\omega)C(\omega,T)D(\omega,T)d\omega\;,$ (9)

and intuitively allows to resolve the contribution of vibrations with
frequency $\omega$ to heat transport through their density of states
$g(\omega)$, the heat carried $C(\omega,T)$, and the diffusion rate
$D(\omega,T)$.

\begin{overpic}[width=216.81pt]{./8_diffusivity.pdf} \put(-2.0,96.0){{a)}}
\put(-2.0,61.0){{b)}} \put(-2.0,40.0){{c)}} \end{overpic}
Figure 8:  Vibrational DOS, specific heat, and diffusivity of vitreous silica.
a) Vibrational density of states of vitreous silica. Red and green are the
vDOS for the 192(D)54 model (192-atom model generated from first-principles
DFT simulations), obtained evaluating the frequencies of such model either
from first principles or from the GAP potential, and using a
$3{\times}3{\times}3$ $\bm{q}$ mesh. The dashed blue line and gray area are
the vDOS of the GAP-generated49 192- and 5184-atom models, computed using GAP
as in Fig. 4 (i.e. relying on a $3{\times}3{\times}3$ $\bm{q}$ mesh in the
former case, and evaluating the diffusivity at $\bm{q}{=}\bm{0}$ only in the
latter case). b) Quantum harmonic specific heat as a function of frequency and
temperature, obtained evaluating the expression (2) (solid lines). The dashed
line represents the classical specific heat (obtained from the quantum
specific heat in the infinite-temperature limit). c) Harmonic Allen-Feldman
thermal diffusivity computed from Eq. (8) using
$\Gamma(\bm{q})_{s}{=}0\;\forall\bm{q},s$ and with the models discusses in
panel a) (same parameters and color code); the inset shows that in $v$-SiO2
(5184(G) model) the AF diffusivity is very similar to the anharmonic
diffusivity (full Eq. (8)).

We report in Fig. 8 all these quantities. In panel a) we show the vDOS for the
192(D) $v$-SiO2 model computed using first principles calculations (red) or
the GAP potential (green). The main difference between these calculations is
that the GAP potential slightly stretches the high-frequency part of the
spectrum towards higher values; the effects on the thermal conductivity of
such a slight stretch of vibrational energies is barely appreciable, as shown
in Fig. 5. The dashed-blue line and the gray area are the vDOS for the 192(G)
and 5184(G) models, respectively; both these vDOS are obtained using GAP, the
former is computed using a $3{\times}3{\times}3$ $\bm{q}$ interpolation mesh,
while the latter is computed at $\bm{q}{=}\bm{0}$ only. The similarity between
these two curves further supports the usage of the $\bm{q}$ interpolation
technique to sample more accurately the vibrations in a glass. Fig. 8a) shows
that the vDOS is slightly affected by the method used to generated the
$v$-SiO2 model, with the 192(G) and 5184(G) models generated using the melt-
quench method 49 leading to vDOS slightly different from the 192(D) model,
which was generated using the bond-switching technique51, 54. In panel b) of
Fig. 8 we report the quantum harmonic specific heat as a function of
frequency, showing how increasing temperature populates vibrational modes of
increasingly larger frequency that consequently contribute to transport. We
highlight how, for all the temperatures considered, the quantum specific heat
differs from the constant classical (equipartition) limit (obtained letting
$T{\to}\infty$ in Eq. (2)); this reiterates the role of quantum statistics of
vibrations in thermal transport 35. Panel c) shows the harmonic AF diffusivity
(Eq. (8) using $\Gamma(\bm{q})_{s}{=}0\;\forall\bm{q},s$) for the models of
$v$-SiO2 discussed in panel a) (same parameters and color code are used).
Considerations analogous to those for the vDOS in panel a) hold: the method
used to generated the $v$-SiO2 model has a small but appreciable effect on the
diffusivity, with the 192(D) model generated using the bond-switching method
having a larger diffusivity at low frequency. Overall, the combined variations
of vDOS and diffusivity due to the technique used to generate the model, or
the approach (first-principles calculations or GAP) used to evaluate
frequencies and velocity operators, yield differences on the thermal
conductivity within 9$\%$, as shown before in Fig. 6. The inset of panel c)
shows that the anharmonic diffusivity (colored lines, computed from Eq. (8))
changes very little with temperature, and is practically very similar to the
temperature-independent AF diffusivity (gray area, obtained using Eq. (8) with
$\Gamma(\bm{q})_{s}{=}0\;\forall\bm{q},s$ and $\eta$ determined from the
convergence test detailed in Fig. 2); the inset shows calculations with GAP on
the 5184(G) model, analogous considerations hold for the other models studied
from first-principles or using GAP.

##  VII Conclusions

We have discussed a computational protocol that allows to determine from
finite-size models of glasses containing less than 200 atoms — thus within the
reach of standard first-principles approaches — the bulk limit of the harmonic
Allen-Feldman conductivity8, 9, as well as of the anharmonic Wigner
conductivity36, 37. To determine the bulk limit of the harmonic AF
conductivity the protocol employs Fourier interpolation to improve the
sampling of the vibrational spectrum of the glass model, and it numerically
represents the Dirac $\delta$ appearing in the AF conductivity with a light-
tailed Gaussian broadening larger than the average vibrational energy-level
spacing (in place of the originally proposed fat-tailed Lorentzian
broadening8, 9). To evaluate the bulk limit of the Wigner conductivity, the
protocol uses a Voigt profile — a two-parameter distribution obtained as a
convolution between the Gaussian used in the AF calculation, and the
Lorentzian with FWHM determined by the linewidths appearing in the Wigner
conductivity. The Voigt profile ensures that the linewidths affect the
conductivity only when they are larger than the broadening used in the AF
calculation; this allows to retain in finite-size models of glasses the
physical property that heat transfer via a wave-like tunneling between
neighboring (quasi-degenerate) vibrational eigenstates can occur even in the
limit of vanishing anharmonicity, provided these eigenstates are not Anderson-
localized, i.e. that they are coupled by non-zero velocity-operator elements.

The protocol has been validated on the paradigmatic glass $v$-SiO2, using a
state-of-the-art GAP potential and atomistic models containing 192 or 5184
atoms49 to compute vibrational states; we have shown that employing the
protocol on a 192-atom model allows to obtain harmonic (AF) and anharmonic
(rWTE) conductivities in perfect agreement with those of this much larger,
5184-atom model generated using the same technique.

After validation, we have used the protocol to predict the AF and rWTE
conductivities of $v$-SiO2 fully from first-principles. We have shown that
anharmonicity does not significantly affect the conductivity of $v$-SiO2, even
at high temperatures, since the AF conductivity is very similar to the rWTE
conductivity over the entire temperature range analyzed ($30{<}T{<}1300$ K).
We have supported this finding by showing that the rWTE conductivity, which
accounts for anharmonicity at the lowest (third) perturbative order 37, is
compatible at high temperatures with the conductivity obtained from first-
principles GKMD from Ercole et. al.27 (we recall that GKMD simulations are
accurate at high temperature, where the quantum specific heat approaches the
classical limit and anharmonic effects are maximized). Our calculations are in
agreement with experiments in the temperature range
$50{\lesssim}T{\lesssim}450$ K, but do not describe the surge of the
conductivity observed at higher temperatures. Future work will aim at
understanding if such a discrepancy can be related, as it seems likely, to
radiative effects98.

The results obtained for $v$-SiO2 have allowed us to gain general insights on
how anharmonicity affects the thermal conductivity of glasses at high
temperature. Specifically, we have shown that the high-temperature trend of
the conductivity is determined by how off-diagonal velocity-operator elements
— which couple pairs of vibrational eigenstates $(\bm{q})_{s}$ and
$(\bm{q})_{s^{\prime}}$, allowing tunnelling between them — vary as a function
of the energy difference
$\hbar\omega_{d}{=}\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}}$.
Velocity-operator elements increasing (decreasing) with $\hbar\omega_{d}$
drive a conductivity increase (decrease) with temperature. In the specific
case of $v$-SiO2, the saturating trend of the conductivity derives from
velocity-operator elements that are constant with respect to the energy
difference between the eigenstates coupled.

Finally, we have interpreted heat conduction in terms of frequency-resolved
thermal diffusivities for vibrations, showing that the harmonic Allen-Feldman
diffusivity characterizes accurately thermal transport in vitreous silica at
the microscopic level. This work paves the way to study the thermal
conductivity of glasses from first principles, and will be particularly
relevant to investigate the thermal properties of amorphous materials for
which developing quantum-accurate interatomic potentials is particularly
challenging or unpractical.

## Acknowledgements

We thank G. Csányi and V. L. Deringer for useful discussions. N. M.
acknowledges funding from the Swiss National Science Foundation under the
Sinergia grant no. 189924. M. S. acknowledges support from Gonville and Caius
College, and from the SNSF project P500PT_203178. Part of the calculations
presented in this work have been performed using computational resources
provided by the Cambridge Tier-2 system operated by the University of
Cambridge Research Computing Service (www.hpc.cam.ac.uk) funded by EPSRC
Tier-2 capital grant EP/T022159/1.

Figure 9: Thermal conductivity of crystalline $\alpha$-cristobalite and
$\alpha$-quartz. For these materials the total bulk conductivity resulting
from the full solution of the WTE (solid lines) is practically
indistinguishable from the bulk conductivity computed relying on the SMA
approximation (dashed lines). The $\alpha$-$\beta$ transition temperature for
cristobalite is around 540 K 102, below this temperature both
$\alpha$-cristobalite and $\alpha$-quartz are to a good approximation simple
crystals, since their total conductivity is predominantly determined by the
particle-like contribution ($\kappa_{{}_{\rm P}}$, dotted lines). Scatter
points are experiments by Eucken 103, Birch et al. 104, Knapp 105, Zeller et
al. 106 and Kanamori et al. 84.

##  Appendix A Accuracy of the SMA approximation in glasses

Glasses feature a low thermal conductivity, originating from strong scattering
of vibrations due to anharmonicity or disorder. Past works73, 107, 108, 37
have shown that as scattering becomes stronger, the relaxation to equilibrium
of the vibrational excitations becomes faster, and therefore the SMA — which
approximates the collision operator as in Eq. 4 (see Ref.37) — becomes more
accurate. These considerations already suggest that the SMA approximation is
accurate in vitreous silica. To further support these expectations, we
investigated numerically the accuracy of the SMA approximation in the
crystalline silica polymorphs $\alpha$-cristobalite and $\alpha$-quartz. These
crystals have a conductivity much larger than vitreous silica, implying a
weaker scattering from disorder or anharmonicity; consequently in
$\alpha$-cristobalite and $\alpha$-quartz the SMA is expected73, 107, 108 to
be less accurate than in vitreous silica. Fig. 9 shows that for
$\alpha$-cristobalite and $\alpha$-quartz solving the WTE in full or employing
the SMA approximation yields conductivities practically indistinguishable;
this demonstrates that the SMA is accurate in $\alpha$-cristobalite and
$\alpha$-quartz, and also suggests that the SMA is even more accurate in
vitreous silica.

##  Appendix B Effects of disorder and temperature on the linewidths

In Fig. 10a) we highlight how the frequency-linewidth distributions for models
of vitreous silica having different size (192 or 108 atoms) are overlapping,
suggesting that the linewidths of these models are not significantly affected
by finite-size effects. In panel b) and c) we show the linewidths of the
crystalline polymorphs $\alpha$-cristobalite (containing 12 atoms per
primitive cell) $\alpha$-quartz (9 atoms per primitive cell). Clearly, the
variation of the anharmonic linewidths is mainly due to temperature, since at
fixed temperature the frequency-linewidth distributions of amorphous and
crystalline silica polymorphs have a similar magnitude. We also note that the
anharmonic linewidths computed with GAP for the 192(G) model (Fig. 3) are very
similar to the anharmonic linewidths computed from first-principles for the
192(D) and 108(D) models (Fig. 10a)).

Figure 10: Effect of disorder and temperature on the linewidths of silica
polymorphs. Increasing temperature yields an increase of the linewidths that
account for third-order anharmonicity 72, 73 and isotopic scattering 76 in
both vitreous and crystalline materials: the left panel shows the 192(D) and
108(D) vitreous-silica model (both computed from first-principles and at
$\bm{q}{=}\bm{0}$ only), the central panel shows $\alpha$-cristobalite, and
the right panel shows $\alpha$-quartz. The temperatures of 50, 300, and 500 K
at which the linewidths distributions are reported are chosen to span the
temperature range over which all these materials are stable. The 1300 K
distribution shows the behavior at high temperature, and is reported only for
vitreous silica since $\alpha$-cristobalite and $\alpha$-quartz are not stable
at this temperature. The insets in panel b) and c) show the linewidths at low
vibrational energies for $\alpha$-cristobalite and $\alpha$-quartz,
respectively. The purple area shows the overdamped regime characterized by
$\hbar\Gamma>\hbar\omega$; the lack of linewidths in the purple region shows
that there are no overdamped vibrations and thus the Wigner formulation can be
employed 37, 42. The gray lines represent the average spacing between the
vibrational energy levels Eq. (6) (in panel a) only the average energy-level
spacing for the 192-atom model is reported).

##  Appendix C Accounting for anharmonicity at a reduced computational cost

In this section we discuss the details of the computation of the analytical
function $\Gamma_{a}[\omega]$, used to approximatively determine the
linewidths as a function of frequency discussed in Fig. 3. The analytical
function $\Gamma_{a}[\omega]$ is determined as

$\Gamma_{a}[\omega]=\frac{1}{\sqrt{\frac{1}{(\Gamma_{1}[\omega])^{2}}+\frac{1}{(\Gamma_{2}[\omega])^{2}}}},$
(10)

where $\Gamma_{1}[\omega]$ and $\Gamma_{2}[\omega]$ are defined as

$\begin{split}&\Gamma_{1}[\omega]{=}\frac{\sum\limits_{\bm{q}=\bm{0},s}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\Big{[}-\frac{\hbar^{2}(\omega(\bm{q})_{s}-\omega)^{2}}{2\sigma^{2}}\Big{]}}{\sum\limits_{\bm{q}=\bm{0},s}\tau(\bm{q})_{s}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\Big{[}-\frac{\hbar^{2}(\omega(\bm{q})_{s}-\omega)^{2}}{2\sigma^{2}}\Big{]}},\\\
&\Gamma_{2}[\omega]{=}p\cdot\omega^{2},\\\
&p=\frac{\sum\limits_{\bm{q}=\bm{0},s}\int_{\omega_{\rm o}}^{2\omega_{\rm
o}}d\omega_{c}\frac{\Gamma(\bm{q})_{s}}{\omega^{2}(\bm{q})_{s}}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\\!\Big{[}{-}\tfrac{\hbar^{2}(\omega(\bm{q})_{s}-\omega_{c})^{2}}{2\sigma^{2}}\Big{]}}{\sum\limits_{\bm{q}=\bm{0},s}\int_{\omega_{\rm
o}}^{2\omega_{\rm
o}}d\omega_{c}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\\!\Big{[}{-}\tfrac{\hbar^{2}(\omega(\bm{q})_{s}-\omega_{c})^{2}}{2\sigma^{2}}\Big{]}}.\end{split}$
(11)

$\omega_{\rm o}$ is the smallest non-zero frequency at $\bm{q}=\bm{0}$ and
$\sigma{=}15$ cm-1 is a broadening chosen sufficiently large to ensure that
the linewidths are averaged in a smooth way. The functional form of the
approximated function $\Gamma_{a}[\omega]$ is inspired by past work 77, 78,
and the specific expressions (10,11) to determine it have been devised and
validated relying on exact calculations performed on the 108(D) $v$-SiO2
model.

Figure 11: Effect of $\bm{q}$-interpolation on the linewidths and
approximation to reduce the computational cost. Panel a), linewidths
distribution of the 108(D) model of $v$-SiO2 computed explicitly on a
$5{\times}5{\times}5$ $\bm{q}$-mesh (light blue is at 30 K, light green is at
200 K, and orange is at 1300 K) or computed at the point $\bm{q}{=}\bm{0}$
only (dark blue is at 30 K, dark green is at 200 K, and red is at 1300 K). The
purple region represents the overdamped regime, where vibrations cannot be
accurately described using the Wigner formulation and spectral-function
approaches have to be employed 37, 42. The horizontal dotted line is the
broadening $\hbar\eta$ used in the Voigt renormalization for the 108-atom
model, all the linewidths below this line are regularized and thus have
negligible effect on the rWTE conductivity. The dashed lines are the
analytical functions $\Gamma_{a}[\omega]$, determined from the distributions
at $\bm{q}{=}\bm{0}$ only as detailed in Eqs. (10,11). Panel b) shows the bare
WTE conductivity of the 108(D) model computed exactly (i.e. using the
linewidths explicitly computed on the $5{\times}5{\times}5$ $\bm{q}$-mesh,
solid red), or using the linewidths approximatively determined using the
function $\Gamma_{a}[\omega]$ (dashed blue); the solid green and dotted black
lines show the rWTE conductivities computed using the exact or approximated
linewidths, respectively. The good agreement between the exact and the
approximated calculations shows that the approximation developed allows to
account for anharmonicity at a reduced computational cost and without
appreciably compromising accuracy, both in the WTE and rWTE calculations.

Specifically, we show in Fig. 11a) that the approximated functions
$\Gamma_{a}[\omega]$ (dashed lines, different colors show different
temperatures) captures the trend of the linewidths explicitly computed over a
dense $5{\times}5{\times}5$ $\bm{q}$ mesh (dense distributions). We recall
that $\Gamma_{a}[\omega]$ is determined from the linewidth distributions at
$\bm{q}{=}\bm{0}$ only (coarse clouds of scatter points). In Fig. 11b) we show
that computing the bare WTE conductivity over a dense $5{\times}5{\times}5$
$\bm{q}$ mesh using the anharmonic linewidths computed exactly (solid red) or
determined approximatively using the function $\Gamma_{a}[\omega]$ (dashed
blue) are practically indistinguishable; consequently, also the rWTE
conductivity computed over the same $5{\times}5{\times}5$ $\bm{q}$ mesh is
practically unchanged when the exact (solid green) or approximated (dotted
black) linewidths are used.

##  Appendix D Quantum harmonic specific heat

We show in Fig. 12 that the theoretical quantum harmonic specific heat at
constant volume ($C_{\rm V}^{\rm Th}(T)=\frac{1}{\rho\mathcal{V}N_{\rm
c}}\sum_{\bm{q},s}C(\bm{q})_{s}$, where $\rho$ is the density) is in close
agreement with the experimental specific heat at constant pressure 109. This
suggests that the renormalization of vibrational energies due to anharmonicity
and temperature 110, 111, 41, 47 are negligibly small for $v$-SiO2 in the
temperature range considered 112, and as such these effects are not considered
in this work.

Figure 12: Specific heat of vitreous silica. Lines are the quantum harmonic
specific heat at constant volume and as a function of temperature computed
from first principles for the 108(D) (green), 144(D) (blue) and 192(D) (red)
vitreous-silica models. Circles are measurements 109 of the specific heat at
constant pressure.

##  Appendix E Velocity operator as a function of frequency

In this section we provide all the details on the plot reported in Fig. 7. In
order to gain insights into the increasing, non-saturating trend of the
conductivity with temperature (Fig. 5), we recast the conductivity formula (1)
as a function of $\omega_{\rm
a}{=}(\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}})/2$ and $\omega_{\rm
d}{=}\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}}$:

$\begin{split}\kappa=&\int\limits_{0}^{\omega_{\rm max}}d\omega_{\rm
a}\int\limits_{-\omega_{\rm max}}^{\omega_{\rm max}}d\omega_{\rm
d}\Bigg{[}\frac{1}{\mathcal{V}N_{\rm
c}}{\sum_{\bm{q},s,s^{\prime}}}\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}}{4}\frac{\rVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\lVert^{2}}{3}\\\
\times&\left(\frac{C(\bm{q})_{s}}{\omega(\bm{q})_{s}}{+}\frac{C(\bm{q})_{s^{\prime}\\!}}{\omega(\bm{q})_{s^{\prime}\\!}}\right)\\!\pi\mathcal{F}_{[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}},\eta]}(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})\\\
\times&{\delta}_{\sigma_{\rm
a}}\left(\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}}{2}{-}\omega_{\rm
a}\right){\delta}_{\sigma_{\rm
d}}\big{(}(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}}){-}\omega_{\rm
d}\big{)}\Bigg{]}\end{split}$ (12)

where the distributions $\delta_{\sigma}$ are Gaussian broadenings of the
Dirac delta:

$\delta_{\sigma}\left(\Omega-\omega\right)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[{-\frac{1}{2\sigma^{2}}\left(\Omega-\omega\right)^{2}}\right].$
(13)

In order to achieve our goal to recast Eq. (12) in terms of physically-
insightful frequency-dependent functions, we approximate the linewidths as
single-valued function of frequency, i.e.
$\Gamma(\bm{q})_{s}=\mathcal{L}(\omega(\bm{q})_{s})$; using this approximation
allows to recast Eq. (12) as follows:

$\begin{split}\kappa=\int\limits_{0}^{\omega_{\rm max}}d\omega_{\rm
a}\int\limits_{-\omega_{\rm max}}^{\omega_{\rm max}}&d\omega_{\rm
d}\mathcal{G}(\omega_{\rm a},\omega_{\rm d})\mathcal{C}(\omega_{\rm
a},\omega_{\rm d})\left<|V^{\rm avg}_{\omega_{\rm a},\omega_{\rm
d}}|^{2}\right>\\\ \times\pi&\mathcal{F}_{[\mathcal{L}(\omega_{\rm
a}+\frac{\omega_{\rm d}}{2}){+}\mathcal{L}(\omega_{\rm a}-\frac{\omega_{\rm
d}}{2}),\eta]}(\omega_{\rm d}),\end{split}$ (14)

where $\mathcal{G}(\omega_{\rm a},\omega_{\rm d})$ is a density of states

$\begin{split}\mathcal{G}(\omega_{\rm a},\omega_{\rm d})=\frac{1}{N_{\rm
at}}\frac{1}{\mathcal{V}N_{\rm
c}}&{\sum_{\bm{q},s,s^{\prime}}}{\delta}_{\sigma_{\rm
a}}\left(\frac{\omega(\bm{q})_{s}+\omega(\bm{q})_{s^{\prime}}}{2}-\omega_{\rm
a}\right)\\\ \times&{\delta}_{\sigma_{\rm
d}}\big{(}(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})-\omega_{\rm
d}\big{)}\;,\end{split}$ (15)

$\mathcal{C}(\omega_{\rm a},\omega_{\rm d})$ is a specific heat

$\mathcal{C}(\omega_{\rm a},\omega_{\rm d})=\frac{\omega_{\rm
a}}{2}\left(\frac{C(\omega_{\rm a}+\frac{\omega_{\rm d}}{2})}{\omega_{\rm
a}+\frac{\omega_{\rm d}}{2}}+\frac{C(\omega_{\rm a}-\frac{\omega_{\rm
d}}{2})}{\omega_{\rm a}-\frac{\omega_{\rm d}}{2}}\right),\;$ (16)

and $\left<|V^{\rm avg}_{\omega_{\rm a},\omega_{\rm d}}|^{2}\right>$ is the
average square modulus of the velocity operator defined in Eq. (7) (whose
Dirac delta must be broadened with the Gaussian $\delta_{\sigma}$ discussed
here) and plotted in Fig. 7. Eq. (14), together with Fig. 7 and Fig. 10, sheds
light on the saturating trend of the conductivity with temperature (Fig. 5a)).
In fact, among the various quantities entering in Eq. (14), the density of
states (15) and the specific heat (16) have a trivial temperature dependence
(the former is temperature-independent, the latter saturates with increasing
temperature). The change of variable performed in Eq. (12) shows that the
temperature-conductivity trend is determined by how the average square
velocity-operator elements vary with the vibrational frequency difference
$\omega_{\rm d}$, because the increase of the linewidths with temperature
(Fig. 10a)) results in a broader Lorentzian distribution (4) that encloses
velocity-operator elements corresponding to increasingly larger frequency
differences $\omega_{\rm d}$. In particular, for vitreous silica, the average
square velocity-operator elements are almost independent from $\omega_{\rm d}$
for all values of $\omega_{\rm a}$ (Fig. 7). It follows that the saturating
trend of the conductivity at high temperature reported in Fig. 5a) is
inherited from the saturating trend of the specific heat.

Fig. 7 is computed relying on the 192(D) model and on first-principles
calculations, using a $3{\times}3{\times}3$ $\bm{q}$-point mesh,
$\hbar\sigma_{\rm a}=15$ cm-1 and $\hbar\sigma_{\rm
d}=\sqrt{\frac{\pi}{2}}\hbar\eta=5$ cm-1 (this latter corresponds to a
Gaussian having height $\frac{1}{\pi\eta}$, with $\hbar\eta=4$ cm-1 equal to
the one used in the computation of the AF conductivity for the 192(D) model).
Increasing the $\bm{q}$-point mesh to $5{\times}5{\times}5$ or multiplying
$\sigma_{\rm a}$ and $\sigma_{\rm d}$ by a factor of 2 does not produce
significant changes. We have verified that plotting the velocity operator for
the other models studied from first principles or using GAP yields results
that are practically indistinguishable from those reported in Fig. 7.

##  Appendix F Effects of anharmonicity using the BKS potential

While our BKS-based rWTE predictions are in agreement with Ref.68, our BKS-
based AF conductivity is significantly different from that reported in Ref.
68, which claimed that anharmonicity enhances appreciably the conductivity in
$v$-SiO2 (specifically, Ref.68 concluded that in $v$-SiO2 in the high-
temperature limit, the anharmonic WTE conductivity is about $30\%$ higher than
the AF conductivity). In Fig. 13 we show that evaluating the AF conductivity
as done in Ref.68, i.e. using the BKS potential, at $\bm{q}=\bm{0}$ only and
using a Lorentzian broadening $\hbar\eta=1.6\;{\rm cm}^{-1}$, yields a value
that significantly underestimates the bulk value for the AF conductivity.
Using the Fourier interpolation with an aptly chosen broadening to determine
the bulk limit ($\hbar\eta=4\;{\rm cm}^{-1}$) yields a larger AF conductivity,
which differs only about 4$\%$ from the bulk limit of the rWTE conductivity.

Figure 13:  Negligible effects of anharmonicity in $v$-SiO2 described with the
BKS potential. Evaluating the conductivity using the protocol to determine the
bulk limit yields an anharmonicity rWTE conductivity (solid purple) that
differs only about 4$\%$ from the bulk limit of the AF conductivity (dotted
purple, evaluated using a Gaussian broadening $\hbar\eta=4\;{\rm cm}^{-1}$
determined from a convergence test analogous to that reported in Fig. 2). In
contrast, evaluating the AF and WTE conductivities at $\bm{q}{=}\bm{0}$ only
(and using a Lorentzian broadening $\hbar\eta=1.6\;{\rm cm}^{-1}$ in the
former case) yield values that strongly underestimates the bulk value for the
AF conductivity (dotted orange) and weakly underestimates the Wigner
conductivity (solid orange), respectively.

##  Appendix G Calculations performed using the GAP potential

The GAP potential for vitreous silica was taken from Ref. 49. The 5184(G) and
192(G) models discussed in Figs. 1,2,3 are taken from Ref. 49, we chose the
192-atom structure number 4 in that reference as the one having the closest
density to the large 5184-atom model. More details on these structures are
reported in Ref. 49. The densities of the models analyzed are reported in Tab.
1. We have checked that the AF conductivity of the other four 192-atom models
of $v$-SiO2 discussed in Ref. 49 are very similar to the bulk AF
conductivities shown in Fig. 1 (differences are smaller than the $9\%$
differences discusses in Fig. 6).

Interatomic forces and stress tensor are computed using lammps113, using the
quip114, 115, 116 interface to call the GAP potential routines. Cell
parameters and atomic positions are relaxed using a threshold of $25$
eV/Angstrom for the atomic forces (i.e. a structure is considered relaxed if
all the Cartesian components of the forces acting on atoms are smaller than
this threshold), and a threshold of 0.1 kBar for pressure. The harmonic
dynamical matrices (which yield the vibrational frequencies and velocity
operators) are computed using phonopy, on a $2{\times}2{\times}2$ supercell
for the 192(G) structure, and on the reference cell ($1{\times}1{\times}1$)
for the 5184(G) model. For the 192(G) model, third-order force constants are
computed in the reference cell using ShengBTE 73, using a cutoff equal to the
$6^{\rm th}$ nearest-neighbor; the resulting force constants are then
converted in phono3py75 format using hiphive117. The anharmonic linewidths for
the 192(G) model (shown in Fig. 3) are computed using phono3py at
$\bm{q}=\bm{0}$ and using Gaussian smearing of $\hbar\sigma=2\;{\rm cm}^{-1}$
for the Dirac delta appearing in the collision operator. The anharmonic
linewidths for the 5184(G) model are determined using the approximation
detailed in Fig. 3 and in Fig. 11 in the Methods.

##  Appendix H First-principles calculations

All the density-functional theory (DFT) calculations have been performed with
the Quantum espresso distribution 118 using the PBE functional with Grimme-D2
corrections (PBE+D2) 119. This choice is motivated by the benchmarks given in
Ref. 120 and accounting for the capability of Quantum espresso to compute
phonons using density-functional perturbation theory (DFPT) 121 with the
PBE+D2 exchange-correlation (XC) functional. This choice is validated by the
agreement between theoretical and experimental densities reported in table 1
and also by the capability of this functional to accurately describe the
thermal properties of $\alpha$-quartz (Fig. 9). We have used pseudopotentials
from the SSSP efficiency library 122, 123 with a cutoff of $50$ Ry and a dual
of 8. In the following we report the details for all the three different
systems studied: vitreous silica, $\alpha$-cristobalite, and $\alpha$-quartz.

Table 1: Densities of silica polymorphs. Comparison between theoretical and experimental densities $\rho$ for the various silica polymorphs analyzed.  Structure |  $\rho$ ($\rm kg/m^{3}$)
---|---
$v$-SiO2, 192(D) PBE+D2 (SiO2 2818 54) | 2241.2
$v$-SiO2, 192(D) GAP (SiO2 2818 54) | 2288.6
$v$-SiO2, 192(D) BKS (SiO2 2818 54) | 2243.8
$v$-SiO2, 192(G) GAP ($\\#4$ 49) | 2188.7
$v$-SiO2, 5184(G) GAP 49 | 2203.6
$v$-SiO2, 144(D) PBE+D2 52 | 2220.6
$v$-SiO2, 108(D) PBE+D2 (SiO2.1586 51) | 2243.9
$v$-SiO2, Experiment 124 | 2203 $\pm$ 3
$v$-SiO2, Experiment 91 | 2200 $\pm$ 10
$\alpha$-cristobalite 125, PBE+D2 | 2383.5
$\alpha$-cristobalite, Experiment 126 | 2326$\pm$12
$\alpha$-quartz 127, PBE+D2 | 2641.9
$\alpha$-quartz, Experiment 128 | 2650

Vitreous silica. The 108(D) structure of vitreous silica is taken from Ref. 51
(“SiO2.1586” model). For this structure, cell parameters and atomic positions
are relaxed with DFT using $\Gamma$-point ($\bm{q}{=}\bm{0}$) sampling, a
threshold of $10^{-4}$ Ry/Bohr for the atomic forces (i.e. a structure is
considered relaxed if all the Cartesian components of the forces acting on
atoms are smaller than this threshold), and a threshold of 0.1 kBar for
pressure. The harmonic dynamical matrices (which yield the vibrational
frequencies and velocity operators) are computed using DFPT on a
$2{\times}2{\times}2$ $\bm{q}$-point mesh and accounting for the non-analytic
term correction due to the dielectric tensor and Born effective charges.
Third-order force constants are computed in the reference cell using ShengBTE
73, and with a cutoff of $0.32$ nm (corresponding to the $6^{\rm th}$ nearest-
neighbor). The linewidths are computed using phono3py 75 on a
$5{\times}5{\times}5$ $\bm{q}$-point mesh using the tetrahedron method, we
checked that computing them using Gaussian smearing of $\hbar\sigma=2\;{\rm
cm}^{-1}$ for the Dirac delta appearing in the collision operator does not
produce appreciable changes. Thermal conductivity calculations are performed
using a $\bm{q}$-interpolation mesh equal to $5{\times}5{\times}5$. The Voigt
profile, used to combine the AF and the WTE conductivities as detailed before,
has been numerically implemented following the prescriptions reported in Ref.
82. We have checked that reducing the $\bm{q}$-interpolation mesh to
$3{\times}3{\times}3$ does not produce appreciable changes.

The 192(D) model is generated using the same techniques discussed in Ref. 51
and is discussed in Ref. 54, while the 144(D) model is taken from Ref. 52 and
is available in the Materials Cloud Archive 53. For all these structures, the
relaxation of the reference cell with DFT, the calculation of the harmonic
dynamical matrices, and the the calculation of third-order anharmonic force
constants are performed using the same parameters used for the 108-atom
structure. These harmonic and anharmonic force constants are then used to
compute the linewidths at the point $\bm{q}=\bm{0}$ only using phono3py 75 and
a Gaussian smearing of $2\;{\rm cm^{-1}}$, and then employed within the
approximated procedure discussed and validated in Fig. 11. We also verified
that the 192(D) model accurately reproduces the experimental bulk modulus of
$v$-SiO2 (the theoretical value we computed is 36.9 GPa, while the
experimental value of Ref. 129 is 36.8 GPa and the experimental value of Ref.
130 is 36.9 GPa).

An analysis of the coordination numbers using the procedure based on the
minimum of the radial distribution function as implemented in the r.i.n.g.s.
software 131 or that based on the position of the Wannier centers discussed in
Ref. 132 have revealed in both cases that the 108(D), 144(D), 192(D) vitreous
structures considered in this work do not have coordination defects or lone
pairs (both before and after the DFT relaxation).

$\alpha$-cristobalite. The crystal structure of $\alpha$-cristobalite is taken
from Ref. 125 (ICSD collection code 47219). In first-principles calculations,
the Brillouin zone is integrated with a Monkhorst–Pack mesh of
$5{\times}5{\times}4$ points, with a (1, 1, 1) shift. Second-order force
constants are computed using DFPT on a $5{\times}5{\times}4$ $\bm{q}$-point
mesh, accounting also for the non-analytic term correction due to the
dielectric tensor and Born effective charges. To obtain the third-order force
constants the finite-difference method implemented in ShengBTE 73 is used,
together with the interconversion software from ShengBTE to Quantum espresso,
available in the d3q package 108, 71. In these third-order force constants
calculations, a $2{\times}2{\times}2$ supercell with a $2{\times}2{\times}2$
k-point sampling is used, and interactions up to the $6^{\rm th}$ nearest
neighbor (corresponding to $\sim 0.39$ nm) are considered. The linewidths
(Fig. 10b)) and thermal conductivity (Fig. 9a)) are computed with the d3q
package 108, 71 using a $17{\times}17{\times}13$ $\bm{q}$-point mesh and a
Gaussian smearing $\hbar\sigma=4\;{\rm cm}^{-1}$.

$\alpha$-quartz. The crystal structure of $\alpha$-quartz is taken from Ref.
127 (Crystallographic Open Database id 1526860). In first-principles
calculations, the Brillouin zone is integrated with a Monkhorst–Pack mesh of
$6{\times}6{\times}5$ points, with a (1, 1, 1) shift. Second-order force
constants are computed using DFPT on a $4{\times}4{\times}4$ $\bm{q}$-point
mesh, accounting also for the non-analytic term correction due to the
dielectric tensor and Born effective charges. Third-order force constants are
computed using the finite difference methods as implemented in ShengBTE 73,
using a $3{\times}3{\times}2$ supercell with a $\Gamma$-only k-point sampling,
and a cutoff for atomic interactions of $0.9$ nm. The linewidths (Fig. 10c))
and thermal conductivity (Fig. 9b)) are computed with the d3q package 108, 71
using a $19{\times}19{\times}15$ $\bm{q}$-point mesh and a Gaussian smearing
$\hbar\sigma=4\;{\rm cm}^{-1}$. Results are compatible to those reported in
Ref. 133.

We conclude by noting that in Fig. 8 the vDOS for all materials are computed
using the same meshes used for the thermal conductivity calculations detailed
above, and a Gaussian smearing of 15 cm-1 for the vibrational energies.

Data availability. Raw data were generated using the SCITAS High Performance
Computing facility at the École Polytechnique Fédérale de Lausanne and the
Cambridge Service for Data-Driven Discovery (CSD3). The atomistic models of
vitreous silica, $\alpha$-quartz, and $\alpha$-cristobalite studied in this
work are available on the Materials Cloud Archive 97, 134.

Code availability. Quantum espresso118 is available at www.quantum-
espresso.org; the scripts related to the computation of the third order force
constants using the finite-difference method are available at
bitbucket.org/sousaw/thirdorder; phonopy and phono3py are available at
github.com/phonopy. The GAP potential for silica polymorphs is available at
49. lammps113 is available at www.lammps.org and the interface for lammps with
the GAP potential114, 115, 116 is available at github.com/libAtoms/QUIP.

## References

  * [1] Pasquarello, A., Hybertsen, M. S. & Car, R.  Interface structure between silicon and its oxide by first-principles molecular dynamics.  _Nature_ 396, 58 (1998).
  * [2] Niu, H., Piaggi, P. M., Invernizzi, M. & Parrinello, M.  Molecular dynamics simulations of liquid silica crystallization.  _Proc. Natl. Acad. Sci. USA_ 115, 5348–5352 (2018).
  * [3] Uyanna, O. & Najafi, H.  Thermal protection systems for space vehicles: A review on technology development, current challenges and future prospects.  _Acta Astronautica_ 176, 341–356 (2020).
  * [4] Hu, Y.-J. _et al._ Predicting densities and elastic moduli of SiO2-based glasses by machine learning.  _Npj Comput. Mater._ 6, 1–13 (2020).
  * [5] Kotz, F. _et al._ Three-dimensional printing of transparent fused silica glass.  _Nature_ 544, 337–339 (2017).
  * [6] Arbab, M. & Finley, J. J.  Glass in architecture.  _Int. J. Appl. Glass Sci._ 1, 118–129 (2010).
  * [7] Freeman, J. J. & Anderson, A. C.  Thermal conductivity of amorphous solids.  _Phys. Rev. B_ 34, 5684–5690 (1986).
  * [8] Allen, P. B. & Feldman, J. L.  Thermal Conductivity of Glasses: Theory and Application to Amorphous Si.  _Phys. Rev. Lett._ 62, 645–648 (1989).
  * [9] Allen, P. B. & Feldman, J. L.  Thermal conductivity of disordered harmonic solids.  _Phys. Rev. B_ 48, 12581–12588 (1993).
  * [10] Allen, P. B., Feldman, J. L., Fabian, J. & Wooten, F.  Diffusons, locons and propagons: Character of atomic vibrations in amorphous Si.  _Philos. Mag. B_ 79, 1715–1731 (1999).
  * [11] Leggett, A. J. & Vural, D. C.  “tunneling two-level systems” model of the low-temperature properties of glasses: Are “smoking-gun” tests possible?  _J. Phys. Chem. B_ 117, 12966–12971 (2013).
  * [12] Paz, A. P., Lebedeva, I. V., Tokatly, I. V. & Rubio, A.  Identification of structural motifs as tunneling two-level systems in amorphous alumina at low temperatures.  _Phys. Rev. B_ 90, 224202 (2014).
  * [13] Müller, C., Cole, J. H. & Lisenfeld, J.  Towards understanding two-level-systems in amorphous solids: insights from quantum circuits.  _Rep. Prog. Phys._ 82, 124501 (2019).
  * [14] Phillips, W. A.  Tunneling states in amorphous solids.  _J. Low Temp. Phys._ 7, 351–360 (1972).
  * [15] Anderson, P. W., Halperin, B. & Varma, C. M.  Anomalous low-temperature thermal properties of glasses and spin glasses.  _Philos. Mag._ 25, 1–9 (1972).
  * [16] Gurevich, V. L., Parshin, D. A. & Schober, H. R.  Anharmonicity, vibrational instability, and the boson peak in glasses.  _Phys. Rev. B_ 67, 094203 (2003).
  * [17] Schirmacher, W.  Thermal conductivity of glassy materials and the boson peak.  _EPL (Europhysics Letters)_ 73, 892 (2006).
  * [18] Lubchenko, V. & Wolynes, P. G.  The origin of the boson peak and thermal conductivity plateau in low-temperature glasses.  _Proc. Natl. Acad. Sci. U.S.A._ 100, 1515–1518 (2003).
  * [19] Shcheblanov, N. S., Povarnitsyn, M. E., Wiles, J. D., Elliott, S. R. & Taraskin, S. N.  Phonon traces in glassy vibrations.  _Phys. Rev. B_ 102, 024202 (2020).
  * [20] Shenogin, S., Bodapati, A., Keblinski, P. & McGaughey, A. J.  Predicting the thermal conductivity of inorganic and polymeric glasses: The role of anharmonicity.  _J. Appl. Phys._ 105, 034906 (2009).
  * [21] Donadio, D. & Galli, G.  Atomistic Simulations of Heat Transport in Silicon Nanowires.  _Phys. Rev. Lett._ 102, 195901 (2009).
  * [22] Larkin, J. M. & McGaughey, A. J. H.  Thermal conductivity accumulation in amorphous silica and amorphous silicon.  _Phys. Rev. B_ 89, 144303 (2014).
  * [23] Lv, W. & Henry, A.  Non-negligible contributions to thermal conductivity from localized modes in amorphous silicon dioxide.  _Sci. Rep._ 6, 35720 (2016).
  * [24] Marcolongo, A., Umari, P. & Baroni, S.  Microscopic theory and quantum simulation of atomic heat transport.  _Nat. Phys._ 12, 80 (2016).
  * [25] Carbogno, C., Ramprasad, R. & Scheffler, M.  Ab Initio Green-Kubo Approach for the Thermal Conductivity of Solids.  _Phys. Rev. Lett._ 118, 175901 (2017).
  * [26] Puligheddu, M., Gygi, F. & Galli, G.  First-principles simulations of heat transport.  _Phys. Rev. Mater_ 1, 060802 (2017).
  * [27] Ercole, L. _et al._ ab initio Simulation of Heat Transport in Silica Glass.  _PhD thesis, SISSA_ (2018).
  * [28] Verdi, C., Karsai, F., Liu, P., Jinnouchi, R. & Kresse, G.  Thermal transport and phase transitions of zirconia by on-the-fly machine-learned interatomic potentials.  _Npj Comput. Mater._ 7, 1–9 (2021).
  * [29] Jund, P. & Jullien, R.  Molecular-dynamics calculation of the thermal conductivity of vitreous silica.  _Phys. Rev. B_ 59, 13707–13711 (1999).
  * [30] Tian, Y. _et al._ Thermal conductivity of vitreous silica from molecular dynamics simulations: The effects of force field, heat flux and system size.  _J. Chem. Phys._ 146, 054504 (2017).
  * [31] Lampin, E., Palla, P. L., Francioso, P. A. & Cleri, F.  Thermal conductivity from approach-to-equilibrium molecular dynamics.  _J. Appl. Phys._ 114, 1–6 (2013).
  * [32] Martin, E., Ori, G., Duong, T.-Q., Boero, M. & Massobrio, C.  Thermal conductivity of amorphous SiO2 by first-principles molecular dynamics.  _J. Non-Cryst. Solids_ 581, 121434 (2022).
  * [33] Isaeva, L., Barbalinardo, G., Donadio, D. & Baroni, S.  Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach.  _Nat. Commun._ 10, 3853 (2019).
  * [34] Feldman, J. L., Kluge, M. D., Allen, P. B. & Wooten, F.  Thermal conductivity and localization in glasses: Numerical study of a model of amorphous silicon.  _Phys. Rev. B_ 48, 12589–12602 (1993).
  * [35] Puligheddu, M., Xia, Y., Chan, M. & Galli, G.  Computational prediction of lattice thermal conductivity: A comparison of molecular dynamics and Boltzmann transport approaches.  _Phys. Rev. Materials_ 3, 085401 (2019).
  * [36] Simoncelli, M. and Marzari, N. and Mauri, F.  Unified theory of thermal transport in crystals and glasses.  _Nat. Phys._ 15, 809–813 (2019).
  * [37] Simoncelli, M., Marzari, N. & Mauri, F.  Wigner formulation of thermal transport in solids.  _Phys. Rev. X, in press; arXiv preprint arXiv:2112.06897_ .
  * [38] Peierls, R. E.  Zur kinetischen theorie der wärmeleitung in kristallen.  _Annalen der Physik_ 395, 1055–1101 (1929).
  * [39] Peierls, R. E.  _Quantum theory of solids_ (Oxford Classics Series, 2001).
  * [40] Here the “diffuson” mechanisms described by Allen-Feldman theory is denoted with “wave-like tunneling” because it bears some analogies to the electronic Zener interband tunneling 59, more details about the terminology can be found in Refs. 36, 37.
  * [41] Xia, Y., Ozoliņš, V. & Wolverton, C.  Microscopic mechanisms of glasslike lattice thermal transport in cubic ${\mathrm{cu}}_{12}{\mathrm{sb}}_{4}{\mathrm{s}}_{13}$ tetrahedrites.  _Phys. Rev. Lett._ 125, 085901 (2020).
  * [42] Caldarelli, G., Simoncelli, M., Marzari, N., Mauri, F. & Benfatto, L.  Many-body green’s function approach to lattice thermal transport.  _Phys. Rev. B_ 106, 024312 (2022).
  * [43] Fiorentino, A. & Baroni, S.  From green-kubo to the full boltzmann kinetic approach to heat transport in crystals and glasses.  _arXiv preprint arXiv:2206.01279_ (2022).
  * [44] Broido, D., Malorny, M., Birner, G., Mingo, N. & Stewart, D.  Intrinsic lattice thermal conductivity of semiconductors from first principles.  _Appl. Phys. Lett._ 91, 231922 (2007).
  * [45] Garg, J., Bonini, N., Kozinsky, B. & Marzari, N.  Role of Disorder and Anharmonicity in the Thermal Conductivity of Silicon-Germanium Alloys: A First-Principles Study.  _Phys. Rev. Lett._ 106, 045901 (2011).
  * [46] Xia, Y., Pal, K., He, J., Ozoliņš, V. & Wolverton, C.  Particlelike Phonon Propagation Dominates Ultralow Lattice Thermal Conductivity in Crystalline ${\mathrm{Tl}}_{3}{\mathrm{VSe}}_{4}$.  _Phys. Rev. Lett._ 124, 065901 (2020).
  * [47] Jain, A.  Multichannel thermal transport in crystalline ${\mathrm{Tl}}_{3}{\mathrm{VSe}}_{4}$.  _Phys. Rev. B_ 102, 201201 (2020).
  * [48] Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ 104, 136403 (2010).
  * [49] Erhard, L. C., Rohrer, J., Albe, K. & Deringer, V. L.  A machine-learned interatomic potential for silica and its relation to empirical models.  _Npj Comput. Mater._ 8, 1–12 (2022).
  * [50] Sun, J., Ruzsinszky, A. & Perdew, J. P.  Strongly Constrained and Appropriately Normed Semilocal Density Functional.  _Phys. Rev. Lett._ 115, 036402 (2015).
  * [51] Charpentier, T., Kroll, P. & Mauri, F.  First-principles nuclear magnetic resonance structural analysis of vitreous silica.  _J. Phys. Chem. C_ 113, 7917–7929 (2009).
  * [52] Giacomazzi, L., Umari, P. & Pasquarello, A.  Medium-range structure of vitreous ${\text{SiO}}_{2}$ obtained through first-principles investigation of vibrational spectra.  _Phys. Rev. B_ 79, 064202 (2009).
  * [53] Giacomazzi, L., Umari, P. & Pasquarello, A.  Medium-range structure of vitreous ${\text{SiO}}_{2}$.  _Materials Cloud Archive_ (2019).  URL https://doi.org/10.24435/materialscloud:2019.0086/v1.
  * [54] Kroll, P.  _Modeling Amorphous Ceramic Structures_ , chap. 2, 39–69 (John Wiley & Sons, Ltd, 2013).
  * [55] van Beest, B. W. H., Kramer, G. J. & van Santen, R. A.  Force fields for silicas and aluminophosphates based on ab initio calculations.  _Phys. Rev. Lett._ 64, 1955–1958 (1990).
  * [56] Carré, A., Berthier, L., Horbach, J., Ispas, S. & Kob, W.  Amorphous silica modeled with truncated and screened Coulomb interactions: A molecular dynamics simulation study.  _J. Chem. Phys_ 127, 114512 (2007).
  * [57] Hinuma, Y., Pizzi, G., Kumagai, Y., Oba, F. & Tanaka, I.  Band structure diagram paths based on crystallography.  _Comput. Mater. Sci._ 128, 140–184 (2017).
  * [58] Here we have exploited the possibility to diagonalize at least one Cartesian components ${\alpha}$ of the velocity operator in the degenerate subspace 108, 37, thus ${\mathsfit{v}^{\alpha}(\bm{q})_{s,s}}$ is the propagation velocity in direction $\alpha$ and ${\mathsfit{v}^{\alpha}(\bm{q})_{s,s^{\prime}}}{=}0$ for $s{\neq}s^{\prime}$ in the degenerate subspace.
  * [59] Kané, G., Lazzeri, M. & Mauri, F.  Zener tunneling in the electrical transport of quasimetallic carbon nanotubes.  _Phys. Rev. B_ 86, 155433 (2012).
  * [60] Simoncelli, M., Marzari, N. & Cepellotti, A.  Generalization of fourier’s law into viscous heat equations.  _Phys. Rev. X_ 10, 011019 (2020).
  * [61] With “structurally stable glass” we mean a glass in which atoms vibrate around equilibrium positions, thus atomic diffusion is negligible. See e.g. refs. 135, 136, 137, 138, 139, 140, 141, 142 for a discussion of how network structure and chemical composition affect the structural stability of glasses and related phenomena.
  * [62] Anderson, P. W.  Absence of Diffusion in Certain Random Lattices.  _Phys. Rev._ 109, 1492–1505 (1958).
  * [63] Maradudin, A. A. & Vosko, S. H.  Symmetry properties of the normal vibrations of a crystal.  _Rev. Mod. Phys._ 40, 1–37 (1968).
  * [64] Fu, C. L. & Ho, K. M.  First-principles calculation of the equilibrium ground-state properties of transition metals: Applications to Nb and Mo.  _Phys. Rev. B_ 28, 5480–5486 (1983).
  * [65] Methfessel, M. & Paxton, A. T.  High-precision sampling for brillouin-zone integration in metals.  _Phys. Rev. B_ 40, 3616–3621 (1989).
  * [66] Marzari, N., Vanderbilt, D., De Vita, A. & Payne, M. C.  Thermal contraction and disordering of the al(110) surface.  _Phys. Rev. Lett._ 82, 3296–3299 (1999).
  * [67] We employ the nomenclature used in statistics, which classifies as “light-tailed” (“fat-tailed”) all the distributions having tails smaller (larger) than the exponential distribution.
  * [68] Zhu, X. & Shao, C.  Effect of anharmonicity on the thermal conductivity of amorphous silica.  _Phys. Rev. B_ 106, 014305 (2022).
  * [69] Marzari, N.  Ab-initio molecular dynamics for metallic systems.  _PhD thesis, Cambridge University_ (1996).
  * [70] de Gironcoli, S.  Lattice dynamics of metals from density-functional perturbation theory.  _Phys. Rev. B_ 51, 6773–6776 (1995).
  * [71] Paulatto, L., Errea, I., Calandra, M. & Mauri, F.  First-principles calculations of phonon frequencies, lifetimes, and spectral functions from weak to strong anharmonicity: The example of palladium hydrides.  _Phys. Rev. B_ 91, 054304 (2015).
  * [72] Paulatto, L., Mauri, F. & Lazzeri, M.  Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene.  _Phys. Rev. B_ 87, 214303 (2013).
  * [73] Li, W., Carrete, J., Katcho, N. A. & Mingo, N.  ShengBTE: A solver of the Boltzmann transport equation for phonons.  _Comput. Phys. Commun._ 185, 1747 – 1758 (2014).
  * [74] Carrete, J. _et al._ almaBTE : A solver of the space–time dependent Boltzmann transport equation for phonons in structured materials.  _Comput. Phys. Commun._ 220, 351 – 362 (2017).
  * [75] Togo, A., Chaput, L. & Tanaka, I.  Distributions of phonon lifetimes in brillouin zones.  _Phys. Rev. B_ 91, 094306 (2015).
  * [76] Tamura, S.-i.  Isotope scattering of dispersive phonons in Ge.  _Phys. Rev. B_ 27, 858–866 (1983).
  * [77] Garg, J.  _Thermal conductivity from first-principles in bulk, disordered, and nanostructured materials_.  Ph.D. thesis, Massachusetts Institute of Technology (2011).
  * [78] Thébaud, S., Berlijn, T. & Lindsay, L.  Perturbation theory and thermal transport in mass-disordered alloys: Insights from green’s function methods.  _Phys. Rev. B_ 105, 134202 (2022).
  * [79] To support this expectation, we note that the sum over the modes appearing in Eq. (1) contains $(3N_{at})^{2}$ term, out of which only $3N_{at}$ are diagonal ($s=s^{\prime}$) and thus vanishing at $\bm{q}{=}\bm{0}$. This implies the difference between a conductivity calculation at $\bm{q}{=}\bm{0}$ only and one on the $n{\times}n{\times}n$ mesh goes to zero with a speed roughly proportional to the ratio between the number of diagonal elements and total number of elements, i.e. $\frac{3N_{at}}{(3N_{at})^{2}}=\frac{1}{3N_{at}}$.
  * [80] Using a denser $5{\times}5{\times}5$ mesh produces practically equivalent results, implying that in the 192-atom model a $3{\times}3{\times}3$ mesh is sufficient to obtain computational convergence.
  * [81] Here “very low temperature” must be interpreted keeping into account that this work is focused on the above-the-plateau temperature regime, i.e. $T>30K$.
  * [82] Ida, T., Ando, M. & Toraya, H.  Extended pseudo-Voigt function for approximating the Voigt profile.  _J. Appl. Crystallogr._ 33, 1311–1316 (2000).
  * [83] Wray, K. L. & Connolly, T. J.  Thermal conductivity of clear fused silica at high temperatures.  _J. Appl. Phys._ 30, 1702–1705 (1959).
  * [84] Kanamori, H., Fujii, N. & Mizutani, H.  Thermal diffusivity measurement of rock-forming minerals from 300 to 1100 K.  _J. Geophys. Res._ 73, 595–605 (1968).
  * [85] Touloukian, Y., Powell, R., Ho, C. & Klemens, P.  Thermophysical properties of matter - the TPRC data series. Volume 2. Thermal conductivity - nonmetallic solids.  Tech. Rep., Thermophysical and Electronic Properties Information Analysis Center (1971).
  * [86] Sergeev, O., Shashkov, A. & Umanskii, A.  Thermophysical properties of quartz glass.  _J. Eng. Phys._ 43, 1375–1383 (1982).
  * [87] Cahill, D. G.  Thermal conductivity measurement from 30 to 750 K: the 3$\omega$ method.  _Rev. Sci. Instrum._ 61, 802–808 (1990).
  * [88] Lee, S.-M. & Cahill, D. G.  Heat transport in thin dielectric films.  _J. Appl. Phys._ 81, 2590–2595 (1997).
  * [89] Abdulagatov, I. _et al._ Thermal conductivity of fused quartz and quartz ceramic at high temperatures and high pressures.  _J. Phys. Chem. Solids_ 61, 779–787 (2000).
  * [90] Yang, S. T., Matthews, M. J., Elhadj, S., Draggoo, V. G. & Bisson, S. E.  Thermal transport in CO2 laser irradiated fused silica: In situ measurements and analysis.  _J. Appl. Phys._ 106, 103106 (2009).
  * [91] Heraeus.  Data sheet: Quartz glass for optics - data and properties (2010).  URL https://www.heraeus.com/media/media/hca/doc_hca/products_and_solutions_8/optics/Data_and_Properties_Optics_fused_silica_EN.pdf.
  * [92] Feldman, J. & Kluge, M. D.  Realistic model calculations based on the kubo theory for the thermal conductivity of amorphous insulators.  _Philos. Mag. B_ 71, 641–647 (1995).
  * [93] In all these cases we used a broadening $\eta=4$ cm-1 for the Voigt distribution, this value was determined from a convergence test analogous to that reported in the upper-left panel of Fig. 2.
  * [94] Obtained using the PBE exchange-correlation functional.
  * [95] Dragoni, D., Daff, T. D., Csányi, G. & Marzari, N.  Achieving dft accuracy with a machine-learning interatomic potential: Thermomechanics and defects in bcc ferromagnetic iron.  _Phys. Rev. Materials_ 2, 013808 (2018).
  * [96] Deringer, V. L.  Bernstein, N., Csányi, G., Mahmoud, C. B., Ceriotti, M., Wilson, M., Drabold, D. A., and Elliott, S. R. Origins of structural and electronic transitions in disordered silicon.  _Nature_ 589, 59–64 (2021).
  * [97] Simoncelli, M., Mauri, F. & Marzari, N.  Supporting data for “Thermal conductivity of glasses above the plateau: first-principles theory and applications”.  _Materials Cloud Archive_ (2022).  URL https://doi.org/10.24435/materialscloud:rw-rs.
  * [98] Bouchut, P., Decruppe, D. & Delrive, L.  Fused silica thermal conductivity dispersion at high temperature.  _J. Appl. Phys._ 96, 3221–3227 (2004).
  * [99] Lee, D. & Kingery, W.  Radiation energy transfer and thermal conductivity of ceramic oxides.  _J. Am. Ceram. Soc._ 43, 594–607 (1960).
  * [100] In this work the velocity operator is considered independent from anharmonicity and temperature, since we employ the standard approximation of considering the force constants to be independent from anharmonicity and temperature (for vitreous silica this approximation is realistic, as discussed in the Methods).
  * [101] This last condition implies that the Voigt distribution analytically reduces to the Lorentzian distribution (4).
  * [102] Hill, V. G. and Roy, Rustum.  Silica Structure Studies: V, The Variable Inversion in Cristobalite.  _J. Am. Ceram. Soc._ 41, 532–537 (1958).
  * [103] Eucken, A.  Über die temperaturabhängigkeit der wärmeleitfähigkeit fester nichtmetalle.  _Ann. Phys._ 339, 185–221 (1911).
  * [104] Birch, A. F. & Clark, H.  The thermal conductivity of rocks and its dependence upon temperature and composition.  _Am. J. Sci._ 238, 529–558 (1940).
  * [105] Knapp, W. J.  Thermal conductivity of nonmetallic single crystals.  _J. Am. Ceram. Soc._ 26, 48–55 (1943).
  * [106] Zeller, R. C. & Pohl, R. O.  Thermal Conductivity and Specific Heat of Noncrystalline Solids.  _Phys. Rev. B_ 4, 2029–2041 (1971).
  * [107] Lindsay, L.  First Principles Peierls-Boltzmann Phonon Thermal Transport: A Topical Review.  _Nanoscale and Microscale Thermophys. Eng._ 20, 67–84 (2016).
  * [108] Fugallo, G., Lazzeri, M., Paulatto, L. & Mauri, F.  Ab initio variational approach for evaluating lattice thermal conductivity.  _Phys. Rev. B_ 88, 045430 (2013).
  * [109] Richet, P., Bottinga, Y., Denielou, L., Petitet, J. & Tequi, C.  Thermodynamic properties of quartz, cristobalite and amorphous SiO2: drop calorimetry measurements between 1000 and 1800 K and a review from 0 to 2000 K.  _Geochim. Cosmochim. Acta_ 46, 2639–2658 (1982).
  * [110] Bianco, R., Errea, I., Paulatto, L., Calandra, M. & Mauri, F.  Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation.  _Phys. Rev. B_ 96, 014111 (2017).
  * [111] Aseginolaza, U. _et al._ Phonon collapse and second-order phase transition in thermoelectric snse.  _Phys. Rev. Lett._ 122, 075901 (2019).
  * [112] Horbach, J., Kob, W. & Binder, K.  Specific heat of amorphous silica within the harmonic approximation.  _J. Phys. Chem. B_ 103, 4104–4108 (1999).
  * [113] Brown, W. M., Wang, P., Plimpton, S. J. & Tharrington, A. N.  Implementing molecular dynamics on hybrid high performance computers–short range forces.  _Comput. Phys. Commun._ 182, 898–911 (2011).
  * [114] Csányi, G. _et al._ Expressive programming for computational physics in fortran 95+.  _IoP Comput. Phys. Newsletter_ Spring 2007 (2007).
  * [115] Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ 104, 136403 (2010).
  * [116] Kermode, J. R.  f90wrap: an automated tool for constructing deep python interfaces to modern fortran codes.  _J. Phys. Condens. Matter_ (2020).
  * [117] Eriksson, F., Fransson, E. & Erhart, P.  The hiphive package for the extraction of high-order force constants by machine learning.  _Adv. Theory Simul._ 2, 1800184 (2019).
  * [118] Giannozzi, P. _et al._ Advanced capabilities for materials modelling with Quantum ESPRESSO.  _J. Phys. Condens. Matter_ 29, 465901 (2017).
  * [119] Grimme, S.  Semiempirical GGA-type density functional constructed with a long-range dispersion correction.  _J. Comput. Chem_ 27, 1787–1799 (2006).
  * [120] Hay, H., Ferlat, G., Casula, M., Seitsonen, A. P. & Mauri, F.  Dispersion effects in ${\mathrm{SiO}}_{2}$ polymorphs: An ab initio study.  _Phys. Rev. B_ 92, 144111 (2015).
  * [121] Baroni, S., de Gironcoli, S., Dal Corso, A. & Giannozzi, P.  Phonons and related crystal properties from density-functional perturbation theory.  _Rev. Mod. Phys._ 73, 515–562 (2001).
  * [122] Prandini, G., Marrazzo, A., Castelli, I. E., Mounet, N. & Marzari, N.  Precision and efficiency in solid-state pseudopotential calculations.  _Npj Comput. Mater._ 4, 1–13 (2018).
  * [123] Lejaeghere, K. _et al._ Reproducibility in density functional theory calculations of solids.  _Science_ 351, aad3000 (2016).
  * [124] Brueckner, R.  Properties and structure of vitreous silica. I.  _J. Non-Cryst. Solids_ 5, 123–175 (1970).
  * [125] Pluth, J. J., Smith, J. V. & Faber, J.  Crystal structure of low cristobalite at 10, 293, and 473 K: Variation of framework geometry with temperature.  _J. Appl. Phys._ 57, 1045–1049 (1985).
  * [126] Downs, R. T. & Palmer, D.  The pressure behavior of $\alpha$ cristobalite.  _Am. Mineral._ 79, 9–14 (1994).
  * [127] Tucker, M., Keen, D. & Dove, M.  A detailed structural characterization of quartz on heating through the $\alpha$–$\beta$ phase transition.  _Mineral. Mag._ 65, 489–507 (2001).
  * [128] Rumble, J.  _CRC handbook of chemistry and physics_ , vol. 100 (CRC press, 2014).
  * [129] Deschamps, T., Margueritat, J., Martinet, C., Mermet, A. & Champagnon, B.  Elastic moduli of permanently densified silica glasses.  _Scientific reports_ 4, 1–7 (2014).
  * [130] Kondo, K.-i., Iio, S. & Sawaoka, A.  Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 gpa.  _J. Appl. Phys._ 52, 2826–2831 (1981).
  * [131] Le Roux, S. & Jund, P.  Ring statistics analysis of topological networks: New approach and application to amorphous GeS2 and SiO2 systems.  _Comput. Mater. Sci._ 49, 70–83 (2010).
  * [132] Silvestrelli, P. L., Marzari, N., Vanderbilt, D. & Parrinello, M.  Maximally-localized wannier functions for disordered systems: Application to amorphous silicon.  _Solid State Commun._ 107, 7–11 (1998).
  * [133] Mizokami, K., Togo, A. & Tanaka, I.  Lattice thermal conductivities of two ${\mathrm{SiO}}_{2}$ polymorphs by first-principles calculations and the phonon Boltzmann transport equation.  _Phys. Rev. B_ 97, 224306 (2018).