MnLargeSymbols’164 MnLargeSymbols’171

#  Vibrational and thermal properties of amorphous alumina from first
principles

Angela F. Harper  Equally contributed  Current address: Fritz-Haber Institut
der Max Planck Gesellschaft, Berlin (DE)  Theory of Condensed Matter Group,
Cavendish Laboratory, University of Cambridge (UK)  Kamil Iwanowski  Equally
contributed  Theory of Condensed Matter Group, Cavendish Laboratory,
University of Cambridge (UK)  William C. Witt  Department of Materials Science
& Metallurgy, University of Cambridge, Cambridge, (UK)  Mike C. Payne  Theory
of Condensed Matter Group, Cavendish Laboratory, University of Cambridge (UK)
Michele Simoncelli<EMAIL_ADDRESS>Theory of Condensed Matter Group,
Cavendish Laboratory, University of Cambridge (UK)

###### Abstract

Amorphous alumina is employed ubiquitously as a high-dielectric-constant
material in electronics, and its thermal-transport properties are of key
relevance for heat management in electronic chips and devices. Experiments
show that the thermal conductivity of alumina depends significantly on the
synthesis process, indicating the need for a theoretical study to elucidate
the atomistic origin of these variations. Here we employ first-principles
simulations to characterize the atomistic structure, vibrational properties,
and thermal conductivity of alumina at densities ranging from 2.28 g/cm3 to
3.49 g/cm3. Moreover, using an interatomic potential trained on first-
principles data, we investigate how system size affects predictions of the
thermal conductivity, showing that simulations containing 120 atoms can
already reproduce the bulk limit of the conductivity. Finally, relying on the
recently developed Wigner formulation of thermal transport, we shed light on
the interplay between atomistic topological disorder and anharmonicity in the
context of heat conduction, showing that the former dominates over the latter
in determining the conductivity of alumina.

##  I Introduction

Alumina (Al2O3) is a deceptively complex material; there are at least nine
known metastable polymorphs [1], an amorphous phase which possesses its own
unique properties as a result of local disorder [2], and a number of predicted
non-stoichiometric structures [3]. Alumina has applications in catalysis [4,
5], energy-storage devices (e.g., Li-ion batteries [6], Li-metal anodes [7],
and solid state electrolytes [8]), as well as in heat-management technologies
[9] and coatings [10]. Moreover, thin film amorphous alumina (am-Al2O3) is
often employed in electronic devices, and its thermal properties play a
crucial role in determining both device efficiency and lifespan [9, 11]. The
structural properties of am-Al2O3 have become subject of extensive research,
since the atomic-layer deposition (ALD) technique allows production of alumina
films with controlled thickness and morphology. There are several open
fundamental questions on how the structure of am-Al2O3 affects its macroscopic
properties. Solid-state 27Al nuclear magnetic resonance has identified three
main aluminum coordination environments in am-Al2O3 [12, 13, 14, 15, 16],
which play a role in the electronic transport as evidenced by X-Ray absorption
spectroscopy and first principles electronic density of states calculations
[17, 18]. The literature on the relationships between structural, vibrational,
and thermal properties of am-Al2O3 at technologically relevant temperatures
(i.e. at $T>30$ K, in the so called above-the-plateau regime [19, 20, 21]) is
less developed. Recently, Li  _et al._ [3] used machine learned force-fields
with molecular dynamics (Green-Kubo (GKMD [22, 23, 24, 25, 26, 27, 28]) and
non-equilibrium (NEMD [29, 30, 31, 32]) approaches) to compute the room-
temperature thermal conductivity of atomistic models of am-AlOx; they studied
models with densities ranging from 2.6 to 3.3 g/cm3 and containing up to 528
atoms. While providing insights on the physics governing thermal transport at
room temperature, this work could not explore how the thermal conductivity
varies with temperature—in fact, GKMD and NEMD are both governed by classical
equipartition and thus cannot be used to investigate how the conductivity
changes as temperature decreases (i.e., when the specific heat deviates
significantly from the classical limit [33]).

Here we employ first-principles calculations alongside a MACE machine learning
potential (MLP) [34, 35] to characterize the structural, vibrational, and
thermal properties of am-Al2O3 across a wide range of densities
($2.28{\leq}\rho{\leq}3.49$ g/cm3). (The MACE architecture utilizes the Atomic
Cluster Expansion [36, 37, 38] and equivariant message passing.) In
particular, in Sec. II we investigate from first principles how the atomic
coordination topology of am-Al2O3 changes with density, and how these changes
affect the vibrational properties. The vibrational properties are then used in
Sec. III as inputs for the Wigner formulation of thermal transport [39, 40],
which can be used to calculate heat conduction in amorphous solids accounting
comprehensively for the effects of structural disorder, anharmonicity, and
quantum Bose-Einstein statistics of vibrations [21]. We employ the
convergence-acceleration protocol discussed in Ref. [21] to compute the bulk
limit of the thermal conductivity of am-Al2O3, highlighting the agreement
between our predictions and experiments at various densities and temperatures.
We describe how the interplay between structural disorder and anharmonicity
affects the thermal conductivity of am-Al2O3, demonstrating that disorder
practically determines thermal transport at all the densities and temperatures
studied ($50{\leq}T{\leq}700$ K). Finally, in Sec. IV we employ the MACE MLP
to study size effects, comparing the structural, vibrational, and thermal
properties of models of am-Al2O3 containing 120, 3240, and 7680 atoms. We
discuss how the convergence-acceleration protocol of Ref. [21] affects the
thermal-conductivity predictions in models having different sizes.

Figure 1: Coordination of oxygen and aluminum in am-Al2O3 at various
densities. The histograms show the probability of oxygen (top) or aluminum
(bottom) to have a certain coordination number in the disordered models
analyzed (different densities are distinguished with colors). Coordination
varies from 2 to 4 for O, and from 4 to 6 for Al.

##  II Structural and vibrational properties

The properties of am-Al2O3 strongly depend on the synthesis process—samples
with densities ranging from 2.1 g/cm3 to 3.6 g/cm3 have been synthesized and
discussed in the literature [41, 42]. To generate structures throughout this
range, we extracted 120-atom structures at 3.17 g/cm3, 3.30 g/cm3 and 3.49
g/cm3 from Ref. [12], and additionally generated two 120-atom am-Al2O3 models
at 2.28 g/cm3 and 2.98 g/cm3 using a melt-quench procedure [12]. Each
structure was generated using ab initio molecular dynamics (AIMD): initially
melted at a temperature of 4000 K, then quenched to a temperature of 300 K,
and finally equilibrated in the NVT ensemble (see Appendix A.1 for details).
In the following, we show that these structures have very different local
atomic environments and coordination topologies, and consequently different
vibrational and thermal properties.

###  II.1 Coordination Topology

With the goal of gaining insights on how the atomic bonding topology affects
the macroscopic conductivity, we start by characterizing the coordination
environments present in our am-Al2O3 models, showing in Fig. 1 the probability
distribution of finding $x$-fold-coordinated oxygen ($\rm O_{x}$) or $y$-fold-
coordinated aluminum ($\rm Al_{y}$) at each density. Six different
coordination environments are present in our models—$\rm O_{x}$ with $x\in
2,3,4$, and $\rm Al_{y}$ with $y\in 4,5,6$—and at least five out of these six
coordination topologies coexist at every model density analyzed.

Figure 2:  Radial distribution function at various densities. Solid lines,
total RDF of the am-Al2O3 generated and studied in this work (a rigid shift of
two is used to distinguish the RDF of different models).

As shown in the top panel of Fig. 1, at increasing densities the concentration
of twofold-coordinated oxygens (O2) generally decreases in favor of fourfold-
coordinated (O4); changes in density have weaker effects on the concentration
of threefold-coordinated oxygens (O3). The bottom panel shows that increasing
density yields a reduction in the concentration of Al4 environments, which is
compensated by an increase of Al6. In contrast, the Al5 environment displays
weaker variations with density; we note that Al5 are characteristic of the
amorphous phase of Al2O3 and absent in the crystalline phases [43]. Finally,
we highlight how the 2.28 g/cm3 model contains a high proportion of O2
bridging oxygens and a majority of Al4 environments, with no Al6. This is
consistent with the trend that at lower densities Al becomes primarily
tetrahedrally coordinated in Al4 environments [10].

From the coexistence of at least five different coordination environments in
all the am-Al2O3 models studied, we expect strong structural disorder. To
validate these expectations and characterize atomistic disorder, we plot in
Fig. 2 the total radial distribution function (RDF) for all AIMD models. We
highlight how the the total RDF converges to one over a distance shorter than
5 Å, indicating the presence of strong short-range disorder (i.e. within the
second coordination shell [44, 45]). The decomposition of the total RDF into
partial RDFs will be discussed later in Fig 11, where we use structures
containing up to 7680 atoms to show that the oscillations in the partial RDF
become very weak (negligible) at distances larger than $\sim~{}10$ Å. We note,
in passing, that this decay distance is shorter than the linear size of all
the 120-atoms structures studied here, and in agreement with that found in
previous work [46]. This suggests that atomistic models containing hundreds of
atoms are sufficiently large to capture the structural properties of strongly
disordered solids [21, 47]. We validate this expectation in Sec. IV using the
MACE MLP and atomistic models containing up to 7680 atoms.

###  II.2 Vibrational properties

The strong variations in the coordination topology observed in the different
structures of am-Al2O3 (Fig. 1) are intuitively expected to affect the atomic
vibrational excitations. Thus, in this section we investigate quantitatively
how the vibrational properties of am-Al2O3 change with the coordination
topology. We used density-functional perturbation theory (DFPT) [48] to
compute the energy $\hbar\omega_{\bm{q}s}$ of each vibrational mode $\bm{q}s$,
and the corresponding displacement pattern $\mathcal{E}_{\bm{q}s}^{b\alpha}$,
which describes how atom $b$ moves in direction $\alpha$ when the vibration
$\bm{q}s$ is excited. We label vibrational modes using both the wavevector
$\bm{q}$ and the mode index $s$ for the sake of generality, keeping in mind
that $\bm{q}$ is necessary only in crystals and in finite-size models of
disordered solids, while in “ideal glasses” (namely, an astronomically large
set of atoms whose arrangement lacks long-range order) one would obtain
$\bm{q}=\bm{0}$, thus $s$ would be sufficient to label the vibrations [21].
Then, we use these quantities to investigate if there is a relationship
between atomic displacements and coordination topology, computing the root
mean square displacement of every atom $b$ [49],

$\text{RMSD}(b,T){=}\sqrt{\frac{\hbar}{N_{c}m_{b}}\sum_{\bm{q}s}\frac{1}{\omega_{\bm{q}s}}\Big{[}\tfrac{1}{2}+\mathsfit{N}_{\bm{q}s}\Big{]}\sum_{\alpha}|\mathcal{E}_{\bm{q}s}^{b\alpha}|^{2}},$
(1)

where
$\mathsfit{N}_{\bm{q}s}{=}[{\exp{(\hbar\omega_{\bm{q}s}/k_{B}T)}{-}1}]^{-1}$
is the Bose-Einstein distribution at temperature $T$, $m_{b}$ is the mass of
atom $b$, and $N_{c}$ is a normalization factor111N is equal to the number of
$\bm{q}$ points used to sample the Brillouin Zone; in the case of an ideal
glass (astronomically large disordered simulation cell) the Brillouin Zone
reduces to the point $\bm{q}=\bm{0}$ only and $N=1$..

