Hydrodynamic finite-size scaling of the thermal conductivity in glasses

In the past few years, the theory of thermal transport in amorphous solids has been substantially extended beyond the Allen-Feldman model. The resulting formulation, based on the Green-Kubo linear response or the Wigner-transport equation, bridges this model for glasses with the traditional Boltzmann kinetic approach for crystals. The computational effort required by these methods usually scales as the cube of the number of atoms, thus severely limiting the size range of computationally affordable glass models. Leveraging hydrodynamic arguments, we show how this issue can be overcome through a simple formula to extrapolate a reliable estimate of the bulk thermal conductivity of glasses from finite models of moderate size. We showcase our findings for realistic models of paradigmatic glassy materials.


INTRODUCTION
The theory of thermal transport in crystalline and disordered solids has witnessed a major advancement in the past few years [1][2][3][4][5][6][7]. Until recently, the only way to simulate heat transport in glasses and complex crystals [8] was through the sheer application of the Green-Kubo (GK) linear-response theory, leveraging the heatflux time series generated by Molecular Dynamics (MD) simulations [9][10][11][12][13]. Among the many merits of this approach is the full account of anharmonic effects of any order; however, the same cannot be said of quantum effects, which may be important at low temperatures, but cannot be captured by any classical method. The quest for a more convenient way of computing the thermal conductivity in glasses and complex crystals led to the unification of the approaches based on the Boltzmann Transport Equation (BTE) for crystalline materials [14,15], and the Allen-Feldman (AF) model of thermal conduction in harmonic glasses [16,17], resulting in the quasi-harmonic GK (QHGK) [2] and Wignertransport-equation (WTE) [3] approaches to heat conduction. The former method consists in an explicit evaluation of the GK expression of thermal conductivity from the anharmonic lattice dynamics of a solid close to mechanical equilibrium [2,5,18], while the latter exploits a different transport equation based on Wigner's quantum dynamics [3,4,6]. Despite their differences, these two methods provide comparable results and are in fact conceptually equivalent, both being approximations of the same order to a general many-body approach based on GK linear-response theory and the Mori memoryfunction formalism [5]. These methods suggest that the behavior of thermal conductivity in anharmonic solids can be understood in terms of the interplay between a quasi-particle kinetic mechanism, whereby traveling phonons propagate and scatter across many interatomic distances, and a diffusive regime, where disorder breaks * Mail should be sent to afiorent@sissa.it the quasi-particle nature of phonons, and heat is transported through short-range hopping mechanisms.
These advances notwithstanding, the numerical simulation of heat transport in glasses remains a formidable task, because it requires the use of large finite models for the bulk systems, comprising up to several thousand atoms [2,7,19,20]. The same would be true even for crystalline materials, the difference being that, in such a case, lattice periodicity and the Bloch theorem can be exploited to map the simulation of a large model onto a number of computations performed for the vibrations of definite wavevector in the Brillouin zone (BZ), k, of a unit cell comprising a much smaller number of atoms. These wavevectors are usually arranged in a regular grid whose number of nodes is the ratio between the total number of atoms in the bulk model and the number of atoms in the unit cell. The effective model that can thus be afforded has a linear dimension L ∼ 2π/k min , k min being the discretization step of the regular grid, and it may encompass up to several hundred thousand atoms. Unfortunately, the spurious crystalline order introduced by the use of periodic boundary conditions (PBCs) in the simulation of finite glass models results in unphysical long-wavelength features in their transport properties [7], as it will be extensively illustrated in the following.
These issues call for a method able to efficiently and accurately extrapolate to infinite size the value of the thermal conductivity of aperiodic solids, such as glasses, without the need to artificially introduce nonphysical normal modes, which may lead to a gross overestimate of the final result. We devise such a method by leveraging arguments based on the hydrodynamics of solids, which exploit the natural partition of glassy normal modes into categories based on their localization in real space [21]: propagons are delocalized, lowfrequency, vibrations whose typical wavelength is much larger than the interatomic distance and that therefore propagate almost freely as plane (sound) waves; diffusons are intermediate-energy vibrations that spread diffusively; finally, locons are localized excitations that populate the highest-frequency portion of the vibrational spectrum and they hardly contribute to the transport arXiv:2303.07010v2 [cond-mat.mtrl-sci] 1 Sep 2023 properties of the system. Being long-wavelength vibrations, propagons are severely affected by the finite size of any glass model. In fact, the size of the system sets a lower bound to the minimum frequency accessible to the calculation, and determines the normal-mode spacing in the low-frequency region, which is already undersampled with respect to other portions of the spectrum, due to the vanishing of the Vibrational Density of States (VDOS) in the zero-frequency limit [2,18,22]. In this article we report on a method, which we dub hydrodynamic extrapolation, able to lift these problems through a combination of two ingredients that can be inexpensively computed from finite models of moderate size: one is the QHGK contribution of diffusons and locons [2,5,18], the other being an effective low-frequency model for propagons. The term "hydrodynamic" is more often used in the context of the collective phonon transport in crystals [23]. Here, we adopt this term in a similar vein, as the underlying equations governing both crystals and glasses at low frequencies obey the same hydrodynamic principles [24].
We underline that the hydrodynamic extrapolation applies above the so-called plateau of the thermal conductivity as a function of temperature which is found in many glassy materials. The temperature dependence of κ in glasses features three universal characteristic trends [25]: at very low temperatures, i.e. T ≲ 2 K, the dominant scattering mechanism is the quantum tunneling between different local minima in the glass energy landscape, and κ ∼ T 2 [26][27][28]; then, up to ≈ 30 K, the thermal conductivity rises and saturates to a plateau value. Despite the absence of an established theoretical agreement in the literature, this phenomenon seems to be related to the crossover between the regime dominated by quantum processes and one where propagating waves are scattered by random disorder [25,[27][28][29]. Above the plateau, the behavior of κ is dictated by the anharmonic decay of nomal modes as prescribed by the QHGK theory. Since we do not have access to quantum-tunneling states that play a crucial role at low temperatures, in this work we focus on the third regime.
The structure of this Article is the following: first, we briefly review the QHGK method for computing the thermal conductivity of crystals and solids alike; then, we discuss the possibility of simulating heat transport in glasses from finite models in PBCs; we then introduce our hydrodynamic extrapolation technique, based on an effective low-frequency model for propagons; next, we illustrate our newly introduced method with a few numerical experiments performed for three paradigmatic classes of glass (amorphous silicon carbide, silica, and silicon); finally, we draw our conclusions. Following is a thorough description of the theoretical methods employed in our work.

RESULTS AND DISCUSSION
The Quasi-Harmonic Green-Kubo Approximation The quantum GK formula for the thermal conductivity tensor, κ, reads [30][31][32]: where V is the system's volume, k B is Boltzmann's constant, T is the temperature, J is any Cartesian component of the energy-flux operator in the Heisenberg representation, and isotropy is assumed throughout this paper, thus allowing us to dispose of Cartesian indices, unless strictly needed for clarity. Both crystalline and amorphous materials are normally described by finitesize models in PBCs, and normal modes can be labeled by wavevectors and band indices. Any Cartesian component of the harmonic energy flux operator reads [2,33]: where (q, ν) labels wavevector and band indices in the BZ, respectively;â † qν andâ qν are the creation and annihilation operators, respectively; ω qν is the angular frequency of the normal mode; v qνν ′ is a generalized group velocity [2,6] (see Eq. (31)) which satisfies Here again, the Cartesian index of the velocity is dropped for notational simplicity. We stress that, while for crystals the classification of normal modes by their wavevectors reflects a physical symmetry of the systems, in glasses it is simply an artifact due to the use of PBCs. When large enough models are used to describe the glass, it is common practice to just sample the BZ center (Γ, q = 0), although ways of leveraging the phonon dispersions resulting from the use of PBCs in glasses have recently been proposed [7]. The first, number-conserving term, in Eq. (2) is called resonant, while the second, anti-resonant, one gives negligible contributions in the small-linewidth limit and will be neglected in the followig. The diagonal term of the generalized velocity is the usual phonon group velocity: v qνν = ∇ q ω qν . Eqs. (1) and 2 lead to the QHGK approximation to the thermal conductivity [2,5,18] where is the generalized two-mode heat capacity, n qν = (e ℏωqν /kBT − 1) −1 the Bose-Einstein distribution, and is the generalized two-mode lifetime. In a nutshell, the QHGK approach entails two different-yet relatedapproximations, The first is the dressed bubble approximation of four-point correlation functions among creation and annihilation operators; it amounts to neglecting vertex corrections to factorize four-point correlation functions into sums of products of two-point correlation functions [5]. This means that normal modes decay independently of one another due to effective interactions with a common heat bath. The second approximation, often referred to as Markovian, is to ignore any memory effects in the interaction between the heat bath and the normal modes [5]. The combination of these two approximations in the QHGK approximation is that four-point correlation functions can be expressed in terms of single-body greater Green's functions, given by g > qν (t) = −i(n qν + 1)e −iωqν t−γqν |t| , where γ qν is the linewidth of the normal mode labeled by (q, ν) as computed from Fermi's Golden Rule. The quasi-harmonic hypothesis requires γ 2 qν /ω 2 qν ≪ 1, implying that only almost-degenerate modes such that |ω qν − ω qν ′ | ≲ γ qν + γ qν ′ are allowed to contribute to the heat conductivity.
Glasses are, by definition, aperiodic. Yet, the effect of periodicity on large-enough models can be regarded as a size effect that vanishes in the thermodynamic limit. It is thus customary to treat glasses as finite models in a large supercell under PBCs sampled only at the zone center of the BZ. Omitting the sum over the (only) wave vector, Γ, and assuming the system is isotropic (thus, the thermal conductivity tensor is proportional to the identity matrix) the thermal conductivity reads: The factor 1/3 comes from the average over the three Cartesian directions, which are equivalent under the hypothesis of isotropy. At the Γ point, the eigenvectors of the dynamical matrix can be chosen to be real, and v νν ′ becomes a real, antisymmetric matrix [2,18]. The harmonic limit is obtained when all the linewidths approach zero. In that case, the QHGK method reduces to the AF model [16,17]: where C ν is the modal specific heat and the modal diffusivity, D ν , is defined as The Dirac-delta function appearing in the diffusivity is due to the harmonic, limit lim γν ,γ ν ′ →0 τ νν ′ , of the Lorentzian generalized lifetimes defined in Eq. (6). For crystalline materials, using a dense mesh in the BZ (q-mesh) allows one to efficiently increase the number of modes, and it is often required to get converged values of physical quantities, a fortiori when the simulation cell is a (usually small) unit cell. In model glasses, the lowest frequency allowed in a cubic simulation cell with side L is of order 2πc/L, c being the sound velocity; sampling a q-mesh mimics an artificial order by adding long-lived modes below this threshold, thus reflecting the periodicity induced by PBCs and determining a spurious increase of the thermal conductivity. Indeed, real glasses below this threshold frequency feature propagating waves whose wavelength is proportional to their inverse frequency, and whose decay depends on disorder [21]. A replicated model glass cannot be disordered on a spatial scale larger than the original simulation cell. Thus, the spurious modes whose wavelength is larger than the simulation cell are unaffected by disorder and would propagate for longer times than the corresponding modes in an infinite, actually disordered, glass. In fact, finite wavevectors other than high-symmetry points at the BZ boundary have, in general, a finite group velocity, resulting in The first term in the above equation is the same obtained by treating a crystal in the Relaxation Time Approximation (RTA) of the Boltzmann Transport Equation (BTE). When only third-order anharmonic effects are included, this term gives rise to the 1/T trend of κ in crystals at high temperatures. Sampling the BZ of a glass model induces the same behavior, as recently observed in Ref. 7, where a 3 × 3 × 3 q-mesh is shown to determine a lowtemperature divergence of the thermal conductivity of amorphous silica systems. The same does not happen for a genuinely disordered system with an equivalent number of atoms (i.e., 3 3 times the number of atoms of the small model) sampled at the Γ point; for such a system the thermal conductivity is a monotonically increasing function of temperature [7]. Thus, replicating a small disordered system by sampling its reciprocal space is qualitatively different from simulating a correspondingly large system in real space, since in the former case disorder is suppressed at scales larger than its size. This behavior was also observed in GKMD calculations on glass samples replicated in real space [34,35].

Hydrodynamic extrapolation of the thermal conductivity
The naive approach to deal with finite-size effects is to compute the heat conductivity for model systems with increasing size, L, up to convergence, κ ∞ (T ) = lim L→∞ κ L (T ). The simplest way to extrapolate κ ∞ (T ) would then be to assume that κ L can be written as a power series in 1/L and perform a linear fit in 1/L. This approach has two key issues: first, the whole procedure is computationally expensive, since it requires the preparation of samples with different sizes-at the very least three, in order to meaningfully fit a straight line, but in practice we observe that more points are usually needed; second, in cases where L is too small for the linear contribution to be dominant, higher powers of 1/L might be necessary, thus increasing the number of points required to perform a reliable fit. In the following, we develop a method to bypass these obstacles making use of the features of the low-energy excitations naturally present in glasses.
As mentioned above, propagons, diffusons, and locons differ by the degree of localization they feature. This can be observed in the Dynamical Structure Factor (DSF) that, for a harmonic system, is defined as [35][36][37]: where the scalar product between normal modes and plane-wave states is α being the Cartesian component; Q = 2π L (n, m, l), with n, m, l ∈ N, is a wavevector in a cubic supercell of side L; R I the position of the Ith atom; ϵ ν the eigenvector of the νth normal mode; and ε b (Q) a polarization (unit) vector (see Methods). The latter can be chosen to be parallel or perpendicular to Q and can be labeled by b = L, T 1 , T 2 , the longitudinal (parallel to Q) and transverse (perpendicular) branches, respectively. Transverse branches are degenerate, so we indicate with T the contributions of both transverse branches. The DSF can be generalized to take into account weak anharmonic effects as where the Dirac-delta in Eq. (11) is replaced by a Lorentzian function centered on the mode's angular frequency, whose spread is given by the mode's linewidth. Eq. (13) gains a temperature dependence through the anharmonic linewidths with respect to Eq. (11). The low-frequency, small-wavevector, portion of each branch of the DSF features an almost linear dispersion, ω ≃ cQ, typical of acoustic phonons. In other words, S b (Q, ω) is a peaked function centered at ω Q = c b Q, c T,L being the transverse/longitudinal speed of sound. For low-enough Q, the profile of S b (ω, Q) can be fitted as a function of the angular frequency with a Lorentzian profile, allowing one to evaluate the speed of sound as well as the wavevector dependence of the effective width, Γ b (Q). Under the assumption of isotropy, Propagons are identified as those low-frequency, low-wavevector, normal modes with linear dispersion that populate the first portion of the DSF. The increasing broadening of the dispersion identifies a cutoff frequency for propagons, ω P , referred to as the Ioffe-Regel limit [38]. The peaked shape of S b (ω, Q) suggests the |Q, b⟩s are long-lived excitations. In the long wavelength limit, they form an almost-orthonormal basis even for aperiodic materials such as glasses [36,39], whose timeevolution for positive times-as inferred from the DSF- Unlike the case of normal modes, which are eigenvectors of the harmonic (disordered) Hamiltonian, plane-wave states decay both through anharmonic interactions and as an effect of harmonic disorder. In fact, even the harmonic DSF, S 0 b , features a finite width for any Q [36] without the need for any anharmonic effect. This does not happen in a periodic system. In order to extrapolate the thermal conductivity of larger systems, we separate the contribution of propagons from that of from diffusons and locons: The first contribution, due to propagons, involves only normal modes below the cutoff frequency, ω P . The second contribution involves normal modes above such frequency, diffusons and locons, and it is mainly due to diffusons. The propagon contribution can be approximately written on the basis of acoustic plane waves |Q, b⟩ [36,40]. The overall idea is to switch from the basis of normal modes, whose decay can only be anharmonic, to a basis of delocalized modes whose decay, encompassed in the definition of Γ b (Q), entails a combination of anharmonicity and disorder. We demonstrate this in the Methods section, while giving here just a sketch of the proof. Let us consider the energy flux of Eq. (2); the resonant propagon contribution to the energy flux reads: whereâ and J bb ′ Q is the matrix element of J in the plane-wave basis. Similarly to what happens in crystals, such matrix elements can be shown to have a block-matrix structure in the wavevector indices The two-point correlation function of the acoustic waves is related to the anharmonic DSF by that can be rewritten as The propagon contribution to Eq. (1) can thus be evaluated from Eq. (20), resulting in The last equation defines the acoustic excitations' life- , and the frequency-dependent heat capacity, C(ω) = ℏω ∂nω ∂T . In principle, an interband contribution should also appear between the transverse and longitudinal branches, with a Lorentzian weight as in Eq. (4). However, longitudinal and transverse branches are energetically well separated, i.e., [41], and because propagons are, by definition, modes with a sharp linear dispersion relation, i.e., c b Q ≫ Γ b (Q). Eq. (21) is reminiscent of the RTA for crystals, and a similar equation appeared many times in the literature on harmonic glasses (e.g., in Refs. 36, 37, and 42), since it is just a statement that, at large-enough length scales, glasses can be effectively described by a continuous elastic model. Here, this fact is derived from first principles (i.e., from the QHGK theory) also including anharmonic effects.
The final step to evaluate the infinite-size limit of κ is to map Eq. (21) onto the kinetic theory of gases. The sum over plane waves can be recast as an integral over angular frequency in the L → ∞ limit using the fact that, for propagons, the linear dispersion relation, ω = cQ, is valid, and introducing per-branch densities of states, . This yields: Acoustic excitations feature a density of states of the Debye form, i.e., ρ T (ω) = ω 2 Once κ P is available, one computes the hydrodynamic extrapolation of the thermal conductivity as where the second term comes from Eq. (7), and the Heaviside-theta limits the sum to the nonpropagonic part of the spectrum. The choice of ω P must satisfy certain requirements. For the model to be valid, all the normal modes below ω P should be bona fide propagons for both polarizations, so as to avoid the mixed transverse-diffusons/longitudinalpropagons regime. Thus, even though the IR limit is different for the two polarizations, ω P should be less than the smallest of the two, which is generally the transverse one. Slightly lower values can also be preferred to guarantee that all the hypotheses regarding linear dispersion are valid. The SI includes a dedicated section that provides practical suggestions on selecting the appropriate value of ω P . Another important point is how to get a good estimate of Γ b as a function of ω, which is needed to compute κ P from Eq. (22). This is done by first fitting the computed DSFs, Eq. (13), as a function of ω, at fixed Q, to a Lorentzian function, Eq. (14). We then use the linear dispersion, ω b = c b Q, to get the functional dependence of of the linewidth upon frequency, Γ b (ω). We finally fit this dependence to a quartic polynomial, Γ b (ω) = aω 2 + bω 4 with a, b > 0, as it can be observed in Fig. 1 for the linewidths of aSi. A model of this form was proposed also in Ref. 43, where the quadratic and quartic terms come, respectively, from the Umklapp and isotopic scattering of long-wavelength phonons in crystalline silicon. More generally, the Γ ∼ ω 2 trend is required by hydrodynamics [24,44], while the Γ ∼ ω 4 behavior-also found in experiments [45]-can be rationalized in the continuous limit through random media theory [46,47], or for atomistic models through harmonic perturbation theory, such as in the case of crystals with mass disorder [48,49] or random spring constants [50]. The crossover between the ∼ ω 2 and ∼ ω 4 behaviors in the sound attenuation can be understood analyzing nonaffine displacements in amorphous solids [51]. A rule of thumb to assess whether the bandwidth can be meaningfully computed from a glass sample of given size is to check whether the minimum frequency available at that size, ω min , is smaller than the crossover frequency between anharmonicity and disorder; namely, ω min ≲ a/b.
In order to validate Eq. (23), we compare its predictions with those of an improved version of the naive procedure outlined at the beginning of this section. Namely, instead of a polynomial fit of κ L (T ) in 1/L, which can easily incur overfitting due to the small number of data points usually at one's disposal, we opt to use a different functional form suggested by the continuous model we are using: where A and B are positive parameters. This ansatz can be derived from Eq. (22) under the following assumptions: i) that the classical limit is valid [i.e., C(ω) = k B )]; ii) that the form of the linewidth is Γ(ω) = aω 2 +bω 4 ; and iii) that the only size effect comes from the value of the minimum frequency available at each size, ω min = 2πc/L. In principle, one would need two sets of parameters, one for each polarization, but this would unnecessarily complicate the fitting procedure. Using just one set of parameters, as in Eq. (24), amounts to saying that we are effectively averaging the contributions due to each polarization. Under these simplifying hypotheses, Eqs. (22) and (23) become: Since only low-energy modes with 0 ≤ Q ≤ 2πc/L are involved, the classical approximation of the modal specific heat is further justified. This ansatz not only reproduces the 1/L trend for large-enough sizes, but it also captures the non-linearity (in particular the convexity) of our data. Fig. 2 summarizes how the hydrodynamic extrapolation scheme outlined above is applied to estimating the bulk thermal conductivity of glass samples.

Numerical experiments
We now showcase the ability of our method to accurately predict the bulk limit of the thermal conductivity for three paradigmatic glasses, chosen because they have been extensively investigated in the literature and because they are representative of different regimes as concerns the role of propagating modes: amorphous silica (aSiO 2 ) is a material whose thermal transport properties are mainly determined by diffusons; in amorphous silicon (aSi) transport is dominated by propagons; the situation is somewhat intermediate in amorphous silicon carbide (aSiC) [37,52].
For each material, the heat conductivity is averaged over ten different samples obtained through a melt-andquench procedure described in the Methods section.
Second-and third-order interatomic force constants as well as normal-mode frequencies and lifetimes have been computed with the κALDo [18] code, using forces obtained by the LAMMPS MD code [53].
For notational convenience, the size-scaling data reported in the following are plotted against the number replicas of the fundamental simulation cell repeated along each Cartesian component, ℓ = L/L 0 , where L is the actual linear size of the simulation cell being employed and L 0 is the size of the smallest cell used for each material, which contains 8, 24, and 8 atoms for aSiC, aSiO 2 , and aSi, respectively.

Amorphous silicon carbide
The top panel of Fig. 3 shows the DSF of aSiC. The dispersion for the longitudinal branch is linear up to approximately 8.5 THz, in agreement with the IR limit found in Ref. 52. From the slope of the linear dispersion, the speeds of sounds are computed to be c L ≈ 8.0 km s −1 and c T ≈ 4.4 km s −1 . We choose a cutoff ω P /2π = 3 THz, in order to stay well within the linear dispersion regime for both polarizations. In the left panel of Fig. 4, we report the resulting κ hydro for samples ranging from 4096 atoms (ℓ = 8) to 13824 atoms (ℓ = 12). Size effects are quite important at 100 K, where the bulk conductivity is almost 35% larger than one of our largest models. The hydrodynamic extrapolation performs well with relatively small samples. At higher temperatures, size effects are less striking; nonetheless, the hydrodynamic extrapolation still improves results with respect to plain QHGK calculations.

Amorphous silica
The method is tested also on amorphous silica, with a ω P /2π = 1.7 THz, which is quite low. The reason for this is that, as shown in the middle panel of Fig. 3, the IR limit for longitudinal waves is reached for frequencies as small as ≈ 2 THz, which is quite smaller than those of amorphous silicon and silicon carbide, whose values are ≈ 8.5 THz and ≈ 9 THz, respectively [52]. This fact suggests that the contribution of propagons to the thermal conductivity of aSiO 2 is minor, and that size effects are rather small. Our calculations confirm these statements: the central panel of Fig. 4 shows that not only is the QHGK thermal conductivity of the 8000-atoms model already large enough to be practically at convergence in size, but that its dependence on inverse size, 1/ℓ, is practically linear for any meaningful smaller size: thus, in this case, there is not much to gain in computing hydrodynamic corrections for aSiO 2 , if compared with other materials.

Amorphous silicon
aSi features a sharp linear dispersion in the lowfrequency region where the propagon linewidth is considerably smaller than in aSiC and aSiO 2 , as one can see from the bottom panel of Fig. 3. This reflects on the large value of the propagonic contribution to the thermal conductivity: in fact, with ω P /2π = 3 THz, κ P constitutes about 60% of the thermal conductivity at room temperature, a fraction that inevitably grows for lower temperatures, due to the increasing importance of lowenergy modes. The consequence of the dominance of the propagons in the value of κ is that size effects become rather crucial in aSi. In particular, the thermal conductivity has a conspicuous nonlinear dependence on 1/ℓ in the size range accessible to numerical simulations. Since such nonlinearity is mainly due to propagons, it is well captured by the fitting procedure discussed above.
In order to verify our results are not due to overfitting, we performed some classical GKMD at high temperatures on systems with larger sizes, up to 64000 atoms (ℓ = 20), which are not affordable through lattice methods such as QHGK. As reported in the Supplementary Information (SI) [39], even with 64000 atoms κ is still far from convergence; however, the computed GKMD thermal conductivities are compatible with the fit on the QHGK data, thus independently confirming the validity of the whole procedure. Notwithstanding, the hydrodynamics extrapolation on samples larger or equal than ℓ = 8 (N = 4096) leads to bulk thermal conductivity values compatible with the fitted ones, as shown in the right panel of Fig. 4.

Comparison with existing methods
A different approach has been recently proposed to extrapolate the value of the bulk thermal conductivity in glasses from simulations performed for small models and leveraging the symmetry properties of the fictitious crystal resulting from the periodic replica of the simulation cell in PBCs [7]. The proposed method consists of two steps: a finite model for a glass under PBCs is first treated as a genuine crystal, i.e., its vibrational properties are sampled on a dense q-mesh in the BZ: we call this pro- cedure a crystalline model for the glass. The 1/T divergence of the heat conductivity, resulting from crystalline order, is then avoided by convolving the Lorentzian lineshape of the unphysical low-frequency modes introduced by the periodic replicas of the fundamental simulation cell with a smearing function whose width, η, is carefully chosen as described below. While choosing the smearing function to be Lorentzian would result in a Lorentzian smeared lineshape whose width would be the sum of the original and smearing linewidths, it was found that using a Gaussian smearing function would lead to a better numerical behavior of the procedure. Of course, the heat conductivities computed with this method depend on the smearing width, and they would in fact diverge, in the low-temperature limit, for vanishing widths. In Ref. 7 the choice of the smearing width relied on the existence of a plateau in the dependence of the heat conductivity on it. When such a plateau exists, the width so chosen should be a compromise between a large enough value allowing to encompass several normal modes and small values necessary not to mimic spuriously large scattering sources.
In a nutshell, this protocol forces the QHGK (or, equivalently, WTE) method for a nonphysical crystal-whose unit cell is the whole simulation representing the glass model-to behave like the AF model for low-frequency modes whose linewidths are smaller than the smearing width-i.e., the modes responsible for the unphysical overestimate of κ-while doing almost nothing to modes with higher frequency/larger linewidth.
This method was applied to amorphous silica (aSiO 2 ) [7]. A computationally affordable aSiO 2 crystalline model treated with a 3 × 3 × 3 q-mesh, once regularized with the above procedure, yielded values of κ compatible with the temperature-dependent thermal conductivity of a genuinely larger system with an equivalent number of atoms, thus validating the whole procedure for aSiO 2 . While we are able to reproduce these results and find a plateau in the dependence of κ on the smearing width for our aSiO 2 samples, we cannot do so for other materials, where no plateau is found, thus making it difficult in general to fix the value of a suitable smearing width. In the SI [39] we test the regularization procedure on amorphous silica, amorphous silicon carbide, and amorphous silicon.
We rationalize these difficulties as follows. The lowtemperature behavior of the heat conductivity of bulk solids is brought about by the interplay of anharmonicity and disorder. The former gives rise to a ω 2 dependence of Γ b at low frequencies and dominates the lowest-frequency portion of the propagon spectrum, while the effects of the latter, contributing as ω 4 , emerge later and dominate at larger frequencies. Sampling a periodically repeated cell of linear size L on regular q-mesh completely misses the effect of disorder for the modes introduced by periodicity of frequency smaller than ∼ 2πc/L. The convolution of the vibrational lineshapes with a smearing function mimics the introduction of an effective boundary-scattering term in the form of a constant linewidth, rather than the correct frequency-dependent scattering due to disorder. The net effect is that the thermal conductivity is forced to quickly go to zero for vanishing temperatures. For materials where the contribution of propagons to the heat conductivity is small, such as aSiO 2 [37], the overdamping of low-frequency modes does not constitute a problem; vice versa, it fails to capture the main damping mechanisms that abate the heat conductivity in the lowtemperature regime in materials where the propagons' contribution is significant, if not dominant. More details on calculations based on Ref. 7 can be found in the SI [39].

Temperature dependence of thermal conductivity
It is instructive to analyze the temperature dependence of aSi since it is the material, among those at our disposal, where propagons contribute the most to the value of the thermal conductivity. The upper panel of Fig. 5 compares the hydrodynamic extrapolation with the QHGK calculations performed on progressively larger systems. The lower panel provides a breakdown of the extrapolated κ into the propagon contribution, κ P , and the diffuson one, κ D . The contribution of diffusons resembles the one of a standard AF calculation [21,36], ranging from small values at low temperatures to some constant value around room temperature. Conversely, the behavior of κ P is rather surprising, as it decreases with increasing temperature. This trend cannot be explained by a purely harmonic theory, where the temperature dependence is solely determined by the heat capacity [2,17], and it is thus an anharmonic effect.
The comparison between the hydrodynamic extrapolation and QHGK calculations highlights the significance of the former, especially at lower temperatures, as illustrated in both the upper panel of Fig. 5 and Fig. 4. This observation can be attributed to two key factors: (i ) at low temperatures, only propagons are populated and contribute significantly to thermal conductivity. Consequently, it becomes essential to employ a continuous model that accounts for frequencies below the minimum frequency allowed by the finite-size sample. (ii ) Numerically accurate QHGK results require the frequency spacing between modes to be smaller or comparable to the anharmonic linewidths, as specified in Eq. (7). As the anharmonic linewidths decrease with temperature, achieving meaningful QHGK thermal conductivity values at low temperatures necessitates a denser VDOS, i.e., larger systems. The lower panel of Fig. 5 also includes experimental measurement of the thermal conductivity of aSi on samples of relatively large sizes obtained with chemical vapor deposition techniques. The data are sourced from Liu et al. [54], Yang et al. [55], and Wieczorek et al. [56]. The respective films have thicknesses of 80 µm, 1.6-2.8 µm, and 2-3.6 µm. The motivation behind selecting these larger samples from the literature is to minimize size effects and obtain an estimate of the bulk thermal conductivity of aSi.
The experimental values of κ exhibit good agreement among themselves above approximately 200 K. However, at lower temperatures, Wieczorek et al. [56] reports slightly smaller values compared to the other sources. Overall, the behavior of the thermal conductivity is well captured by κ hydro , including the relatively high value observed around 100 K. Nevertheless, our model appears to slightly underestimate the thermal conductivity across the entire temperature range. We identify two primary reasons for this discrepancy.
One factor is the choice of the force field, specifically the Tersoff force field, which is not explicitly designed to accurately reproduce thermal conductivity. As a result, achieving numerical agreement with experimental results using this force field may be challenging. To address this issue, ab inito methods could be employed.
Another factor is the refinement of fitting schemes to estimate the parameters of Γ b . More efficient techniques could be utilized to compute the linewidths of the VDSF from larger systems without the need of explicitly diagonalizing large matrices [57].
Future work will address these issues to improve the accuracy of our method and enhance the quantitative agreement with experimental results [58].

Conclusions
We have developed an effective model for the lowfrequency excitations in amorphous solids that allows one to accurately compute the bulk limit of the thermal conductivity from glass models of moderate size. Our method, which stands on a combination of the QHGK approximation with various ideas that have been floating around in the literature for some time now, naturally accounts for the interplay of anharmonicity and disorder in determining the transport properties of glasses. The resulting protocol gets around the need for mock crystalline models for the glass, which introduce spurious modes whose group velocities bear little physical meaning. This fact is exemplified by the extreme case of a 2×2×2 q-mesh, where group velocities vanish identically as all the points of the mesh lie on the surface of the BZ and differ by a reciprocal-lattice vector from their negatives. As a result, the thermal conductivity computed on a 2 × 2 × 2 q-mesh is practically indistinguishable from the one computed at the zone center. Moreover, sampling the full BZ of a crystalline model for a glass only affects the first, intraband RTA-like, term in Eq. (10), while the second is already at convergence at the zone center for rather small glass models. The RTA-like term grows with the size of the q-mesh until it converges to some value, whose magnitude may be very large and devoid of any significance.
The only requirement of our extrapolation technique is a finite sample whose size can range from hundreds up to a few thousand atoms, depending on the characteristics of the material. We have tested our model on three paradigmatic glassy materials, aSiC, aSiO 2 , and aSi, which display different convergence properties to the bulk limit. When the effect of disorder is prominent (e.g., in amorphous silica), size effects are small, and hydrodynamic extrapolation contributes little to the converged value of κ. On the other hand, when the glass is less disordered (e.g., amorphous silicon), tens of thousands of atoms are required to obtain a converged value of the bulk thermal conductivity, while the hydrodynamic extrapolation gives satisfactory results with an order of magnitude fewer atoms.

Acoustic-sound-waves basis
For a system of N atoms, |Q, b⟩ is a 3N -dimensional vector whose projection on the Ith atomic site in the α direction is: where ε b (Q), with b = L, T 1 , T 2 , are three orthonormal polarization unit vectors. The scalar product between two plane-wave states is: where the last sum is proportional to the Fourier Transform of the atomic number density, ρ(r) = I δ(r− R I ), i.e.ρ Assuming that the material is homogeneous at length scales larger than a certain wavelength, λ, implies that ρ(k < 2π/λ) tends to a Dirac-delta function; consequently: Therefore, the subset of vectors |Q, b⟩ with |Q| < 2π/λ is effectively orthonormal. Regarding the completeness problem, we are only interested in describing the propagons, not all the normal modes, which on the other hand would require a basis of 3N vectors. It can be argued that this set is sufficient for this purpose [40], as it can also be qualitatively understood from the plot of the DSF, where low-frequency modes are decomposed in small-Q plane waves only.
We are now in a position to write the propagon contribution to thermal conductivity in the plane-wave basis. The first term in the right-hand side of Eq. 15 in the main text couples only pairs of propagons. This is not a property of the energy flux operator, which has nonzero components for all different pairs of normal modes, but it is a consequence of the narrow Lorentzian functions appearing in the QHGK theory, Eq. (7). Taking into account this separation in frequency thanks to the Heaviside step function, let us consider the energy flux operator where the P subscript means that only pairs of propagons are involved; the anti-resonant part of the flux, which gives a negligible contribution to thermal conductivity [2], is omitted. The generalized velocity matrix is defined as [2,18]: where α, β, γ are Cartesian coordinates and Φ Jα is the Dynamical Matrix element between atoms I and J, whose masses and positions are respectively M I , M J and R I , R J , and V is the potential energy. To unclutter the notation, from now on we omit the sums over Cartesian-coordinate labels, which is understood when indices are repeated. In order to get to J bb ′ Q , we first expand the normal modes in plane waves and we observe that for propagons: The important point here is that we are able to factor ω ν out of the sum due to the almost-linear dispersion of the DSF. This translates into an orthogonality relation of normal modes whose frequency is far from the dispersion line, i.e., |⟨Q, b|ν⟩| ≈ 0 if |ω n − c b Q| ≪ Γ b (Q). Plugging the last equations into the energy flux operator of the propagons yields: Under the assumption that the material is practically homogeneous above the λ scale, the dynamical matrix can only depend on R I − R J ; therefore, with a change of variables (R I , and then: where From the DSF, we know that |Q, b⟩ is practically a linear combination of almost-degenerate normal modes, that are eigenvectors of the dynamical matrix with eigenvalue c 2 b Q 2 ; therefore, for b = b ′ , we have that: Thus, the energy flux becomes The mixed-polarization terms do not contribute to the value of thermal conductivity, since the two polarizations are well separated in frequency. As such, these terms are not included in our calculations.

Melt-and-quench generation of glass samples
Amorphous samples are obtained through a melt-andquench procedure starting from a crystalline conventional cubic cell replicated ℓ times along each Cartesian direction. The molecular dynamics simulations are carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS [53]). The meltand-quench procedures are performed according to the following recipes: aSiO 2 . The interatomic forces are described with the Vashishta force-field [59]. Simulations are carried out with a timestep of 1 fs. Starting from the β-cristobalite cubic conventional unit cell with a mass density of 2.20, g cm −3 , the crystal is melted at 7000 K. The molten sample is quenched from 7000 K to 500 K in 10 ns [10,60]. The sample is then thermalized at 500 K for 400 ps, and finally equilibrated for 100 ps in the N V E ensemble. For each size, we randomized the seed of the thermostat to obtain ten different samples. The average mass density of the amorphous samples is 2.43 g cm −3 , with a standard deviation across sizes of 0.01 g cm −3 .
aSiC. The interatomic forces are described with the Vashishta force-field [61,62]. Simulations are carried out with a timestep of 1 fs. For each size, the starting configuration is a crystalline cubic zincblend structure with a mass density of 3.22 g cm −3 . Following Ref. 61, the system is gradually heated from 300 K to 4000 K at constant null pressure. The solid/liquid transition is characterized by a sharp increase in the volume between 3000 and 3500 K [61]. Ten different liquid configurations are extracted every 2 ps from an N V E trajectory. The molten configurations are then quenched to 1000 K in 300 ps, and thermalized for 80 ps at constant pressure. The system is further cooled down to 500 K with the same procedure. The average mass density of the amorphous samples is 2.98 g cm −3 , with a standard deviation across sizes of 0.01 g cm −3 .
aSi. The interatomic forces are described with the Tersoff force-field [63]. Simulations are carried out with a timestep of 0.5 fs. For each size, the starting configuration is a crystalline diamond conventional unit cell replicated ℓ times along each Cartesian direction. Following Ref. 64, the crystal is melted at 6000 K and brought to 3000 K in 2 ns at fixed zero pressure. Then the system is equilibrated at 3000 K for another 2 ns. The molten sample is successively quenched from 3000 K to 2000 K in 10 ns, and finally annealed from 2000 K to 300 K at fixed volume in another 10 ns. Each glassy sample is equilibrated at 300 K for 10 ns. For each size, ten different samples are prepared according to this recipe beginning the quenching procedure from liquid configurations obtained initializing the atomic velocities with different random seeds. The average mass density of the amorphous samples is 2.275 g cm −3 , with a standard deviation across sizes of 0.003 g cm −3 .

Computation of the dynamical structure factors
The dynamical structure factors (DSFs) are computed according to Eq. (13) in the main text. For low frequencies, the anharmonic linewidths become smaller than the finite-size frequency spacing. When this happens, the Lorentzian functions in Eq. (13) are so narrow that few to no data points fall within its width. Thus, the contribution to the DSF coming from low-frequency modes is undersampled and affected by numerical noise. Under the assumption that the DSF is a Lorentzian function, one can overcome this issue by exploiting the closure of Lorentzian functions under convolutions. In fact, introducing the shorthand notation A ν (Q) ≡ |⟨ν|Q, b⟩| 2 , and L(ω, η) ≡ 1 π η ω 2 + η 2 , a modified version of the DSF can be introduced as: S η b (ω, Q) = ν |⟨ν|Q, b⟩| 2 1 π γ ν + η (γ ν + η) 2 + (ω − ω ν ) 2 = ν A n (Q)L(ω − ω ν , γ ν + η) (38) Then, by the closure of Lorentzians under convolution: Since, for propagons, we assume that S b (ω, Q) is a Lorentzian function centered on ω b (Q) = c b Q and whose width is Γ b (Q), we get Due to this property, one can compute Γ b (Q) subtracting η from the width of S η b (ω, Q) numerically calculated from Eq. (38) with a large-enough η, since S η b (ω, Q) is much less affected by undersampling than S b (ω, Q).

Interpolation scheme for the anharmonic linewidths
The computation of third-order anharmonic linewidths in glasses generally constitutes the bottleneck of QHGK calculations [18], as it scales as N 3 atoms . Thus, this computation by itself would severely limit our ability to study size effects. It is thus customary to adopt some interpolation scheme for obtaining the linewidths of a large model starting from a smaller one [2,7,18]. First, in order to smoothen the data, we apply a Gaussian filter to the computed anharmonic linewidths at each temperature, γ ν (T ): Then, we spline-interpolate the smoothed function imposing a quadratic behavior at vanishing frequencies, which is compatible with the hydrodynamics of solids [24]. We average γ(ω, T ) over disorder by computing the mean of the interpolated linewidths over different same-size samples of each material. The spline functions are finally evaluated on the frequencies of larger samples in order to obtain their approximated anharmonic linewidths.

DATA AVAILABILITY
The data that support the plots and relevant results within this paper are available on the Materials Cloud platform [65]. See DOI:10.24435/materialscloud:k2-0n.

CODE AVAILABILITY
The codes that support the relevant results within this paper are available to the respective developers. Analysis scripts are available on GitHub [66] and on the Materials Cloud platform [65]. See DOI:10.24435/materialscloud:k2-0n.