Vibrational and thermal properties of amorphous alumina from first principles

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 (Al 2 O 3 ) 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-Al 2 O 3 ) 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-Al 2 O 3 have become subject of extensive research, since the atomiclayer 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-Al 2 O 3 affects its macroscopic properties.Solid-state 27 Al nuclear magnetic resonance has identified three main aluminum coordination environments in am-Al 2 O 3 [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-Al 2 O 3 at technologically relevant temperatures (i.e. at T > 30 K, in the so called above-theplateau 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 nonequilibrium (NEMD [29][30][31][32]) approaches) to compute the room-temperature thermal conductivity of atomistic models of am-AlO x ; they studied models with densities ranging from 2.6 to 3.3 g/cm 3 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-Al 2 O 3 across a wide range of densities (2.28≤ρ≤3.49g/cm 3 ).(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-Al 2 O 3 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-Al 2 O 3 , 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-Al 2 O 3 , demonstrating that disorder practically determines thermal transport at all the densities and temperatures studied (50≤T ≤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-Al 2 O 3 containing 120, 3240, and 7680 atoms.We discuss how the convergence-acceleration protocol of Ref. [21] affects the thermal-conductivity predic-arXiv:2303.08637v2[cond-mat.mtrl-sci]24 Dec 2023 tions in models having different sizes.

II. STRUCTURAL AND VIBRATIONAL PROPERTIES
The properties of am-Al 2 O 3 strongly depend on the synthesis process-samples with densities ranging from 2.1 g/cm 3 to 3.6 g/cm 3 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/cm 3 , 3.30 g/cm 3 and 3.49 g/cm 3 from Ref. [12], and additionally generated two 120-atom am-Al 2 O 3 models at 2.28 g/cm 3 and 2.98 g/cm 3 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.

A. 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-Al 2 O 3 models, showing in Fig. 1 the probability distribution of finding x-fold-coordinated oxygen (O x ) or y-fold-coordinated aluminum (Al y ) at each density.Six different coordination environments are present in our models-O x with x ∈ 2, 3, 4, and Al y with y ∈ 4, 5, 6-and at least five out of these six coordination topologies coexist at every model density analyzed.
As shown in the top panel of Fig. 1, at increasing densities the concentration of twofold-coordinated oxygens (O 2 ) generally decreases in favor of fourfold-coordinated (O 4 ); changes in density have weaker effects on the concentration of threefold-coordinated oxygens (O 3 ).The bottom panel shows that increasing density yields a reduction in the concentration of Al 4 environments, which is compensated by an increase of Al 6 .In contrast, the Al 5 environment displays weaker variations with density; we note that Al 5 are characteristic of the amorphous phase of Al 2 O 3 and absent in the crystalline phases [43].Finally, we highlight how the 2.28 g/cm 3 model contains a high proportion of O 2 bridging oxygens and a majority of Al 4 environments, with no Al 6 .This is consistent with the trend that at lower densities Al becomes primarily tetrahedrally coordinated in Al 4 environments [10].
From the coexistence of at least five different coordination environments in all the am-Al 2 O 3 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.2.28 g/cm 3   2.98 g/cm 3   3.17 g/cm 3   3.30 g/cm 3   3.49 g/cm 3 FIG.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).
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 ∼ 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.

B. Vibrational properties
The strong variations in the coordination topology observed in the different structures of am-Al 2 O 3 (Fig. 1) are intuitively expected to affect the atomic vibrational excitations.Thus, in this section we investigate quantitatively how the vibrational properties of am-Al 2 O 3 change with the coordination topology.We used densityfunctional perturbation theory (DFPT) [48] to compute the energy ℏω qs of each vibrational mode qs, and the corresponding displacement pattern E bα qs , which describes how atom b moves in direction α when the vibration qs is excited.We label vibrational modes using both the wavevector q and the mode index s for the sake of generality, keeping in mind that 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 q = 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], (1) where N qs =[exp (ℏω qs /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 factor [50].
To see if a relationship between RMSD and coordination topology exists, we compute the average RMSD(b, T ) over atoms with equal coordination.Fig. 3 shows that the average RMSD of O x atoms decreases as coordination increases in all the atomistic models analyzed.The average RMSD of Al y is anticorrelated with coordination in models with ρ≤3.17 g/cm 3 , and displays deviations from this anticorrelated trend in models with ρ = 3.30 g/cm 3 (where all Al atoms have a very similar average RMSD, regardless of their coordination) and with ρ = 3.49 g/cm 3 (where Al 6 has average RMSD slightly larger than Al 5 ).Interestingly, as density increases the average RMSD of O 4 approaches the RMSD of Al 4 -considering the significant differences between the mass of oxygen and aluminum (m Al /m O ≈ 1.7), we note that this behavior departs from the trend RMSD ∝ √ 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(ω)=(VN c ) −1 q,s δ(ω−ω qs ) (here, 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: where δ b,tx 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 tx (ω) allows us to resolve how the coordination topology affects the VDOS.Results for g tx (ω) are shown in Fig. 4. We highlight how the PDOS for different coordination environments have similar shapes at different densities.In particular,

III. THERMAL PROPERTIES A. 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 [40,52,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]: where ℏω qs denotes the harmonic energy of the vibration qs, and ℏΓ qs its linewidth (broadening of the harmonic energy level due to anharmonicity [55][56][57][58][59][60][61] and isotopic-mass disorder [62]); ∥v (q) s,s ′ ∥ 2 = 3 α=1 v α (q) s,s ′ v α (q) s ′ ,s is the square modulus of the velocity operator [40] between eigenstates s and s ′ at fixed q (α 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 κ αβ [21]); is the quantum specific heat of the mode qs; V, N c , and 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 η in the low-temperature (harmonic) limit where anharmonicity phases out (Γ qs +Γ qs ′ → 0 ∀ q, s).Setting η to a value slightly larger than the average energy-level spacing ℏ∆ω avg = ℏωmax 3•Nat (ℏω max is the maximum vibrational energy, 3N 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 [63] 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 η, the Voigt profile becomes a Lorentzian with FWHM Γ qs +Γ qs ′ , 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 ℏ∆ω avg > Γ qs , 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 at .Thus, the Brillouin Zone (BZ) corresponding to the (large) finitesize model does not reduce to the aforementioned "ideal glass" limit (q=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 q=0 only correspond to considering periodic boundary conditions at the boundaries of the simulation box; calculations at q=(1, 0, 0)π/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 q=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-Al 2 O 3 , 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 [70].

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-Al 2 O 3 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).
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 η 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.49g/cm 3 ) model.These differences become smaller as the density decreases, and are practically invisible in the lowest-density (2.28 g/cm 3 ) model.The good agreement between AF and rWTE shows that in am-Al 2 O 3 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 temperatureconductivity curve; such a saturating trend is inherited from that of the specific heat (more on this later in Sec.III B 2).Our calculations shed light on the thermal properties of am-Al 2 O 3 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/cm 3 ) and Ref. [71] (DC-sputtered samples with density 3.51 g/cm 3 ).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 ∼ 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 roomtemperature 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 (κ(ρ) 300K ≈ a • ρ + b, where a=0.637±0.004 We compare our predictions with the experiments by Gorham et al. [73], who character-0 100 200 300 400 500 600 700 2.28 g/cm 3   2.98 g/cm 3   3.49 g/cm 3   Lee 3.51 g/cm 3  Monachon, 3.00 g/cm 3  Gorham, 3.00 g/cm ized the room-temperature thermal conductivity of ALD am-Al 2 O 3 grown on Si substrate [77] at densities ranging from 2.66 g/cm 3 to 3.12 g/cm 3 .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.

Thermal diffusivity and effects of anharmonicity
To gain microscopic, fundamental insights on why anharmonicity has negligible effects on the thermal conductivity of am-Al 2 O 3 , 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] where ω max is the maximum vibrational frequency of the system, g(ω) is the VDOS discussed in Sec.II B), C(ω, T ) is the specific heat for a vibration with frequency ω at temperature T (see Eq. ( 4)), and D(ω, T ) is the rWTE diffusivity [21], The expression for D qs is determined by factorizing the specific heat C qs in Eq. ( 3) and by the requirement that in the coupling between two vibrations qs and qs ′ each contributes to the coupling with a weight equal to the relative specific heat [40] (e.g. for vibration qs the weight is

Cqs
Cqs+C qs ′ , and correspondingly for vibration qs ′ the weight is C qs ′ Cqs+C qs ′ ).In the harmonic AF limit η≫Γ qs +Γ qs ′ →0, thus the Voigt distribution in Eq. ( 7) reduces to the Gaussian representation of the Dirac δ, 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(ω, T ), implying that differences between AF and rWTE conductivities derive exclusively from differences between the AF and rWTE diffusivities.
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 qs = C(ω qs , T )→k B and Γ qs ≫η ∀ qs).Specifically, the Voigt distribution in Eqs.(6,7) implies that in the high-temperature limit the rWTE diffusivity is a Lorentzian-weighted average of velocityoperator elements ∥v (q) s,s ′ ∥ 2 , and in practice such average is mainly determined by elements satisfying |ω qs −ω qs ′ |<[Γ qs +Γ qs ′ ]; increasing temperature implies an increase of the linewidths, and therefore velocityoperator 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 ℏ|ω qs −ω qs ′ | 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 (ℏ|ω qs −ω qs ′ | → 0) and therefore it does not depend on temperature.Applying this reasoning to am-Al 2 O 3 , 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 ℏω a = ℏ[ω qs +ω qs ′ ]/2 and difference ℏω d = ℏ|ω qs −ω qs ′ | of the modes coupled [21]: 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/cm 3 , computed from Eq. ( 8).We see that the average velocity operators for all densities are relatively unchanged at increasing differences, ω d .These small differences imply that the effects of anharmonicity on the conductivity are negligible (see text).
where G(ω a , ω d ) is a density of states that serves as normalization The plots confirm our expectations, i.e., in the am-Al 2 O 3 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 ω 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(ω) to simplify the discussion (i.e.Eqs.(6,7) evaluated with Γ qs = 0 ∀qs and η 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 ω is determined by their density of states g(ω), amount of heat carried C(ω, T ) and rate of diffusion D(ω).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.Analyzing the frequency-resolved AF diffusivity D(ω) reported in Fig. 8 for all our models of am-Al 2 O 3 allows us to address these questions.We see that an overall increase of diffusivity with density is visible when comparing low-density (ρ =2.28 g/cm 3 ), mediumdensity (ρ =2.98 g/cm 3 ) and high-density (ρ ≥ 3.17 g/cm 3 ) models, especially at low frequencies (from 0 to 100 cm −1 ).More precisely, comparing the diffusivity of the highest-density ρ = 3.49 g/cm 3 model with that of the lowest-density ρ = 2.28 g/cm 3 model, we find that the highest-density model has a diffusivity that is a factor ∼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(ω, 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.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 cm −1 .Inset, quantum specific heat as a function of 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-Al 2 O 3 , 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-Al  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/cm 3 model and the AF diffusivity of the 2.28 g/cm 3 model; (ii) using the VDOS of the 2.28 g/cm 3 model and the AF diffusivity of the 3.49 g/cm 3 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 ∼100 atoms, a result supported by other recent studies [47,79].In this section we study how the size of the atomistic model affects the theoretical predictions for the conductivity of am-Al 2 O 3 .To overcome the limitations posed by the computational cost of first-principles calculations, we generated a MACE MLP [34,35] for am-Al 2 O 3 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-Al 2 O 3 models containing 120, 3240, and 7680 atoms before finally computing the vibrational properties of these models, as well as the AF and rWTE conductivities.

A. Model generation via melt-quench
We used MACE to generate large models of am-Al 2 O 3 .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/cm 3 (hereafter referred to as AIMD120), we generated 3×3×3 and 4×4×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/cm 3 .We also relaxed with MACE the AIMD120 model; using the same threshold mentioned above and without any constraint on the ge-ometry of the box, obtaining a 'MACE120' model with density 2.88 g/cm 3 .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-Al 2 O 3 .Then, using the 120-, 3240-, and 7680-atom MACE models, we investigate how these properties depend on the size of the model.

B. Structural properties
We start by comparing coordination histograms of MACE structures with the AIMD120 structure.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.
tions 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 Total and partial RDF of the 120 atom Al2O3 model generated using AIMD at 2.98 g/cm 3 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.
(negligible) at distances larger than the linear size of our first-principles 120-atom models (∼ 11 Å, where the solid light-green and dashed dark-green lines stop).This suggests that atomistic models with linear size of ∼ 11 Å are sufficiently large to capture the most important features of structural disorder in am-Al 2 O 3 .Importantly, these tests also validate the capability of the MACE MLP to describe the structural properties of am-Al 2 O 3 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.

C. Vibrational properties
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.

D. 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.
In Fig. 13 we compare the bare WTE conductivities (i.e., Eq. (3) computed at q = 0 with η = 0, hereafter referred to as 'bWTE') against the regularized rWTE conductivity (computed relying on the q-mesh interpolation and using η 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 reg- 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 q=0 only); MACE120 and AIMD120 (both computed on a 5x5x5 q mesh).
ularization protocol has negligible effect on them.Importantly, the rWTE conductivity for the small AIMD120 and MACE120 models computed on converged 5x5x5 q-mesh and using η = 6 cm −1 yields a conductivity compatible [80] 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 A 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 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 antiperiodic boundary conditions employed by Feldman et al. [65] to accelerate computational convergence.To better understand the effect of the q-mesh interpolation, we show in Fig. 18 that the velocity-operator as a function of frequency for MACE120 on a 5 × 5 × 5 q mesh, approaches the velocity operator as a function of frequency computed at q = 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 ℏ∆ω avg , i.e. a broadening comparable to the value η used in the evaluation of the AF conductivity.In the context of Eq. ( 3), employing a broadening of the order of ℏ∆ω 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 η, 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 0 100 200 300 400 500 600 700 ) bWTE AIMD120 q = 0 bWTE MACE120 q = 0 bWTE MACE3240 q = 0 bWTE MACE7680 q = 0 rWTE AIMD120 q 5x5x5 rWTE MACE120 q 5x5x5 rWTE MACE3240 q 3x3x3 rWTE MACE7680 q 3x3x3 FIG. 13. rWTE and bWTE thermal conductivities of am-Al2O3 models at density ∼3.0 g/cm 3 .The solid lines (scatter points) are bWTE (rWTE) conductivities: MACE7680 is purple, MACE3240 is orange, MACE120 is blue and AIMD is green.
as am-Al 2 O 3 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.
To further validate our assertion that structural disorder is the dominant limiting factor for heat conduction in am-Al 2 O 3 over the range 50-700 K, we selected the MACE7680 model-which is large enough to directly compute the bulk conductivity, without regularizationand 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-Al 2 O 3 between 50 and 700 K, showing that in this system structural disorder dominates over anharmonicity 0 100 200 300 400 500 600 700 T (K) 0.0 0.5 MACE7680 q = 0 only MACE7680 q = 0 only -double linewidth MACE7680 q = 0 only -half linewidth FIG.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.
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-Al 2 O 3 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-Al 2 O 3 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-Al 2 O 3 from AIMD with densities ranging from 2.28 g/cm 3 to 3.49 g/cm 3 , describing how the atomic coordination topology varies with density.We have shown that at least five different atomic coordination environments coexist in am-Al 2 O 3 , 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-Al 2 O 3 , 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 theoryevaluated 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-Al 2 O 3 containing 3240 and 7680 atoms at density ∼ 3 g/cm 3 , 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 > ∼ 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.The 120-atom am-Al 2 O 3 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 q = 0 only is used as starting point.Therefore, to compute the vibrational properties of am-Al 2 O 3 on a mesh denser than q = 0 only and to use the most accurate densityfunctional perturbation theory (DFPT) technique [48], the vibrational properties are computed using Quantum ESPRESSO [84,85] on a 2 × 2 × 2 q-mesh (at present the DFPT implementation in VASP is restricted to calculations at q=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-Al 2 O 3 models were relaxed using a threshold for forces of 2×10 −4 Ry/Bohr and of 0.01 kBar for pressure (vc-relax command).

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 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.
FIG. 15.Convergence of the AF conductivity for am-Al2O3 with respect to the broadening η for the Dirac δ, The three upper panels show convergence plateaus for 5 × 5 × 5 (lines) and 7 × 7 × 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 q = 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 η 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 × 5 × 5 q mesh is dense enough to achieve computational convergence with respect to mesh size, since results obtained on a 7 × 7 × 7 mesh are practically indistinguishable from those obtained using a 5×5×5 mesh.Calculations at q = 0 only are far from computational convergence and consequently underestimate the conductivity.In order to calculate the bulk limit of the thermal conductivity of strongly disordered solids such as am-Al 2 O 3 , 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 η 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 A for details).
All the am-Al 2 O 3 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×5×5 or 7×7×7 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 q = 0 only for a 120-atom model is far from computational convergence and underestimates the thermal conductivity.The values of η that we determined from the convergence test discussed in Fig. 15 and that we employed in our calculations are reported in Table I. ρ (g/cm 3 ) 2.28 2.98 3.17 3.30 3.49 ℏη (cm −1 ) 6.0 6.0 6.0 6.0 8.34

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

Anharmonic Linewidths
The third-order interatomic force constants are computed in the 120 atom cells using ShengBTE [90]  The scatter points represent linewidths computed at q = 0, and the solid lines are coarse grained functions Γa[ω, 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 [40,88,89].The purple region denotes the overdamped regime Γ > ω [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 8 th 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 [33,[92][93][94][95]); (ii) the linewidths were computed considering exclusively the cubic terms in the Taylor expansion of the interatomic potential [21,40] and the contribution due to isotopicmass disorder [62].
Fig. 16 shows the linewidths for AIMD structures calculated at q = 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 Γ a [ω] that approximately describe the anharmonic linewidths as a single-valued functions of frequency, obtained following the approaches discussed in Ref. [40] (see also Refs.[47,88,89] for similar approximated treatments of anharmonicity).The approximated function Γ a [ω] 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].

rWTE conductivity calculation
The quantities needed to evaluate the rWTE conductivity (3) were computed as follows: (i) the η 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 q mesh; (iii) the linewidths were evaluated from the frequencies determined at the previous point using the function Γ a (ω) discussed in Sec.C 3.
We checked that increasing the size of the q mesh from 5×5×5 to 7×7×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 × 5 × 5 q mesh.We limited our calculations to temperatures higher than 50 K, since at temperatures lower than 50 K the thermal properties of am-Al 2 O 3 are dominated by lowfrequency vibrational modes that are likely to feature glassy anomalies [96][97][98]; accurately sampling these lowfrequencies vibrational modes requires using atomistic models containing thousands of atoms, and it is therefore beyond the scope of the present work.D(ω) data for the plot in Fig. 8 was calculated using Eq. ( 6) with δ-function smeared to a Gaussian with variance π 2 η 2 0 , with ℏη 0 = 8.34 cm −1 .
Appendix D: Thermal conductivity calculations using MACE 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 q-mesh that is Convergence of the AF conductivity for MACE models of am-Al2O3 with respect to the broadening η for the Dirac δ, 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 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 η 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 η 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 q-mesh or q=0 point only, and using larger models extends the plateau to lower values of η, as in vitreous silica [21].
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 q=0 point only that is compatible with the plateau at 3x3x3 q-mesh indicating achievement of computational convergence.
Given that Fig. 14 shows that anharmonicity has negligible effects on the conductivity of am-Al 2 O 3 , 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/cm 3 AIMD model (central panel in Fig. 16).
Appendix E: Velocity operator for MACE and AIMD models at densities close to 3.0 g/cm 3   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 (⟨|ν avg ωaω d | 2 ⟩, 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 q=0 and 5x5x5 q-mesh calculations of velocity operator elements for 120-atom MACE model.We see that q=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 q-mesh, with those at q=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., 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 q-mesh.The middle panel shows that for MACE120, the velocity operator at q=0 (dashed) is overall smaller than the velocity operator computed over a 5x5x5 q-mesh (solid).The bottom panel shows that the velocity operator for MACE120 computed over a 5x5x5 q-mesh (solid) is in reasonable agreement with the velocity operator of MACE7680 computed at q=0 (dashed).The delta functions used for calculation of ⟨|ν avg ωaω d | 2 ⟩ were replaced with Gaussian with variance related to the smearing parameter from the convergence plateau: σ 2 = η 2 π/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-Al 2 O 3 at densities close to 3.0 g/cm 3 , influence of anharmonicity on the magnitude of conductivity is very weak, which is exemplified by very good agreement between WTE q=0 and AF 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 q-mesh in small models give the same trend and very similar magnitude as predictions of AF conductivity using q=0 in large models for very disordered solids, which am-Al 2 O 3 is an example of.The largest differences of value of AF conductivity at 700K are between MACE7680 at q=0 and on 3x3x3 q-mesh , and is equal to approximately 0.04 W/mK (2 %). ) bWTE MACE3240 q = 0 bWTE MACE7680 q = 0 AF MACE3240 q = 0 AF MACE7680 q = 0 AF AIMD120 q 5x5x5 AF MACE120 q 5x5x5 AF MACE3240 q 3x3x3 AF MACE7680 q 3x3x3 FIG.19.AF thermal conductivity for models of am-Al2O3 with densities close to 3.0 g/cm 3 .Solid (dashed) lines are AF (WTE) calculations done at 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).

FIG. 1 .
FIG. 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.
FIG. 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/cm 3 , computed from Eq. (8).We see that the average velocity operators for all densities are relatively unchanged at increasing differences, ω d .These small differences imply that the effects of anharmonicity on the conductivity are negligible (see text).
FIG.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 cm −1 .Inset, quantum specific heat as a function of temperature.

2 O 3
at technologically relevant temperatures.
FIG. 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.

FIG. 11 .
FIG. 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/cm 3 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.

98 FIG. 12 .
FIG.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 q=0 only); MACE120 and AIMD120 (both computed on a 5x5x5 q mesh).
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 EP-SRC 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 C: First-principles thermal conductivity calculations 1 .
Convergence of the Allen-Feldman theory

FIG. 16 .
FIG.16.Effect of temperature on anharmonic linewidths of am-Al2O3 at various densities and temperatures.The scatter points represent linewidths computed at q = 0, and the solid lines are coarse grained functions Γa[ω, 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[40,88,89].The purple region denotes the overdamped regime Γ > ω[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].
FIG. 17. Convergence of the AF conductivity for MACE models of am-Al2O3 with respect to the broadening η for the Dirac δ, 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 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 η 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 η 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 q-mesh or q=0 point only, and using larger models extends the plateau to lower values of η, as in vitreous silica[21].

1 MACE120 2 .FIG. 18 .
FIG.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 ⟨|ν avg ωaω d | 2 ⟩ is represented as a function of energy differences ℏω d = ℏ(ω(qs) − ω(q s ′ )) and averages ℏωa = ℏ ω(qs)+ω(q s ′ ) 2 [51] , and Al 6 are similar in shape.We note that this is in sharp contrast to the crystalline θ-Al 2 O 3 phase, which contains only Al 4 tetrahedra and Al 6 octahedra and displays a clear distinction between high-frequency breathing modes associated with Al 4 environments, and low-frequency bending modes associated with Al 6 environments[51].Varying the density of am-Al 2 O 3 causes a variation in the relative magnitude of the Al x PDOS: at 2.28 g/cm 3 , the Al 4 PDOS is always larger in magnitude than Al 5 ; with increasing density, the relative magnitude of the PDOS of Al 4 and Al 5 progressively reverses, with the highest-density 3.49 g/cm 3 model featuring a PDOS for Al 5 always larger than that for Al 4 .We note that the O 3 and Al 5 coordination environments have a significant PDOS in all the models analyzed, regardless of the models' density.In contrast, the PDOS of O 2 , O 4 , Al 4 , Al 6 , display much stronger variations with density-increasing density yields a progressive suppression of the O 2 and Al 4 environments, compensated by the progressive appearance of O 4 and Al 6 environments.
[78]mal 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]ts: 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.[71]in1995 (red triangle, DC sputtering).
Microscopic mechanisms underlying conductivity increase with density.Solid lines are exact AF conductivity calculation for the highest-density 3.49 g/cm 3 model (orange, g3.49D3.49)andlowest-density2.28 g/cm 3 model (cyan, g2.28D2.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, g3.49D2.28),or the VDOS of the lowest-density model and the diffusivity of the highestdensity model (grey, g2.28D3.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.

TABLE I .
Broadening parameters η used for the Gaussian representation of the Dirac δ function appearing in the AF conductivity expression, and for the Voigt distribution appearing in the rWTE expression.