Figure 3: Average atomic root mean square displacement in am-Al2O3 at various
densities, calculated at $T=300K$ using Eq. (1). The average RMSD of oxygen
decreases as coordination increases; in contrast, aluminum displays a much
weaker anticorrelation between RMSD and coordination number in models with
$\rho\leq 3.17$ g/cm3, such an anticorrelation disappears in models with
$\rho\geq 3.30$ g/cm3. Error bars represent the standard deviation of the
average RMSD, and are not reported in the cases in which only a single
coordination environment was detected.

To see if a relationship between RMSD and coordination topology exists, we
compute the average ${\rm RMSD}(b,T)$ over atoms with equal coordination. Fig.
3 shows that the average RMSD of Ox atoms decreases as coordination increases
in all the atomistic models analyzed. The average RMSD of Aly is
anticorrelated with coordination in models with $\rho{\leq}3.17$ g/cm3, and
displays deviations from this anticorrelated trend in models with $\rho=3.30$
g/cm3 (where all Al atoms have a very similar average RMSD, regardless of
their coordination) and with $\rho=3.49$ g/cm3 (where Al6 has average RMSD
slightly larger than Al5). Interestingly, as density increases the average
RMSD of O4 approaches the RMSD of Al4—considering the significant differences
between the mass of oxygen and aluminum ($m_{Al}/m_{O}\approx 1.7$), we note
that this behavior departs from the trend RMSD $\propto\sqrt{{m_{b}}^{-1}}$
ubiquitously observed in solids with simple coordination topology. This shows
that the presence of complex coordination topologies can have non-trivial
effects on the vibrational properties.

Having investigated how the constraints imposed by the coordination topology
affect the average RMSDs, we now study how coordination influences the
vibrational frequencies. To this aim, we compute the vibrational density of
states (VDOS), $g(\omega){=}({\mathcal{V}N_{\rm
c}})^{-1}\sum_{\bm{q},s}\delta(\omega{-}\omega_{\bm{q}s})$ (here,
$\mathcal{V}$ is the volume of the cell used to simulate the disordered
system), we use the eigenvectors to decompose the VDOS into partial (single-
atom) contributions (PDOS), and finally we integrate the single-atom PDOS
using an indicator function that allows us to resolve different coordination
environments:

$g_{t_{x}}(\omega)=\frac{1}{\mathcal{V}N_{c}}\sum_{{\bm{q}s}}\delta(\omega-\omega_{\bm{q}s})\sum_{b,\alpha}|\mathcal{E}^{b\alpha}_{{\bm{q}s}}|^{2}\delta_{b,t_{x}},$
(2)

where $\delta_{b,t_{x}}$ is indicator function equal to one if the atom $b$ is
of type $t$ and has coordination $x$ (e.g., $t_{x}=O_{3}$ for threefold-
coordinated oxygen), and zero otherwise. Thus, $g_{t_{x}}(\omega)$ allows us
to resolve how the coordination topology affects the VDOS.

Figure 4:  Vibrational density of states of am-Al2O3 at various densities. The
total VDOS is solid black. The colored solid lines distinguish coordination
environments for Al atoms: green is Al4, red is Al5, and yellow is Al6. Dashed
colored lines are used for coordination environments of O atoms: cyan for O2,
blue for O3, and purple for O4.

Results for $g_{t_{x}}(\omega)$ are shown in Fig. 4. We highlight how the PDOS
for different coordination environments have similar shapes at different
densities. In particular, different coordination environments for oxygen have
fingerprints in the coordination-resolved PDOS: O3 is characterized by a
bimodal PDOS with a local minimum around 500 cm-1, while O2 and O4 have a
unimodal PDOS with peaks at about 150 and 500 cm-1, respectively. In contrast,
the coordination-resolved PDOS for Al4, Al5, and Al6 are similar in shape. We
note that this is in sharp contrast to the crystalline $\theta$-Al2O3 phase,
which contains only Al4 tetrahedra and Al6 octahedra and displays a clear
distinction between high-frequency breathing modes associated with Al4
environments, and low-frequency bending modes associated with Al6 environments
[51].

Varying the density of am-Al2O3 causes a variation in the relative magnitude
of the Alx PDOS: at 2.28 g/cm3, the Al4 PDOS is always larger in magnitude
than Al5; with increasing density, the relative magnitude of the PDOS of Al4
and Al5 progressively reverses, with the highest-density 3.49 g/cm3 model
featuring a PDOS for Al5 always larger than that for Al4.

We note that the O3 and Al5 coordination environments have a significant PDOS
in all the models analyzed, regardless of the models’ density. In contrast,
the PDOS of O2, O4, Al4, Al6, display much stronger variations with
density—increasing density yields a progressive suppression of the O2 and Al4
environments, compensated by the progressive appearance of O4 and Al6
environments.

##  III Thermal properties

###  III.1 Thermal conductivity of glasses

In this section we summarize the salient features of the Wigner formulation of
thermal transport [39], which describes heat conduction accounting for the
interplay between structural disorder, anharmonicity, and quantum Bose-
Einstein statistics of atomic vibrations. This allows us to describe the
thermal conductivity of solids ranging from crystals [52, 40, 53] to glasses
[21, 54].

When applied to atomistic models of amorphous systems, the Wigner formulation
of thermal transport yields the following ‘rWTE’ conductivity expression[21]:

$\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}_{[\eta;\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}]}(\omega_{\bm{q}s}-\omega_{\bm{q}s^{\prime}})\;,\end{split}$
(3)

where $\hbar\omega_{\bm{q}s}$ denotes the harmonic energy of the vibration
$\bm{q}s$, and $\hbar\Gamma_{\bm{q}s}$ its linewidth (broadening of the
harmonic energy level due to anharmonicity [55, 56, 57, 58, 59, 60, 61] and
isotopic-mass disorder [62]);
$\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}$
is the square modulus of the velocity operator [40] between eigenstates $s$
and $s^{\prime}$ at fixed $\bm{q}$ ($\alpha$ denotes a Cartesian direction,
and since amorphous solids are in general isotropic, the scalar conductivity
(3) is computed as the average trace of the tensor $\kappa^{\alpha\beta}$
[21]);

$C_{\bm{q}s}{=}C[\omega_{\bm{q}s},T]{=}\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{)}$
(4)

is the quantum specific heat of the mode $\bm{q}s$; $\mathcal{V}$, $N_{\rm
c}$, and $\mathsfit{N}_{\bm{q}s}$ are, in order, the simulation’s cell volume,
normalization factor and Bose-Einstein distribution already discussed in Sec.
II.2. Finally,
$\mathcal{F}_{[\eta;\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}]}(\omega_{\bm{q}s}-\omega_{\bm{q}s^{\prime}})$
is a Voigt distribution, i.e. a two-parameter distribution that reduces to a
Lorentzian with full width at half maximum (FWHM)
$\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}$ when
$\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}\gg\eta$, and to a Gaussian
representation of the Dirac delta with variance $\eta^{2}\pi/2$ in the
opposite limit ($\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}\ll\eta$).

As discussed in detail in Ref. [21], Eq. (3) describes thermal transport as
originating from couplings between vibrations that have an energy difference
smaller than the broadening of the Voigt profile. Such a broadening is
determined by $\eta$ in the low-temperature (harmonic) limit where
anharmonicity phases out ($\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}\to
0\;\forall\;\bm{q},s$). Setting $\eta$ to a value slightly larger than the
average energy-level spacing $\hbar\Delta\omega_{\rm
avg}=\frac{\hbar\omega_{\rm max}}{3\cdot N_{\rm at}}$ ($\hbar\omega_{\rm max}$
is the maximum vibrational energy, $3N_{\rm at}$ is the number of vibrational
energy levels, i.e. three times the number of atoms in the system’s reference
cell) accounts for the physical property that heat transfer via a wave-like
tunneling between neighboring (quasi-degenerate) vibrational eigenstates can
occur 222A necessary condition for this to happen is to have vibrations that
are not localized in the Anderson sense [99], i.e. to have non-zero velocity
operator elements in Eq. (3). even in the limit of vanishing anharmonicity,
implying that in such a limit Eq. (3) reduces to the harmonic Allen-Feldman
(AF) thermal conductivity for glasses [20, 64]. In contrast, at temperatures
where the anharmonic linewidths are much larger than the computational
broadening $\eta$, the Voigt profile becomes a Lorentzian with FWHM
$\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}$, effectively reducing to the
anharmonic Wigner conductivity expression [40]. In practice, when applied to
finite-size atomistic models of amorphous solids at temperatures for which
$\hbar\Delta\omega_{\rm avg}>\Gamma_{\bm{q}s}$, the rWTE conductivity
expression (3) ensures that the low-temperature harmonic Allen-Feldman limit
is correctly described, and the effects of anharmonicity are considered only
when they are not spuriously altered by finite-size effects [40].

In actual calculations, an amorphous solid is approximately described as a
crystal having a primitive cell containing a large but finite number of atoms
$N_{\rm at}$. Thus, the Brillouin Zone (BZ) corresponding to the (large)
finite-size model does not reduce to the aforementioned “ideal glass” limit
($\bm{q}{=}\bm{0}$ only), but has a (small) finite volume. Recent work [21]
has shown that when the lengthscale of structural disorder is shorter than the
size of the simulation cell, periodic boundary conditions and Fourier
interpolation over the small BZ of the glass can be exploited to improve the
sampling of the vibrational properties, and extrapolate the bulk limit of the
thermal conductivity. In practice, using the Fourier interpolation in a
disordered model corresponds to averaging over different possible boundary
conditions. For example, the vibrational properties computed at
$\bm{q}{=}\bm{0}$ only correspond to considering periodic boundary conditions
at the boundaries of the simulation box; calculations at
$\bm{q}{=}(1,0,0)\pi/L$ correspond to considering antiperiodic boundary
conditions along direction $x$ and periodic along the other directions
(Feldman  _et al._ [65]). Of course, increasing the size of the atomistic
model, one expects the choice of the boundary conditions to become less and
less relevant, and in practice this can be verified by comparing a calculation
at $\bm{q}{=}\bm{0}$ only with one performed relying on Fourier
interpolation—when differences are negligible, one can conclude that the
boundary conditions are irrelevant and the system is large enough to allow a
brute-force simulation of the bulk limit. In practice it is not always
possible to simulate this, due to exceedingly large computational cost. In the
next section we show that in finite-size models of strongly disordered systems
such as am-Al2O3, the bulk limit of the conductivity can be computed using the
rWTE (3), which ensures that: (i) heat transfer mediated by tunneling between
neighboring eigenstates can always occur, preserving a physical property that
would otherwise emerge only in the thermodynamic limit; (ii) the sampling of
the vibrational properties is improved by averaging over different boundary
conditions. These statements are validated in Sec. IV using the MACE MLP for a
brute-force calculation of the bulk limit in models containing up to 7680
atoms.

Finally, it is worth mentioning that the Wigner thermal conductivity
expression (3) can be derived also from a many-body Green-Kubo approach [66,
67], and such an expression has been recently employed [68, 69], in
combination with interatomic potentials, to describe the thermal properties of
several glasses 333More precisely , Refs. [68, 69] evaluated the bulk limit of
the conductivity relying on empirical interatomic potentials and atomistic
models containing thousands of atoms, i.e. having a size large enough to
achieve computational convergence by evaluating Eq. (3) at $\bm{q}=\bm{0}$
only and without relying on the Voigt regularization. See Ref. [21] for
details on the conditions under which evaluating Eq. (3) yields equivalent
results with or without relying on the Fourier interpolation and Voigt
regularization..

###  III.2 Numerical results

####  III.2.1 Thermal conductivity from first principles

In this section we evaluate the rWTE conductivity from first principles (see
Appendix A.1 for details) for each am-Al2O3 structure, in the temperature
range from 50 K to 700 K. We focus on this temperature range because it is the
most relevant for technological applications related to electronics [6, 7, 8],
and can be studied with the 120-atom models at our disposal [21] (this last
statement is discussed in detail and validated in Sec. IV, where we use a MLP
to study models containing up to 7680 atoms).

Figure 5: Wigner vs Allen-Feldman thermal conductivity of am-Al2O3.
Predictions from anharmonic rWTE are solid, from harmonic AF are dashed;
different colors distinguish different densities 2.28 g/cm3 (blue and cyan),
2.98 g/cm3 (dark and light green) and 3.49 g/cm3 (red and orange). The effects
of anharmonicity are overall small, and become weaker as density decreases.
Scatter points are experiments: red triangles are taken from Lee  _et al._
[71] (DC sputtering), green circles from Monachon and Weber [72] (ALD on Si
substrate), and the green square is from Gorham  _et al._ [73] (ALD on Si
substrate). Theory and experiments are in reasonably good agreement over the
entire temperature range.

We show in Fig. 5 the predictions for the thermal conductivity obtained using
the harmonic AF or anharmonic rWTE formulations (see Appendix C.1 for details
on the convergence test for the AF broadening parameter $\eta$ and for the
calculation of the anharmonic linewidths). We find that the effects of
anharmonicity are in general weak—AF and rWTE differ at most by 10$\%$ in the
highest-density (3.49 g/cm3) model. These differences become smaller as the
density decreases, and are practically invisible in the lowest-density (2.28
g/cm3) model. The good agreement between AF and rWTE shows that in am-Al2O3
the vibrations’ damping due to disorder is strong enough to dominate over
anharmonicity. In addition, we highlight how both AF and rWTE predict an
increasing-up-to-saturation trend for the temperature-conductivity curve; such
a saturating trend is inherited from that of the specific heat (more on this
later in Sec. III.2.2). Our calculations shed light on the thermal properties
of am-Al2O3 below room temperature, where the quantum Bose-Einstein statistics
of vibrations has a major effect on thermal transport. This is a step forward
compared to previous studies based on molecular dynamics (MD) [3], which were
governed by classical equipartition and thus had to be limited to the high-
temperature regime [33] (i.e. at temperature large enough to have a quantum
specific heat effectively very close to the constant classical limit).

In Fig. 5 we compare our calculations with the experimental measurements from
Refs. [72, 73] (ALD samples grown on Si substrate and having density 3.0
g/cm3) and Ref. [71] (DC-sputtered samples with density 3.51 g/cm3). Thermal
conductivity experiments have uncertainties that depend on many factors,
including, e.g., the sample’s size, shape, and measurement method [74, 75,
76]; as shown by the error bars from Refs. [72, 73], these uncertainties are
typically around $\sim 10-20\%$. We see that our predictions at different
densities are overall compatible (within the error bars or within 10$\%$) with
the corresponding experiments.

We now turn our attention to the dependence of the conductivity on density. In
Fig. 6 we plot the room-temperature rWTE conductivity as a function of density
(we note from Fig. 5 that at 300 K the rWTE and AF conductivities are
practically indistinguishable across the entire density range analyzed),
finding an approximately linear increase of the conductivity with density
($\kappa(\rho)_{300K}\approx a\cdot\rho+b$, where
$a{=}0.637{\pm}0.004\tfrac{\rm W\cdot cm^{3}}{\rm m\cdot K\cdot g}$,
$b{=}{-}0.495{\pm}0.014\tfrac{\rm W}{\rm m\cdot K}$). We compare our
predictions with the experiments by Gorham  _et al._ [73], who characterized
the room-temperature thermal conductivity of ALD am-Al2O3 grown on Si
substrate 444We note, in passing, that Refs. [72, 78, 73] did not observe a
significant dependence of the thermal conductivity from the substrate on which
the am-Al2O3 sample was grown. at densities ranging from 2.66 g/cm3 to 3.12
g/cm3. To give an idea about the conductivity variability found when comparing
independent experiments, we also report the room-temperature conductivities
extracted from the dataset of Lee  _et al._ [71] (1995, DC-sputtered) and Lee
_et al._ [78] (ALD on Si, 2017) already presented in Fig. 5. We note that
there is some variance between the conductivity observed in independent
experiments, still they show a conductivity that overall increases with
density, and such a trend is reproduced in our calculations.

Figure 6:  Thermal conductivity as a function of density. The solid green line
is the theoretical rWTE conductivity computed at 300 K. Empty scatter points
are experiments performed at 300 K by Gorham  _et al._ [73] in 2014 (green
squares, ALD samples grown on Si), by Lee  _et al._ [78] in 2017 (blue
diamond, ALD samples grown on Si), and by Lee  _et al._ [71] in 1995 (red
triangle, DC sputtering).

####  III.2.2 Thermal diffusivity and effects of anharmonicity

To gain microscopic, fundamental insights on why anharmonicity has negligible
effects on the thermal conductivity of am-Al2O3, it is useful to resolve the
amount of heat carried by each atomic vibration and its diffusion rate. This
information can be obtained by recasting the rWTE conductivity expression (3)
as [21]

$\kappa(T)=\int_{0}^{\omega_{\rm
max}}g(\omega)C(\omega,T)D(\omega,T)d\omega\;,$ (5)

where $\omega_{\rm max}$ is the maximum vibrational frequency of the system,
$g(\omega)$ is the VDOS discussed in Sec. II.2), $C(\omega,T)$ is the specific
heat for a vibration with frequency $\omega$ at temperature $T$ (see Eq. (4)),
and $D(\omega,T)$ is the rWTE diffusivity [21],

$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}),$ (6)
$\displaystyle D_{\bm{q}s}{=}$
$\displaystyle{\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}$
(7) $\displaystyle\hskip
14.22636pt\times\pi\mathcal{F}_{[\eta;\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}]}(\omega_{\bm{q}s}-\omega_{\bm{q}s^{\prime}})\;.$

The expression for $D_{\bm{q}s}$ is determined by factorizing the specific
heat $C_{\bm{q}s}$ in Eq. (3) and by the requirement that in the 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 [40] (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}}}$). In the
harmonic AF limit
$\eta{\gg}\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}{\to}0$, thus the Voigt
distribution in Eq. (7) reduces to the Gaussian representation of the Dirac
$\delta$, accounting only for couplings between quasi-degenerate vibrational
eigenstates and effectively reducing to the temperature-independent AF
diffusivity [20, 64]. We note that in the context of Eq. (5) the only
difference between AF and rWTE originates from the diffusivity $D(\omega,T)$,
implying that differences between AF and rWTE conductivities derive
exclusively from differences between the AF and rWTE diffusivities.

Figure 7:  Velocity operator of am-Al2O3 as a function of vibrational energy
differences and averages, for the models having density 2.28, 2.98, and 3.49
g/cm3, computed from Eq. (8). We see that the average velocity operators for
all densities are relatively unchanged at increasing differences,
$\omega_{d}$. These small differences imply that the effects of anharmonicity
on the conductivity are negligible (see text).

Eqs. (6,7) allow us to obtain additional, microscopic information on why the
harmonic AF and anharmonic rWTE conductivities in Fig. 5 display small
(negligible) differences in the high-temperature limit
($C_{\bm{q}s}=C(\omega_{\bm{q}s},T){\to}k_{B}$ and $\Gamma_{\bm{q}s}{\gg}\eta$
$\forall\;\bm{q}s$). Specifically, the Voigt distribution in Eqs. (6,7)
implies that in the high-temperature limit the rWTE diffusivity is a
Lorentzian-weighted average of velocity-operator elements
$\lVert\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}\lVert^{2}\\!$, and in practice
such average is mainly determined by elements satisfying
$|\omega_{\bm{q}s}{-}\omega_{\bm{q}s^{\prime}}|{<}[\Gamma_{\bm{q}s}{+}\Gamma_{\bm{q}s^{\prime}}]$;
increasing temperature implies an increase of the linewidths, and therefore
velocity-operator elements with increasingly larger frequency difference
contribute to such average. Therefore, the trend of the velocity-operator
elements as a function of the energy difference
$\hbar|\omega_{\bm{q}s}{-}\omega_{\bm{q}s^{\prime}}|$ determines the trend of
the diffusivity as a function of temperature: elements increasing (decreasing)
with frequency difference imply a conductivity increasing (decreasing) with
temperature. In contrast, in the harmonic AF limit, the diffusivity is always
determined by quasi-degenerate velocity-operator elements
($\hbar|\omega_{\bm{q}s}{-}\omega_{\bm{q}s^{\prime}}|\to 0$) and therefore it
does not depend on temperature. Applying this reasoning to am-Al2O3, where the
rWTE conductivity saturates to a temperature-independent value at high
temperature, we expect the velocity-operator elements to be roughly
independent from the frequency difference. To validate this reasoning, we show
in Fig. 7 the velocity operator as a function of the energy average
$\hbar\omega_{a}=\hbar[\omega_{\bm{q}s}{+}\omega_{\bm{q}s^{\prime}}]/2$ and
difference
$\hbar\omega_{d}=\hbar|\omega_{\bm{q}s}{-}\omega_{\bm{q}s^{\prime}}|$ of the
modes coupled [21]:

$\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}$
(8)

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

$\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}\left(\frac{\omega(\bm{q})_{s}+\omega(\bm{q})_{s^{\prime}}}{2}-\omega_{\rm
a}\right)\\\
\times&{\delta}\big{(}(\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}})-\omega_{\rm
d}\big{)}\;.\end{split}$ (9)

The plots confirm our expectations, i.e., in the am-Al2O3 models studied the
saturating trend of the rWTE conductivity, and the negligible effect of
anharmonicity, derive from having microscopic velocity-operator elements that
do not vary appreciably with $\omega_{\rm d}$ across the range of densities
studied. Keeping the negligible differences between AF and rWTE diffusivity in
mind, in the following we focus on the temperature-independent AF limit of the
diffusivity $D(\omega)$ to simplify the discussion (i.e. Eqs. (6,7) evaluated
with $\Gamma_{\bm{q}s}=0\;\forall\bm{q}s$ and $\eta$ determined from the
convergence plateau as discussed in the Appendix C.1).

Eq. (5) shows that the contributions to heat transport of atomic vibrational
modes with frequency $\omega$ is determined by their density of states
$g(\omega)$, amount of heat carried $C(\omega,T)$ and rate of diffusion
$D(\omega)$. In previous sections we showed that: (i) VDOS increases with
density as discussed in Fig. 4; (ii) conductivity increases (linearly) with
density as evidenced in Fig. 6. Therefore, it is natural to ask to what extent
the conductivity increase observed in Fig. 6 derives from the increase in the
VDOS with density, and how density affects the diffusivity of vibrations.

Figure 8: AF diffusivity of am-Al2O3 at various densities, computed using Eq.
6 in the AF limit. The diffusivity tends to increase with density, especially
at low frequencies (0 - 100 cm-1 range) and between 400 and 600 ${\rm
cm}^{-1}$. Inset, quantum specific heat as a function of temperature.

Analyzing the frequency-resolved AF diffusivity $D(\omega)$ reported in Fig. 8
for all our models of am-Al2O3 allows us to address these questions. We see
that an overall increase of diffusivity with density is visible when comparing
low-density ($\rho=$2.28 g/cm3), medium-density ($\rho=$2.98 g/cm3) and high-
density ($\rho\geq 3.17$ g/cm3) models, especially at low frequencies (from 0
to 100 cm-1). More precisely, comparing the diffusivity of the highest-density
$\rho=3.49$ g/cm3 model with that of the lowest-density $\rho=2.28$ g/cm3
model, we find that the highest-density model has a diffusivity that is a
factor $\sim$2 larger than that of the lowest-density model in the low-
frequency region (0-100 cm-1) and a up to a factor 1.5 larger in the region
between 400 and 600 cm-1. The inset of Fig. 8 shows the quantum specific heat
$C(\omega,T)$ at various temperatures, whose frequency dependence at fixed
temperature $T$ is indicative of the portion of vibrational spectrum that
significantly contributes to heat transport at that temperature.

We highlight how the low-frequency vibrational modes which have density-
dependent diffusivity are significantly populated at all the temperatures
considered, implying that the increase of the thermal conductivity with
density observed in Fig. 6 receives contributions also from increases in the
diffusivity. We also note that the saturation of the specific heat shown in
the inset of Fig. 8 drives the saturation of the AF thermal conductivity at
high temperature (Fig. 5). Given that the effects of anharmonicity are
unimportant for thermal transport in am-Al2O3, the saturation of the rWTE
conductivity has to be attributed to the saturation of the specific heat.
Finally, we highlight how the classical limit (dashed line in the inset of
Fig. 8) is not yet reached even at temperatures as high as 1200 K; this
underscores the importance of correctly accounting for the quantum Bose-
Einstein statistics of vibrations to describe thermal transport in am-Al2O3 at
technologically relevant temperatures.

Figure 9: Microscopic mechanisms underlying conductivity increase with
density. Solid lines are exact AF conductivity calculation for the highest-
density 3.49 g/cm3 model (orange, $g_{3.49}D_{3.49}$) and lowest-density 2.28
g/cm3 model (cyan, $g_{2.28}D_{2.28}$). The dashed lines show results of
artificial conductivity calculations, performed using Eq. (5) with the VDOS of
the highest-density model and the diffusivity of the lowest-density model
(purple, $g_{3.49}D_{2.28}$), or the VDOS of the lowest-density model and the
diffusivity of the highest-density model (grey, $g_{2.28}D_{3.49}$). The
artificial calculations yield conductivities lying approximately halfway
between the exact limits, indicating that the increase of conductivity with
density is determined, in similar proportion, by both an increases in vDOS and
by an increase in diffusivity.

In summary, we have found that both VDOS (Fig. 4) and diffusivity (Fig. 8)
increase with density. In order to estimate how much the increase of
conductivity with density shown in Figs. 5,6 depends on increase of the VDOS
and how much on the increase in diffusivity, we computed the conductivity
artificially combining the VDOS and diffusivity of the highest and lowest-
density models in the following ways: (i) using the VDOS of the 3.49 g/cm3
model and the AF diffusivity of the 2.28 g/cm3 model; (ii) using the VDOS of
the 2.28 g/cm3 model and the AF diffusivity of the 3.49 g/cm3 model. The
comparison between these artificial conductivities and the exact ones (taken
from Fig. 5) are reported in Fig. 9, and show that these two artificial
calculations yield a conductivity that lies approximately halfway between the
actual conductivities of the lowest-density and highest-density models over a
wide temperature range, demonstrating that variations of the thermal
conductivity with density are determined by both variations in the VDOS and in
the diffusivity.

##  IV Size effects

The AF and rWTE conductivities presented in the previous section were
evaluated following the protocol presented in Ref. [21], which showed that the
thermal conductivity of strongly disordered solids such as vitreous silica can
be accurately reproduced using models containing ${\sim}100$ atoms, a result
supported by other recent studies [79, 47]. In this section we study how the
size of the atomistic model affects the theoretical predictions for the
conductivity of am-Al2O3. To overcome the limitations posed by the
computational cost of first-principles calculations, we generated a MACE MLP
[34, 35] for am-Al2O3 using the first-principles dataset released by Li et al.
[3]; computational details are reported in Appendix A.2. Then, we used this
potential to produce am-Al2O3 models containing 120, 3240, and 7680 atoms
before finally computing the vibrational properties of these models, as well
as the AF and rWTE conductivities.

###  IV.1 Model generation via melt-quench

Figure 10:  Coordination histogram for AIMD and MACE models of am-Al2O3. The
histograms show the proportion of oxygen (top) and aluminum (bottom) with a
certain coordination number. The atomistic structures above the plot show
differences in the linear size of the models studied.

We used MACE to generate large models of am-Al2O3. Specifically, we performed
molecular-dynamics melt-quench simulations, using the same protocol employed
in [12] to generate the AIMD 120-atom models discussed in the previous
sections. Starting from the AIMD 120-atom model with density 2.98 g/cm3
(hereafter referred to as AIMD120), we generated $3\times 3\times 3$ and
$4\times 4\times 4$ supercells, containing 3240 and 7680 atoms, respectively.
These supercells were used as the initial configuration for a melt-quench
simulation: they were first heated from 0 to 4000 K in 10 ps; then melted for
10 ps at 4000 K to ensure randomness of the structure; and then quenched to
300 K over 10 ps. After the melt-quench simulation, each structure was
equilibrated in an NVT ensemble at 300 K for 10 ps. Finally, the resulting
models were relaxed to a pressure lower than $0.001$ kBar and to interatomic
forces lower than 1 meV/Å (without imposing any constraint on the geometry of
the simulation box). After this final relaxation, both the 3240- and 7680-atom
models (hereafter referred to as MACE3240 and MACE7680) displayed a density of
2.92 g/cm3. We also relaxed with MACE the AIMD120 model; using the same
threshold mentioned above and without any constraint on the geometry of the
box, obtaining a ‘MACE120’ model with density 2.88 g/cm3. The densities of
structures obtained from first principles and MLP are in acceptable agreement,
being within $3\%$. In the next section we validate the capability of MACE to
reproduce within few-percent accuracy the first-principles results for the
structural, vibrational, and thermal properties of am-Al2O3. Then, using the
120-, 3240-, and 7680-atom MACE models, we investigate how these properties
depend on the size of the model.

###  IV.2 Structural properties

Figure 11:  Radial distribution function for models of am-Al2O3. Total and
partial RDF of the 120 atom Al2O3 model generated using AIMD at 2.98 g/cm3
compared with the two models generated using the MACE potential. The RDFs for
the 120-atom models only extends to $r{=}10.9$Å, a distance equal to the
linear size of these models.

We start by comparing coordination histograms of MACE structures with the
AIMD120 structure. Fig. 10 shows that both Al and O coordination distributions
are in good agreement between MACE120 and AIMD120. Importantly, these
distributions are practically independent from the size of the model, since
MACE120, MACE3240, and MACE7680 show extremely similar distributions.

Next, to resolve possible more subtle difference between the various
structures, we compare in Fig. 11 the total and partial RDFs. The total RDF of
MACE120, MACE3240, and MACE7680 are practically indistinguishable, and they
are in remarkable agreement with the RDF of AIMD120. Similarly, all the
partial RDFs (Al-Al, Al-O, and O-O) are in remarkable agreement between all
the four aforementioned models. We highlight how oscillations in the partial
RDFs become very weak (negligible) at distances larger than the linear size of
our first-principles 120-atom models ($\sim 11$ Å, where the solid light-green
and dashed dark-green lines stop). This suggests that atomistic models with
linear size of $\sim 11$ Å are sufficiently large to capture the most
important features of structural disorder in am-Al2O3. Importantly, these
tests also validate the capability of the MACE MLP to describe the structural
properties of am-Al2O3 with first-principles accuracy. Therefore, we continue
our validation tests for MACE, discussing in the next section its capability
to reproduce the first-principles vibrational properties.

###  IV.3 Vibrational properties

Figure 12:  Vibrational density of states for AIMD and MACE am-Al2O3 models.
The total VDOS is solid black. The colored solid lines distinguish
coordination environments for Al atoms: green is Al4, red is Al5, and yellow
is Al6. Dashed colored lines are used for coordination environments of O
atoms: cyan for O2, blue for O3, and purple for O4. From top to bottom:
MACE7680, MACE3240 (both computed at $\bm{q}{=}0$ only); MACE120 and AIMD120
(both computed on a 5x5x5 $\bm{q}$ mesh).

In Fig. 12 we test the capability of MACE to reproduce the vibrational
properties previously computed from first principles, and we also check how
these depend on the size of the model. Starting from the bottom, we see that
the AIMD120 and MACE120 models show very similar full vDOS (black line). Also,
the decomposition of the vDOS into PDOS contributions from different
coordination environments (colored lines) is compatible between MACE120 and
AIMD120.

Turning our attention to size effects on vibrational properties, we see that
the large MACE3240 and MACE7680 models display a vDOS similar in magnitude and
shape to the 120-atom models but differing in the roughness, i.e., larger
models have vDOS, which is a smoothed version of the vDOS computed for small
MACE120 and AIMD120 models. This is also true for the PDOS, where general
shape is maintained but roughness is reduced as the model size increases. We
conclude by noting that, in all these models, the vibrational frequencies are
non-negative, confirming that our models were correctly relaxed to an energy
minimum and therefore are structurally stable.

###  IV.4 Thermal conductivity

In this section we use the vibrational properties of the MACE models to
evaluate the Wigner thermal conductivity (3), and thus test how it converges
with respect to system size.

Figure 13:  rWTE and bWTE thermal conductivities of am-Al2O3 models at density
$\sim$3.0 g/cm3. The solid lines (scatter points) are bWTE (rWTE)
conductivities: MACE7680 is purple, MACE3240 is orange, MACE120 is blue and
AIMD is green.

In Fig. 13 we compare the bare WTE conductivities (i.e., Eq. (3) computed at
$\bm{q}=\bm{0}$ with $\eta=0$, hereafter referred to as ‘bWTE’) against the
regularized rWTE conductivity (computed relying on the $\bm{q}$-mesh
interpolation and using $\eta$ determined from the beginning of the
convergence plateau, see Figs. 15,17 in Appendix for details). We see that for
the large MACE7680 and MACE3240 models the bWTE and rWTE conductivities are
practically indistinguishable, confirming that these models are sufficiently
large to realistically represent a bulk system at temperatures above 50 K,
thus the convergence-acceleration protocol based on the Voigt regularization
protocol has negligible effect on them. Importantly, the rWTE conductivity for
the small AIMD120 and MACE120 models computed on converged 5x5x5 $\bm{q}$-mesh
and using $\eta=6$ cm-1 yields a conductivity compatible 555The largest
difference between the converged conductivities is of about 3.3$\%$ ($\approx$
0.057 W/mK) observed at 700 K between the rWTE computed over a 5x5x5
$\bm{q}$-mesh for MACE120, and rWTE computed over a 3x3x3 $\bm{q}$-mesh for
MACE7680. with that of the large MACE7680 and MACE3240 models. We highlight
that it is crucial to employ the rWTE to extrapolate the bulk limit of the
conductivity from small models—Fig. 13 shows that in AIMD120 and MACE120 the
bWTE underestimates the conductivity by approximately $20\%$.

It is worth commenting on how the rWTE protocol extrapolates successfully to
the bulk limit for the conductivity. As anticipated in Sec. III.1 and in Ref.
[21], the rWTE protocol enforces the physical property that couplings between
different modes can always occur in a truly disordered system, and exploits
$\bm{q}$-mesh interpolation to average the vibrational properties over many
possible different boundary conditions—this last operation is analogous to the
averaging over periodic and anti-periodic boundary conditions employed by
Feldman  _et al._ [65] to accelerate computational convergence. To better
understand the effect of the $\bm{q}$-mesh interpolation, we show in Fig. 18
that the velocity-operator as a function of frequency for MACE120 on a
$5\times 5\times 5$ $\bm{q}$ mesh, approaches the velocity operator as a
function of frequency computed at $\bm{q}=\bm{0}$ for the larger MACE7680
model. We note that to represent the velocity operator as a function of
frequency using Eq. (8) in a finite-size model, the Dirac delta must be
broadened with a finite-width (Gaussian) distribution having width of the
order of few energy-level spacings $\hbar\Delta\omega_{\rm avg}$, i.e. a
broadening comparable to the value $\eta$ used in the evaluation of the AF
conductivity. In the context of Eq. (3), employing a broadening of the order
of $\hbar\Delta\omega_{\rm avg}$ ensures in practice that such an expression
behaves analogously to that describing a strongly disordered system in the
bulk limit. The vibrational eigenstates repel each other significantly [81];
therefore, vibrational modes are dense but never degenerate and, for
arbitrarily small values for the broadening $\eta$, couplings between
neighboring vibrational eigenstates are allowed. The model size necessary to
achieve computational convergence with the rWTE protocol depends on the degree
of disorder present: for strongly disordered systems such as am-Al2O3 and
vitreous silica [21, 79], models containing hundreds of atoms are sufficient
to reproduce the bulk limit; future work will investigate systems with lower
degree of disorder.

Figure 14:  Negligible effect of anharmonicity on thermal conductivity of am-
Al2O3. Solid, bWTE conductivity of MACE7680, calculated with physical
linewidths. The dashed and dotted lines are bWTE conductivities calculated by
artificially doubling and halving the anharmonic linewidths, respectively. We
see that these artificial rescalings negligibly affect the conductivity,
indicating that structural disorder dominates over anharmonicity in
determining heat transfer in am-Al2O3.

To further validate our assertion that structural disorder is the dominant
limiting factor for heat conduction in am-Al2O3 over the range 50-700 K, we
selected the MACE7680 model—which is large enough to directly compute the bulk
conductivity, without regularization—and calculated the bWTE conductivity in
three ways: (i) with ‘physical’ linewidths computed from first principles
using the AIMD120 model (solid line in Fig. 14); (ii) using artificially
enlarged linewidths, obtained multiplying by a factor 2.0 the physical
linewidths (dashed line); (iii) using artificially reduced linewidths,
obtained multiplying the physical ones by a factor 0.5 (dotted line).
Intuitively, when structural disorder dominates over anharmonicity in
determining heat conduction, we expect the conductivity to be practically
indistinguishable in the three cases above. This behavior is in sharp contrast
with that observed in ordered simple crystals with well separated phonon
bands—in this case anharmonicity determines the conductivity, and in the high-
temperature limit doubling the linewidth directly implies a reduction of the
conductivity by a factor of two (one can verify this analytically by simply
rescaling the linewidths appearing in the thermal conductivity expression for
crystals, e.g. Eq. (49) of Ref. [40]). Fig. 14 shows that artificially
rescaling the linewidths produces practically unnoticeable changes in the
thermal conductivity of am-Al2O3 between 50 and 700 K, showing that in this
system structural disorder dominates over anharmonicity in determining heat
conduction. In other words, keeping in mind that Eq. (3) shows that heat
conduction in amorphous solids is mediated by couplings between vibrational
modes, Fig. 14 shows that in strongly disordered solids such as am-Al2O3
anharmonicity only serves to allow couplings between modes, and as soon as
anharmonicity is large enough (in the thermodynamic limit an infinitesimal
value is sufficient), its exact magnitude is unimportant and does not
influence the value of the conductivity. This behavior is related to the
presence of velocity-operator elements between pairs of vibrational
eigenstates that do not significantly depend on the energy difference between
the eigenstates coupled. These findings are consistent with the analysis
reported in Fig. 7, as well as with the negligible differences observed
between the AF and rWTE conductivities between 50 and 700 K.

##  V Conclusions

The ubiquitous use of am-Al2O3 in electronic devices and the several open
fundamental questions on how its atomistic structural and vibrational
properties determine its macroscopic thermal properties prompted us to study
this material from a first principles level of theory. We generated and
characterized atomistic models of am-Al2O3 from AIMD with densities ranging
from 2.28 g/cm3 to 3.49 g/cm3, describing how the atomic coordination topology
varies with density. We have shown that at least five different atomic
coordination environments coexist in am-Al2O3, and these lead to significant
structural disorder already at the sub-nanometre lengthscale. We have
discussed how the atomic coordination topology affects the vibrational
properties, showing that different coordination environments for oxygen and
aluminium have fingerprints on the coordination-resolved PDOS. We have
described the thermal properties using the recently introduced Voigt-
regularized Wigner formulation (rWTE) [21], accounting comprehensively for the
effects of structural disorder, anharmonicity, and quantum Bose-Einstein
statistics. We have shed light on the microscopic physics underlying thermal
transport in am-Al2O3, discussing the dominant role played by strong
structural disorder (emerging from having at least five coexisting different
atomic coordination topologies) and the negligible role played by
anharmonicity. Specifically, we showed that the harmonic Allen-Feldman
theory—evaluated using the convergence-acceleration protocol discussed in Ref.
[21]—yields predictions in close agreement with the anharmonic rWTE protocol
and with experiments. We have validated these first-principles calculations
using a MACE MLP, generating models of am-Al2O3 containing 3240 and 7680 atoms
at density $\sim 3$ g/cm3, showing that their thermal conductivity is
compatible with that of the 120-atom first-principles models. We discussed how
the increase in the thermal conductivity observed with density derives from an
increase of the vibrational density of states with density, as well as from an
increase of the diffusivity with density. Importantly, we have investigated
the thermal properties also below room temperature ($T\gtrsim 50K$), where the
quantum Bose-Einstein statistics of vibrations yields a specific heat
significantly different from the classical limit, providing information on the
thermal conductivity in a regime inaccessible by molecular-dynamics-based
methods [3, 33], which are governed by classical equipartition [33] and thus
limited to high temperatures. Ultimately, this study further validates the
capability of the recently developed rWTE protocol [21] to describe the
thermal properties of strongly disordered glasses using atomistic models
containing hundreds of atoms, and thus within the reach of first-principles
techniques.

## Acknowledgements

A. F. H. acknowledges the financial support of the Gates Cambridge Trust and
the Winton Programme for the Physics of Sustainability, University of
Cambridge. K.I. acknowledges support from Winton & Cavendish Scholarship at
the Department of Physics, University of Cambridge. W.C.W. and M.C.P.
acknowledge support from the EPSRC (Grant EP/V062654/1). M. S. acknowledges
support from Gonville and Caius College, and from the SNSF project
P500PT_203178. The calculations presented in this work have been performed
using computational resources provided by: (i) the Sulis Tier 2 HPC platform
hosted by the Scientific Computing Research Technology Platform at the
University of Warwick (Sulis is funded by EPSRC Grant EP/T022108/1 and the HPC
Midlands+ consortium); (ii) 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, and GPU resources were obtained through a
University of Cambridge EPSRC Core Equipment Award, EP/X034712/1.

##  Appendix A Computational details

###  A.1 First-principles calculations

The 120-atom am-Al2O3 structures were obtained from the open-source dataset of
Ref. [82]. These models were generated via first-principles molecular dynamics
simulations, using the melt quench procedure described in Ref. [12], with VASP
v5.4 [83]. In order to apply the computational protocol of Ref. [21] to study
the thermal properties, the vibrational frequencies have to be interpolated in
Fourier space. Applying Fourier interpolation to the vibrational frequencies
of disordered atomistic models containing less than 200 atoms yields more
accurate results when a mesh denser than the point $\bm{q}=\bm{0}$ only is
used as starting point. Therefore, to compute the vibrational properties of
am-Al2O3 on a mesh denser than $\bm{q}=\bm{0}$ only and to use the most
accurate density-functional perturbation theory (DFPT) technique [48], the
vibrational properties are computed using Quantum ESPRESSO [84, 85] on a
$2\times 2\times 2$ $\bm{q}$-mesh (at present the DFPT implementation in VASP
is restricted to calculations at $\bm{q}{=}\bm{0}$). Quantum ESPRESSO
calculations were carried out using PBE functional, with pseudopotentials from
the standard solid-state pseudopotential libraries (SSSP) precision library
[86]. To compute the vibrational properties, the cell and atomic positions of
all the am-Al2O3 models were relaxed using a threshold for forces of
2$\times$10-4 Ry/Bohr and of $0.01$ kBar for pressure (vc-relax command).

###  A.2 Generation of the MACE potential

The MACE MLP was trained using the first-principles (PBE functional) dataset
from [3]. We used a two-layer MACE and a per-layer cutoff of 4.5 Å, resulting
in a total receptive field of 9 Å, as well as 128 embedding channels and
$L_{\mathbf{max}}=1$. Both energies and forces were used in training and ten
percent of the full dataset was randomly held out for validation; the
remaining ninety percent comprised the training set. The training proceeds
until the onset of overfitting is observed: we monitor the errors on the
validation set during training and exit when those errors begin to increase.
Following this procedure, the training terminated with a validation error of
5.4 meV/atom for the energies and 123.1 meV/Å for the forces.

##  Appendix B Coordination environment

To determine the coordination topology we calculated the number of atoms in a
sphere of radius equal to the first minimum of the radial distribution
function (Fig. 2), which ranged from e.g., 2.2 Å for the lowest-density
structure, to 2.35 Å for the highest-density structure.

Figure 15: Convergence of the AF conductivity for am-Al2O3 with respect to the
broadening $\eta$ for the Dirac $\delta$, The three upper panels show
convergence plateaus for $5\times 5\times 5$ (lines) and $7\times 7\times 7$
(scatter points) meshes of models with densities 2.28, 2.98 and 3.49
$g/cm^{3}$. Black (red) correspond to evaluating the AF conductivity using the
Gaussian (Lorentzian) representation of the Dirac delta function. The three
lower panels contain calculations at $\bm{q}=\bm{0}$ only for models with
densities 2.28, 2.98 and 3.49 $g/cm^{3}$. Black and red denote Gaussian and
Lorentzian, as in the panels above. The broadening $\eta$ is chosen as the
value approximately determining the beginning of the convergence plateau [21],
and it is set to 6.0 cm-1 for models with density up to 3.30 $g/cm^{3}$ and
8.34 for the model with density $3.49g/cm^{3}$. We see that the $5\times
5\times 5$ $\bm{q}$ mesh is dense enough to achieve computational convergence
with respect to mesh size, since results obtained on a $7\times 7\times 7$
mesh are practically indistinguishable from those obtained using a $5\times
5\times 5$ mesh. Calculations at $\bm{q}=\bm{0}$ only are far from
computational convergence and consequently underestimate the conductivity.

##  Appendix C First-principles thermal conductivity calculations

###  C.1 Convergence of the Allen-Feldman theory

In order to calculate the bulk limit of the thermal conductivity of strongly
disordered solids such as am-Al2O3, we rely on the convergence-acceleration
protocol discussed in Ref. [21] both for the AF and rWTE conductivities. The
capability of such a protocol to accurately extrapolate the bulk limit of the
thermal conductivity of strongly disordered glasses from finite-size models
containing hundreds of atoms is validated in Sec. IV, and in Ref. [21]. The
protocol requires determining the broadening parameters $\eta$ for the Voigt
profile appearing in Eq. (3) as a value determining the beginning of the
convergence plateau shown in Fig. 15 (see Sec. III.1 for details).

All the am-Al2O3 models analyzed display a clear and broad convergence plateau
for the AF conductivity. The three upper panels in Fig. 15 show results
obtained employing a $5\times 5\times 5$ or $7\times 7\times 7$ $\bm{q}$-mesh;
the good agreement between these two calculations indicates that computational
convergence has been achieved.

The three bottom panels of Fig. 15 show that a calculation performed at
$\bm{q}=0$ only for a 120-atom model is far from computational convergence and
underestimates the thermal conductivity. The values of $\eta$ that we
determined from the convergence test discussed in Fig. 15 and that we employed
in our calculations are reported in Table 1.

$\rho$ ($g/cm^{3}$) | 2.28 | 2.98 | 3.17 | 3.30 | 3.49
---|---|---|---|---|---
$\hbar\eta$ (cm-1) | 6.0 | 6.0 | 6.0 | 6.0 | 8.34
Table 1: Broadening parameters $\eta$ used for the Gaussian representation of
the Dirac $\delta$ function appearing in the AF conductivity expression, and
for the Voigt distribution appearing in the rWTE expression.

###  C.2 Radial Distribution Functions

All the radial distribution functions are computed following Ref. [87], using
a Gaussian distribution with standard deviation equal to 0.075 Å.

###  C.3 Anharmonic Linewidths

Figure 16: Effect of temperature on anharmonic linewidths of am-Al2O3 at
various densities and temperatures. The scatter points represent the
linewidths computed at $\bm{q}=\bm{0}$, and the solid lines are coarse grained
functions $\Gamma_{a}[\omega,T]$ used to approximately describe the anharmonic
linewidths as a single-valued functions of frequency, thus to estimate the
effects of anharmonicity at a reduced computational cost[88, 89, 40]. The
purple region denotes the overdamped regime $\Gamma>\omega$ [40, 66]. The gray
dashed lines show the average spacing between vibrational energy levels. We
note that the linewidths of am-Al2O3 are similar to those found other oxide
glasses, e.g., vitreous silica [21].

The third-order interatomic force constants are computed in the 120 atom cells
using ShengBTE [90] up to the 8th nearest neighbor. The linewidths were then
computed using phono3py [57, 91], with a Gaussian smearing of 0.18 THz or ~6
cm-1, and the standard perturbative treatment of anharmonicity, i.e. (i)
vibrational frequencies were considered to be independent from temperature
(i.e. it neglects thermal expansion and the renormalization of frequencies due
to anharmonicity [92, 93, 94, 95, 33]); (ii) the linewidths were computed
considering exclusively the cubic terms in the Taylor expansion of the
interatomic potential [40, 21] and the contribution due to isotopic-mass
disorder [62].

Fig. 16 shows the linewidths for AIMD structures calculated at $\bm{q}=\bm{0}$
at various densities and temperatures. We see that at 50K a significant
proportion of vibrational modes have linewidths smaller than the average
energy-level spacing. As mentioned in section III, these linewidths are
employed within the Voigt distribution, which ensures that heat transfer
between neighboring vibrational eigenstates can always occur, implying that
the effects of anharmonicity are accounted for only when they are not altered
by finite-size effects [21].

The solid lines in Fig. 16 are the functions $\Gamma_{a}[\omega]$ that
approximately describe the anharmonic linewidths as a single-valued functions
of frequency, obtained following the approaches discussed in Ref. [40] (see
also Refs. [88, 89, 47] for similar approximated treatments of anharmonicity).
The approximated function $\Gamma_{a}[\omega]$ is employed to compute the
anharmonic linewidths as a function of frequency when the Fourier
interpolation is used to extrapolate the bulk limit of the thermal
conductivity 3, following the protocol discussed in Ref. [40].

###  C.4 rWTE conductivity calculation

The quantities needed to evaluate the rWTE conductivity (3) were computed as
follows: (i) the $\eta$ parameter was computed as discussed in Sec. C.1; (ii)
the software phono3py [57, 91] was used to evaluate frequencies and velocity
operators on a $\bm{q}$ mesh; (iii) the linewidths were evaluated from the
frequencies determined at the previous point using the function
$\Gamma_{a}(\omega)$ discussed in Sec. C.3.

We checked that increasing the size of the $\bm{q}$ mesh from $5\times 5\times
5$ to $7\times 7\times 7$ produced practically indistinguishable results in
the RMDS, VDOS, and thermal conductivities. Therefore, all the results
discussed in the main text are evaluated on a $5\times 5\times 5$ $\bm{q}$
mesh.

We limited our calculations to temperatures higher than 50 K, since at
temperatures lower than 50 K the thermal properties of am-Al2O3 are dominated
by low-frequency vibrational modes that are likely to feature glassy anomalies
[96, 97, 98]; accurately sampling these low-frequencies vibrational modes
requires using atomistic models containing thousands of atoms, and it is
therefore beyond the scope of the present work.

$D(\omega)$ data for the plot in Fig. 8 was calculated using Eq. (6) with
$\delta$-function smeared to a Gaussian with variance
$\frac{\pi}{2}\eta_{0}^{2}$, with $\hbar\eta_{0}=8.34$ cm-1.

##  Appendix D Thermal conductivity calculations using MACE

Figure 17: Convergence of the AF conductivity for MACE models of am-Al2O3 with
respect to the broadening $\eta$ for the Dirac $\delta$, The three upper
panels show convergence for 5x5x5 and 3x3x3 meshes of MACE models containing
120, 3240 and 7680 atoms. Black (red) correspond to evaluating the AF
conductivity using the Gaussian (Lorentzian) representation of the Dirac
delta. The three lower panels contain calculations at $\bm{q}{=}0$ for MACE
models with 120, 3240 and 7680 atoms. Black and red denote Gaussian and
Lorentzian, as in the panels above. The broadening $\eta$ is chosen as the
value approximately determining the beginning of the convergence plateau [21],
and it is set to 6.0 cm-1 for MACE120, 2.0 cm-1 for MACE3240, and 1.0 cm-1 for
MACE7680; the values of $\eta$ are shown as dotted vertical black lines. We
see that a model size of 3240 atoms is sufficient to obtain equivalent results
of AF thermal conductivity using the $\bm{q}$-mesh or $\bm{q}{=}0$ point only,
and using larger models extends the plateau to lower values of $\eta$, as in
vitreous silica [21].

For the MACE models, we present their convergence tests in Fig. 17. All three
structures display a broad convergence plateau for AF conductivity on a
$\bm{q}$-mesh that is increased in length as one increases the system size. We
further note that 3240-atom and 7680-atom structures are both big enough to
exhibit a convergence plateau at $\bm{q}{=}0$ point only that is compatible
with the plateau at 3x3x3 $\bm{q}$-mesh indicating achievement of
computational convergence.

Given that Fig. 14 shows that anharmonicity has negligible effects on the
conductivity of am-Al2O3, in the calculation of bWTE and rWTE conductivity for
120, 3240 and 7680-atom MACE models we used the coarse-grained functions for
the anharmonic linewidths derived from the 2.98 g/cm3 AIMD model (central
panel in Fig. 16).

##  Appendix E Velocity operator for MACE and AIMD models at densities close
to 3.0 g/cm3

To compare predictions of AIMD and MACE for the velocity operator, we plot in
Fig. 18 the velocity operator represented as a function of frequency
difference and average
($\langle|\nu^{\mathrm{avg}}_{\omega_{a}\omega_{d}}|^{2}\rangle$, Eq. 8). In
the upper panel, we can see comparison between AIMD and MACE structures with
120 atoms each. We note the agreement between different methods is
satisfactory, especially between 200 and 700 cm-1, where vibrational DOS is at
its maximum. In the middle panel of Fig. 18 we compare $\bm{q}{=}\bm{0}$ and
5x5x5 $\bm{q}$-mesh calculations of velocity operator elements for 120-atom
MACE model. We see that $\bm{q}{=}\bm{0}$ calculation underestimates velocity
operator elements for frequencies below 150 cm-1, ultimately leading to
underestimate of thermal conductivity as discussed in Fig. 13. In the lower
panel of Fig. 18 we compare the velocity operator elements for a small
120-atom model averaged over a $\bm{q}$-mesh, with those at $\bm{q}{=}\bm{0}$
for a large 7680-atom model: these are overall similar, explaining the
compatible predictions for rWTE conductivities shown in Fig. 13, and
validating the idea of obtaining velocity operator elements averaging over
different boundary conditions in small models of disordered solids.

Figure 18:  Velocity operator of am-Al2O3 for 120-atom AIMD model, and 120-
and 7680-atom MACE models. The square modulus of the velocity operator
$\langle|\nu^{\mathrm{avg}}_{\omega_{a}\omega_{d}}|^{2}\rangle$ is represented
as a function of energy differences
$\hbar\omega_{d}=\hbar(\omega(\mathbf{q}_{s})-\omega(\mathbf{q}_{s^{\prime}}))$
and averages
$\hbar\omega_{a}=\hbar\frac{\omega(\mathbf{q}_{s})+\omega(\mathbf{q}_{s^{\prime}})}{2}$,
following Eq. (8). The upper panel shows a reasonable agreement between
velocity operator elements of the 120-atom models generated with MACE (solid)
and AIMD (dashed), both evaluated on a 5x5x5 $\bm{q}$-mesh. The middle panel
shows that for MACE120, the velocity operator at $\bm{q}{=}0$ (dashed) is
overall smaller than the velocity operator computed over a 5x5x5 $\bm{q}$-mesh
(solid). The bottom panel shows that the velocity operator for MACE120
computed over a 5x5x5 $\bm{q}$-mesh (solid) is in reasonable agreement with
the velocity operator of MACE7680 computed at $\bm{q}{=}0$ (dashed). The delta
functions used for calculation of
$\langle|\nu^{\mathrm{avg}}_{\omega_{a}\omega_{d}}|^{2}\rangle$ were replaced
with Gaussian with variance related to the smearing parameter from the
convergence plateau: $\sigma^{2}=\eta^{2}\pi/2$.

##  Appendix F Allen-Feldman conductivity of MACE models

To supplement the comparison of rWTE conductivities in Fig. 13, we present AF
conductivities for AIMD and MACE models in Fig. 19. We first note that for am-
Al2O3 at densities close to 3.0 g/cm3, influence of anharmonicity on the
magnitude of conductivity is very weak, which is exemplified by very good
agreement between WTE $\bm{q}{=}0$ and AF $\bm{q}{=}0$ predictions for 3240
and 7680-atom MACE models between 50 and 700 K. Our second conclusion is that
predictions of AF conductivity using a $\bm{q}$-mesh in small models give the
same trend and very similar magnitude as predictions of AF conductivity using
$\bm{q}{=}0$ in large models for very disordered solids, which am-Al2O3 is an
example of. The largest differences of value of AF conductivity at 700K are
between MACE7680 at $\bm{q}{=}0$ and on 3x3x3 $\bm{q}$-mesh , and is equal to
approximately $0.04$ W/mK ($2$ %).

Figure 19:  AF thermal conductivity for models of am-Al2O3 with densities
close to 3.0 g/cm3. Solid (dashed) lines are AF (WTE) calculations done at
$\bm{q}{=}0$ for 3240-atom (orange) and 7680-atom (purple) MACE models. The
scatter points are AF harmonic conductivities calculated on a mesh for
120-atom AIMD model (green, empty squares, 5x5x5 mesh), 120-atom MACE model
(blue, filled diamonds, 5x5x5 mesh), 3240-atom MACE model (orange, inverted,
filled triangles, 3x3x3 mesh) and 7680-atom MACE model (empty, purple circles,
3x3x3 mesh).

## References

  * Levin and Brandon [1998] I. Levin and D. Brandon, Metastable alumina polymorphs: Crystal structures and transition sequences, Journal of the American Ceramic Society 81, 1995 (1998).
  * Tavakoli  _et al._ [2013] A. H. Tavakoli, P. S. Maram, S. J. Widgeon, J. Rufner, K. Van Benthem, S. Ushakov, S. Sen, and A. Navrotsky, Amorphous alumina nanoparticles: structure, surface energy, and thermodynamic phase stability, The Journal of Physical Chemistry C 117, 17123 (2013).
  * Li  _et al._ [2020] W. Li, Y. Ando, and S. Watanabe, Effects of density and composition on the properties of amorphous alumina: A high-dimensional neural network potential study, The Journal of Chemical Physics 153, 164119 (2020).
  * Shafiq  _et al._ [2022] I. Shafiq, S. Shafique, P. Akhter, W. Yang, and M. Hussain, Recent developments in alumina supported hydrodesulfurization catalysts for the production of sulfur-free refinery products: A technical review, Catalysis Reviews 64, 1 (2022).
  * Oh  _et al._ [2019] J. Oh, H. B. Bathula, J. H. Park, and Y.-W. Suh, A sustainable mesoporous palladium-alumina catalyst for efficient hydrogen release from n-heterocyclic liquid organic hydrogen carriers, Communications chemistry 2, 68 (2019).
  * Jin  _et al._ [2019] Y. Jin, H. Yu, Y. Gao, X. He, T. A. White, and X. Liang, Li4Ti5O12 coated with ultrathin aluminum-doped zinc oxide films as an anode material for lithium-ion batteries, Journal of Power Sources 436, 226859 (2019).
  * Qu  _et al._ [2019] S. Qu, W. Jia, Y. Wang, C. Li, Z. Yao, K. Li, Y. Liu, W. Zou, F. Zhou, Z. Wang, _et al._ , Air-stable lithium metal anode with sputtered aluminum coating layer for improved performance, Electrochimica Acta 317, 120 (2019).
  * Randau  _et al._ [2021] S. Randau, F. Walther, A. Neumann, Y. Schneider, R. S. Negi, B. Mogwitz, J. Sann, K. Becker-Steinberger, T. Danner, S. Hein, _et al._ , On the additive microstructure in composite cathodes and alumina-coated carbon microwires for improved all-solid-state batteries, Chemistry of Materials 33, 1380 (2021).
  * Paterson  _et al._ [2020] J. Paterson, D. Singhal, D. Tainoff, J. Richard, and O. Bourgeois, Thermal conductivity and thermal boundary resistance of amorphous Al${}_{\textrm{2}}$O${}_{\textrm{3}}$ thin films on germanium and sapphire, Journal of Applied Physics 127, 245105 (2020).
  * Mavrič  _et al._ [2019] A. Mavrič, M. Valant, C. Cui, and Z. M. Wang, Advanced applications of amorphous alumina: From nano to bulk, Journal of Non-Crystalline Solids 521, 119493 (2019).
  * Scott  _et al._ [2018] E. A. Scott, J. T. Gaskins, S. W. King, and P. E. Hopkins, Thermal conductivity and thermal boundary resistance of atomic layer deposited high-k dielectric aluminum oxide, hafnium oxide, and titanium oxide thin films on silicon, APL Materials 6, 058302 (2018).
  * F. Harper  _et al._ [2023] A. F. Harper, S. P. Emge, P. C. M. M. Magusin, C. P. Grey, and A. J. Morris, Modelling amorphous materials via a joint solid-state NMR and X-ray absorption spectroscopy and DFT approach: application to alumina, Chemical Science 14, 1155 (2023).
  * Lee  _et al._ [2010] S. K. Lee, S. Y. Park, Y. S. Yi, and J. Moon, Structure and disorder in amorphous alumina thin films: Insights from high-resolution solid-state NMR, The Journal of Physical Chemistry C 114, 13890 (2010).
  * Sarou-Kanian  _et al._ [2013] V. Sarou-Kanian, A. N. Gleizes, P. Florian, D. Samélor, D. Massiot, and C. Vahlas, Temperature-dependent 4-, 5- and 6-fold coordination of aluminum in mocvd-grown amorphous alumina films: A very high field 27Al-NMR study, The Journal of Physical Chemistry C 117, 21965 (2013).
  * Hashimoto  _et al._ [2022] H. Hashimoto, Y. Onodera, S. Tahara, S. Kohara, K. Yazawa, H. Segawa, M. Murakami, and K. Ohara, Structure of alumina glass, Scientific Reports 12, 516 (2022), number: 1 Publisher: Nature Publishing Group.
  * Shi  _et al._ [2019] C. Shi, O. L. Alderman, D. Berman, J. Du, J. Neuefeind, A. Tamalonis, J. R. Weber, J. You, and C. J. Benmore, The structure of amorphous and deeply supercooled liquid alumina, Frontiers in Materials 6, 38 (2019).
  * Dicks  _et al._ [2019] O. A. Dicks, J. Cottom, A. L. Shluger, and V. V. Afanas’ev, The origin of negative charging in amorphous Al2O3 films: the role of native defects, Nanotechnology 30, 205201 (2019).
  * Leung [2021] K. Leung, First principles, explicit interface studies of oxygen vacancy and chloride in alumina films for corrosion applications, Journal of The Electrochemical Society 168, 031511 (2021).
  * Cahill and Pohl [1987] D. G. Cahill and R. O. Pohl, Thermal conductivity of amorphous solids above the plateau, Physical Review B 35, 4067 (1987), publisher: American Physical Society.
  * Allen and Feldman [1989] P. B. Allen and J. L. Feldman, Thermal Conductivity of Glasses: Theory and Application to Amorphous Si, Phys. Rev. Lett. 62, 645 (1989).
  * Simoncelli  _et al._ [2023] M. Simoncelli, F. Mauri, and N. Marzari, Thermal conductivity of glasses: first-principles theory and applications, npj Computational Materials 9, 106 (2023).
  * Ercole  _et al._ [2016] L. Ercole, A. Marcolongo, P. Umari, and S. Baroni, Gauge Invariance of Thermal Transport Coefficients, Journal of Low Temperature Physics 185, 79 (2016).
  * Knoop  _et al._ [2023a] F. Knoop, M. Scheffler, and C. Carbogno, Ab initio Green-Kubo simulations of heat transport in solids: Method and implementation, Physical Review B 107, 224304 (2023a), publisher: American Physical Society.
  * Knoop  _et al._ [2023b] F. Knoop, T. A. Purcell, M. Scheffler, and C. Carbogno, Anharmonicity in Thermal Insulators: An Analysis from First Principles, Physical Review Letters 130, 236301 (2023b), publisher: American Physical Society.
  * Marcolongo  _et al._ [2016] A. Marcolongo, P. Umari, and S. Baroni, Microscopic theory and quantum simulation of atomic heat transport, Nat. Phys. 12, 80 (2016).
  * Eriksson  _et al._ [2023] F. Eriksson, E. Fransson, C. Linderälv, Z. Fan, and P. Erhart, Tuning the Through-Plane Lattice Thermal Conductivity in van der Waals Structures through Rotational (Dis)ordering, ACS Nano , acsnano.3c09717 (2023).
  * Lv and Henry [2016] W. Lv and A. Henry, Non-negligible contributions to thermal conductivity from localized modes in amorphous silicon dioxide, Sci. Rep. 6, 35720 (2016).
  * Donadio and Galli [2009] D. Donadio and G. Galli, Atomistic Simulations of Heat Transport in Silicon Nanowires, Phys. Rev. Lett. 102, 195901 (2009).
  * R. Hahn  _et al._ [2022] K. R. Hahn, C. Melis, and L. Colombo, Thermal conduction and rectification phenomena in nanoporous silicon membranes, Physical Chemistry Chemical Physics 24, 13625 (2022), publisher: Royal Society of Chemistry.
  * Jund and Jullien [1999] P. Jund and R. Jullien, Molecular-dynamics calculation of the thermal conductivity of vitreous silica, Phys. Rev. B 59, 13707 (1999).
  * Shenogin  _et al._ [2009] S. Shenogin, A. Bodapati, P. Keblinski, and A. J. McGaughey, Predicting the thermal conductivity of inorganic and polymeric glasses: The role of anharmonicity, J. Appl. Phys. 105, 034906 (2009).
  * Felix and Pereira [2018] I. M. Felix and L. F. C. Pereira, Thermal Conductivity of Graphene-hBN Superlattice Ribbons, Scientific Reports 8, 2737 (2018).
  * Puligheddu  _et al._ [2019] M. Puligheddu, Y. Xia, M. Chan, and G. Galli, Computational prediction of lattice thermal conductivity: A comparison of molecular dynamics and Boltzmann transport approaches, Phys. Rev. Materials 3, 085401 (2019).
  * Batatia  _et al._ [2022] I. Batatia, D. P. Kovacs, G. Simm, C. Ortner, and G. Csányi, Mace: Higher order equivariant message passing neural networks for fast and accurate force fields, Advances in Neural Information Processing Systems 35, 11423 (2022).
  * Kovacs  _et al._ [2023] D. P. Kovacs, I. Batatia, E. S. Arany, and G. Csanyi, Evaluation of the mace force field architecture: from medicinal chemistry to materials science, arXiv preprint arXiv:2305.14247 (2023).
  * Drautz [2019] R. Drautz, Atomic cluster expansion for accurate and transferable interatomic potentials, Physical Review B 99, 014104 (2019), publisher: American Physical Society.
  * Dusson  _et al._ [2022] G. Dusson, M. Bachmayr, G. Csányi, R. Drautz, S. Etter, C. van der Oord, and C. Ortner, Atomic cluster expansion: Completeness, efficiency and stability, Journal of Computational Physics 454, 110946 (2022).
  * Witt  _et al._ [2023] W. C. Witt, C. van der Oord, E. Gelžinytė, T. Järvinen, A. Ross, J. P. Darby, C. H. Ho, W. J. Baldwin, M. Sachs, J. Kermode, N. Bernstein, G. Csányi, and C. Ortner, ACEpotentials.jl: A Julia implementation of the atomic cluster expansion, The Journal of Chemical Physics 159, 164101 (2023).
  * Simoncelli, M. and Marzari, N. and Mauri, F. [2019] Simoncelli, M. and Marzari, N. and Mauri, F., Unified theory of thermal transport in crystals and glasses, Nat. Phys. 15, 809 (2019).
  * Simoncelli  _et al._ [2022] M. Simoncelli, N. Marzari, and F. Mauri, Wigner formulation of thermal transport in solids, Phys. Rev. X 12, 041011 (2022).
  * Bhatia  _et al._ [1989] C. Bhatia, G. Guthmiller, and A. Spool, Alumina films by sputter deposition with ar/o2: Preparation and characterization, Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 7, 1298 (1989).
  * Koski  _et al._ [1999] K. Koski, J. HoÈlsaÈ, and P. Juliet, Properties of aluminium oxide thin ®lms deposited by reactive magnetron sputtering, Thin Solid Films (1999).
  * Paz  _et al._ [2014] A. P. Paz, I. V. Lebedeva, I. V. Tokatly, and A. Rubio, Identification of structural motifs as tunneling two-level systems in amorphous alumina at low temperatures, Physical Review B 90, 224202 (2014).
  * Elliott [1984] S. R. Elliott,  _Physics of amorphous materials_ (Longman Group Ltd., 1984).
  * Young  _et al._ [2020] M. J. Young, N. M. Bedford, A. Yanguas-Gil, S. Letourneau, M. Coile, D. J. Mandia, B. Aoun, A. S. Cavanagh, S. M. George, and J. W. Elam, Probing the Atomic-Scale Structure of Amorphous Aluminum Oxide Grown by Atomic Layer Deposition, ACS Applied Materials & Interfaces 12, 22804 (2020).
  * Lizárraga  _et al._ [2011] R. Lizárraga, E. Holmström, S. C. Parker, and C. Arrouvel, Structural characterization of amorphous alumina and its polymorphs from first-principles xps and nmr calculations, Phys. Rev. B 83, 094201 (2011).
  * Fiorentino  _et al._ [2023] A. Fiorentino, P. Pegolo, and S. Baroni, Hydrodynamic finite-size scaling of the thermal conductivity in glasses, npj Computational Materials 9, 1 (2023), number: 1 Publisher: Nature Publishing Group.
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Reviews of Modern Physics 73, 515 (2001).
  * Wallace [1998] D. Wallace,  _Thermodynamics of Crystals_, Dover books on physics (Dover Publications, 1998).
  * Note [1] N is equal to the number of $\bm{q}$ points used to sample the Brillouin Zone; in the case of an ideal glass (astronomically large disordered simulation cell) the Brillouin Zone reduces to the point $\bm{q}=\bm{0}$ only and $N=1$.
  * Łodziana and Parliński [2003] Z. Łodziana and K. Parliński, Dynamical stability of the $\alpha$ and $\theta$ phases of alumina, Physical Review B 67, 174106 (2003).
  * Simoncelli  _et al._ [2020] M. Simoncelli, N. Marzari, and A. Cepellotti, Generalization of fourier’s law into viscous heat equations, Phys. Rev. X 10, 011019 (2020).
  * Di Lucente  _et al._ [2023] E. Di Lucente, M. Simoncelli, and N. Marzari, Crossover from boltzmann to wigner thermal transport in thermoelectric skutterudites, arXiv.2303.07019 (2023).
  * Liu  _et al._ [2023] Y. Liu, H. Liang, L. Yang, G. Yang, H. Yang, S. Song, Z. Mei, G. Csányi, and B. Cao, Unraveling Thermal Transport Correlated with Atomistic Structures in Amorphous Gallium Oxide via Machine Learning Combined with Experiments, Advanced Materials n/a, 2210873 (2023).
  * Paulatto  _et al._ [2013] L. Paulatto, F. Mauri, and M. Lazzeri, Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene, Phys. Rev. B 87, 214303 (2013).
  * Fugallo  _et al._ [2013] G. Fugallo, M. Lazzeri, L. Paulatto, and F. Mauri, Ab initio variational approach for evaluating lattice thermal conductivity, Phys. Rev. B 88, 045430 (2013).
  * Togo  _et al._ [2015] A. Togo, L. Chaput, and I. Tanaka, Distributions of phonon lifetimes in brillouin zones, Phys. Rev. B 91, 094306 (2015).
  * Tadano  _et al._ [2014] T. Tadano, Y. Gohda, and S. Tsuneyuki, Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations, J. Phys. Condens. Matter 26, 225402 (2014).
  * Cepellotti  _et al._ [2022] A. Cepellotti, J. Coulter, A. Johansson, N. S. Fedorova, and B. Kozinsky, Phoebe: a high-performance framework for solving phonon and electron Boltzmann transport equations, Journal of Physics: Materials 5, 035003 (2022), publisher: IOP Publishing.
  * Carrete  _et al._ [2017] J. Carrete, B. Vermeersch, A. Katre, A. van Roekeghem, T. Wang, G. K. Madsen, and N. Mingo, almaBTE : A solver of the space–time dependent Boltzmann transport equation for phonons in structured materials, Computer Physics Communications 220, 351 (2017).
  * Barbalinardo  _et al._ [2020] G. Barbalinardo, Z. Chen, N. W. Lundgren, and D. Donadio, Efficient anharmonic lattice dynamics calculations of thermal transport in crystalline and disordered solids, Journal of Applied Physics 128, 135104 (2020).
  * Tamura [1983] S.-i. Tamura, Isotope scattering of dispersive phonons in Ge, Physical Review B 27, 858 (1983).
  * Note [2] A necessary condition for this to happen is to have vibrations that are not localized in the Anderson sense [99], i.e. to have non-zero velocity operator elements in Eq. (3).
  * Allen and Feldman [1993] P. B. Allen and J. L. Feldman, Thermal conductivity of disordered harmonic solids, Phys. Rev. B 48, 12581 (1993).
  * Feldman  _et al._ [1993] J. L. Feldman, M. D. Kluge, P. B. Allen, and F. Wooten, Thermal conductivity and localization in glasses: Numerical study of a model of amorphous silicon, Physical Review B 48, 12589 (1993).
  * Caldarelli  _et al._ [2022] G. Caldarelli, M. Simoncelli, N. Marzari, F. Mauri, and L. Benfatto, Many-body Green’s function approach to lattice thermal transport, Physical Review B 106, 024312 (2022).
  * Fiorentino and Baroni [2023] A. Fiorentino and S. Baroni, From green-kubo to the full boltzmann kinetic approach to heat transport in crystals and glasses, Phys. Rev. B 107, 054311 (2023).
  * Isaeva  _et al._ [2019] L. Isaeva, G. Barbalinardo, D. Donadio, and S. Baroni, Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach, Nat. Commun. 10, 3853 (2019).
  * Lundgren  _et al._ [2021] N. W. Lundgren, G. Barbalinardo, and D. Donadio, Mode localization and suppressed heat transport in amorphous alloys, Physical Review B 103, 024204 (2021).
  * Note [3] More precisely , Refs. [68, 69] evaluated the bulk limit of the conductivity relying on empirical interatomic potentials and atomistic models containing thousands of atoms, i.e. having a size large enough to achieve computational convergence by evaluating Eq. (3) at $\bm{q}=\bm{0}$ only and without relying on the Voigt regularization. See Ref. [21] for details on the conditions under which evaluating Eq. (3) yields equivalent results with or without relying on the Fourier interpolation and Voigt regularization.
  * Lee  _et al._ [1995] S.-M. Lee, D. G. Cahill, and T. H. Allen, Thermal conductivity of sputtered oxide films, Phys. Rev. B 52, 253 (1995).
  * Monachon and Weber [2015] C. Monachon and L. Weber, Influence of a nanometric Al2O3 interlayer on the thermal conductance of an al/(si, diamond) interface, Advanced Engineering Materials 17, 68 (2015).
  * Gorham  _et al._ [2014] C. S. Gorham, J. T. Gaskins, G. N. Parsons, M. D. Losego, and P. E. Hopkins, Density dependence of the room temperature thermal conductivity of atomic layer deposition-grown amorphous alumina (Al2O3), Applied Physics Letters 104, 253107 (2014).
  * Fournier  _et al._ [2020] D. Fournier, M. Marangolo, and C. Fretigny, Measurement of thermal properties of bulk materials and thin films by modulated thermoreflectance (MTR), Journal of Applied Physics 128, 241101 (2020).
  * Koh and Cahill [2007] Y. K. Koh and D. G. Cahill, Frequency dependence of the thermal conductivity of semiconductor alloys, Physical Review B 76, 075207 (2007).
  * Zhou  _et al._ [2020] W.-X. Zhou, Y. Cheng, K.-Q. Chen, G. Xie, T. Wang, and G. Zhang, Thermal conductivity of amorphous materials, Advanced Functional Materials 30, 1903829 (2020).
  * Note [4] We note, in passing, that Refs. [72, 78, 73] did not observe a significant dependence of the thermal conductivity from the substrate on which the am-Al2O3 sample was grown.
  * Lee  _et al._ [2017] S.-M. Lee, W. Choi, J. Kim, T. Kim, J. Lee, S. Y. Im, J. Y. Kwon, S. Seo, M. Shin, and S. E. Moon, Thermal conductivity and thermal boundary resistances of ald Al2O3 films on Si and sapphire, International Journal of Thermophysics 38, 176 (2017).
  * Ndour  _et al._ [2023] M. Ndour, P. Jund, and L. Chaput, Practical approach to thermal conductivity calculations of small SiO2 samples, Journal of Non-Crystalline Solids 621, 122618 (2023).
  * Note [5] The largest difference between the converged conductivities is of about 3.3$\%$ ($\approx$ 0.057 W/mK) observed at 700 K between the rWTE computed over a 5x5x5 $\bm{q}$-mesh for MACE120, and rWTE computed over a 3x3x3 $\bm{q}$-mesh for MACE7680.
  * Simkin and Mahan [2000] M. V. Simkin and G. D. Mahan, Minimum Thermal Conductivity of Superlattices, Physical Review Letters 84, 927 (2000).
  * Harper  _et al._ [2023] A. Harper, S. P. Emge, P. Magusin, C. P. Grey, and A. J. Morris, Research data supporting" modelling amorphous materials via a joint solid-state NMR and X-ray absorption spectroscopy and DFT approach: application to alumina",  (2023).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Physical review B 47, 558 (1993).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Prandini  _et al._ [2018] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations, npj Computational Materials 4, 72 (2018).
  * Le Roux and Jund [2010] S. Le Roux and P. Jund, Ring statistics analysis of topological networks: New approach and application to amorphous GeS2 and SiO2 systems, Computational Materials Science 49, 70 (2010).
  * Thébaud  _et al._ [2022] S. Thébaud, T. Berlijn, and L. Lindsay, Perturbation theory and thermal transport in mass-disordered alloys: Insights from green’s function methods, Phys. Rev. B 105, 134202 (2022).
  * Garg  _et al._ [2011] J. Garg, N. Bonini, B. Kozinsky, and N. Marzari, Role of Disorder and Anharmonicity in the Thermal Conductivity of Silicon-Germanium Alloys: A First-Principles Study, Phys. Rev. Lett. 106, 045901 (2011).
  * Li  _et al._ [2014] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, ShengBTE: A solver of the Boltzmann transport equation for phonons, Comput. Phys. Commun. 185, 1747 (2014).
  * Togo [2023] A. Togo, First-principles Phonon Calculations with Phonopy and Phono3py, Journal of the Physical Society of Japan 92, 012001 (2023).
  * Monacelli  _et al._ [2021] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects, Journal of Physics: Condensed Matter 33, 363001 (2021).
  * Tadano and Saidi [2022] T. Tadano and W. A. Saidi, First-Principles Phonon Quasiparticle Theory Applied to a Strongly Anharmonic Halide Perovskite, Physical Review Letters 129, 185901 (2022).
  * Jain [2020] A. Jain, Multichannel thermal transport in crystalline Tl 3 VSe 4, Physical Review B 102, 201201 (2020).
  * Feng  _et al._ [2017] T. Feng, L. Lindsay, and X. Ruan, Four-phonon scattering significantly reduces intrinsic thermal conductivity of solids, Physical Review B 96, 161201 (2017).
  * Schirmacher [2006] W. Schirmacher, Thermal conductivity of glassy materials and the boson peak, EPL (Europhysics Letters) 73, 892 (2006).
  * Lubchenko and Wolynes [2003] V. Lubchenko and P. G. Wolynes, The origin of the boson peak and thermal conductivity plateau in low-temperature glasses, Proc. Natl. Acad. Sci. U.S.A. 100, 1515 (2003).
  * Wang  _et al._ [2019] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Low-frequency vibrational modes of stable glasses, Nature Communications 10, 26 (2019).
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).