## Thermal conductivity in modified oxide glasses is governed by modal phase changes

Philip Rasmussen 1 and Søren S. Sørensen ∗ , 1 Department of Chemistry and Bioscience, Aalborg University, Fredrik Bajers Vej 7H, 9220 Aalborg,

Denmark

(*Electronic mail: soe@bio.aau.dk)

(Dated: 5 August 2024)

The thermal conductivity of glasses is well-known to be significantly harder to theoretically describe compared to crystalline materials. Because of this fact, the fundamental understanding of thermal conductivity in glasses remain extremely poor when moving beyond the case of simple glasses, e.g., glassy SiO2, and into so-called 'modified' oxide glasses, that is, glasses where other oxides (e.g. alkali oxides) have been added to break up the network and alter e.g. elastic and thermal properties. This lack of knowledge is apparent despite how modified glasses comprise the far majority of known glasses. In the present work we study an archetypical series of sodium silicate ( x Na2O- 100-( x ) SiO2) glasses. Analyses of modal contributions reveal how increasing Na2O content induces increasing vibrational localization with a change of vibrations to be less ordered, and a related general decrease in modal contributions to thermal conductivity. We find the vibrational phases (acoustic vs. optical) of sodium vibrations to be relatively disordered compared to the network-forming silicon and oxygen species, explaining how increasing Na2O content decreases thermal conductivity. Our work sheds new light on the fundamentals of glassy heat transfer as well as the interplay between thermal conduction and modal characteristics in glasses.

## I. INTRODUCTION

Thermal conductivity ( κ ) is a fundamental materials property with relation to both specific applications in e.g. insulating materials and thermoelectrics, as well as fundamental linkages to the underlying vibrational character of materials 1-4 . However, for the property of thermal conductivity, big differences between solid materials families exist. The perhaps biggest difference being that between crystalline and amorphous/glassy materials, where crystals generally feature larger κ than their amorphous counterparts 2,3,5 . While glasses (i.e., amorphous materials exhibiting a glass transition 6 ) is a somewhat diverse materials family, the most industrially relevant glasses are the oxides with practical applications in bulk applications such as construction and screens, yet they also see use in more specialized applications as e.g. electrical insulation materials in microchip manufacturing 7,8 due to their electrically-insulating properties. In detail, oxide glasses are chemically distinct by being composed of network formers (e.g., SiO 2 , GeO2, B2O3, P2O5) 9 , often in combination with so-called network modifiers (alkali, alkaline earth, and transition metal oxides) 8,10 . Structurally, the network formers comprise the back-bone (e.g. Si-O-Si) of the glass, while the introduction of the network modifier breaks up the backbone structure by introducing so-called non-bridging oxygens (e.g., Si-O -Na + ) where the metal-ion charge balances the created negatively charged non-bridging oxygens 8,10 .

cal rules for tweaking κ in oxide glasses are quite vast 11-13 , the underlying fundamental understanding of the heat propagation mechanism is much less studied, with the majority of studies focusing on pure SiO2 14,15 and, to our knowledge, only very few studies on modified oxides exist 16,17 . This lack of fundamental understanding is largely due to how the the lack of long-range structural order complicates the theoretical description of heat propagation which is otherwise applicable for crystals (and has been for decades) 18 as well as the lack of reliable interatomic potentials capable of describing the vibrational characteristicss of these systems. However, very recently, two theories of estimating thermal conductivity were unified to comprehend both amorphous and crystalline materials as well as structures in between these limits 1,4,19,20 based on lattice dynamics to the third order. This has allowed a significant broadening in the possible studies capable of being performed. Yet despite the improved theoretical description, studies are still limited and a general description beyond that of the simplest single-oxide case is lacking.

Generally, the range of thermal conductivity in both organic and inorganic glass materials is rather narrow (typically 0.1-2.0 W m -1 K 1 ) compared to that of crystals (0.1--10 4 W m -1 K 1 ). -For oxides specifically, both pure network formers as well as modified oxide glasses are experimentally rather well-studied, typically covering a κ -range of 0.4-1.4 W m -1 K 111 , -with pure network formers being in the range of 0.5-1.4 W m -1 K 111 . -However, while empiri-

To overcome this, in the present work we showcase how an existing classical force field is capable of reproducing the structural and vibrational features of several archetypical sodium silicate glasses (in the series of x Na2O- 100-( x ) SiO2 with x =[20-35]). Using the recent advances in modelling thermal conductivity, we furthermore show how the used force field provides a superior estimation of κ in the given systems in very good agreement with experimental findings. We supplement this estimation with analyses of mode localization and mode character to provide general knowledge on how addition of modifying oxides to pure network forming glasses affect the fundamental modal types and how these contribute to the thermal conductivity in modified oxide glass systems. Ultimately, our results provide new knowledge on the fundamentals of the connection between heat transfer and structure in archetypical modified glass systems and may aid in future development of glasses of ultra-high or -low thermal conduc-

tivity.

## II. METHODS

## A. Glass preparation

The archetypal sodium oxide glass system x Na2O- 100-( x ) SiO2 with x = {0, 10, 15, 20, 25, 30, 33.3, 35, 40} has been investigated by employing the well-established potential of Teter with the potential form,

$$U ( r _ { i j } ) = \frac { q _ { i } q _ { j } } { 4 \pi r _ { i j } \epsilon _ { 0 } } + A _ { i j } \exp \left ( - \frac { r _ { i j } } { \rho _ { i j } } \right ) - \frac { C _ { i j } } { r _ { i j } ^ { 6 } } + \frac { D _ { i j } } { r _ { i j } ^ { 2 4 } }, \quad ( 1 )$$

where U is the interatomic potential energy, q is atomic charge, r is the interatomic separation, ε 0 is the permittivity of vacuum, while i and j denote atomic species. Here, the first term is the Coulombic term which accounts for the charge-charge interactions between atoms. The second and third term are the short range Buckingham potential with parameters Aij , ρ i j and Cij used as fitted by Teter 21,22 . The last term is a short-repulsive term added to prevent the "Buckingham catastrophe" at very short distances and high temperatures 23 . The correction is illustrated in Figure S1 in the Supplementary Material showing the correction provided to remove the infinitely deep potential well at short separating and hence prevent unphysical forces and ultimately crashing simulations. Potential parameters used in simulations are provided in Table I. The short-range interactions were evaluated with a cutoff of 8 Å while the long-range Coulombic interactions were evaluated in real space below 12 Å and using the PPPM algorithm using an accuracy of 10 -5 above 12 Å. Glasses were prepared using the melt-quench procedure using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package 24 with a fixed timestep of 1.0 fs

Box sizes of 3000 atoms were employed in simulations to study structural features, elastic properties and various vibrational characteristics to limit errors from finite size effects. The number of Si, O and Na atoms in each glass were fixed according the corresponding molar composition of the glasses investigated, and box sizes were estimated according to experimental glass densities (see Figure S2 in the Supplementary Material) adapted from various studies 25-27 . Detailed simulation cell information is provided in Table S1 in the Supplementary Material.

TABLE I. Interatomic potential parameters used in simulation of sodium silicate glasses.

| Interaction           |   A (eV) |    ρ (Å) |   C (eV Å 6 |   ) D (eV Å 24 ) |
|-----------------------|----------|----------|-------------|------------------|
| Na + 0 . 6 -O - 1 . 2 |  4383.76 | 0.243838 |      30.7   |                1 |
| O - 1 . 2 -O - 1 . 2  |  2029.22 | 0.343645 |     192.58  |              113 |
| Si + 2 . 4 -O - 1 . 2 | 13702.9  | 0.193817 |      54.681 |               29 |

Initially, 3000 atoms were placed randomly in a simulation box separated by a minimum of 1.6 Å to avoid unphysical dynamics at the initiating MD steps, followed by an energy mini- mization with an energy and force tolerance of 10 -10 . Glasses were then equilibrated at 5000 K for 2 ns to ensure complete memory loss of the initial atomic configuration. A subsequent linear cooling was performed to 300 K with a cooling rate of 10 K ps -1 . Finally, glasses were relaxed for 100 ps where atomic trajectories were averaged and analyzed over 100 ps. All simulations were performed in the NVT ensemble.

Short range structural features such as coordination numbers, angular distribution functions and bond lengths were estimated from the final part of the quenching procedure employing the partial radial distribution functions computed with the formula,

$$g _ { i j } ( r ) = \frac { \langle n _ { i j } ( r ) \rangle d r } { 4 \pi r ^ { 2 } d r \rho _ { 0 } C _ { j } },$$

where ρ 0 is the average atomic number density and nij is the average number of particles between r + dr . The angular distribution and coordination numbers were determined by investigating the first coordination shells of different atomic species where a cutoff of 2.2 Å and ∼ 3.3 Å were used for the Si-O and Na-O interactions, respectively.

## B. Structural characterization

In order to investigate the ability of the potential of Teter to predict realistic medium-range structure in sodium silicate glasses, we computed the total structure factor to compare it to experimental scattering data. To do so, first, the partial structure factors were computed,

$$\stackrel { \dots } { 2 \, \mathbb { A }. } } \\ \underset { \text{ising} } { \sin } \quad S _ { i j } ( Q ) = 1 + \rho _ { 0 } \int _ { 0 } ^ { R } 4 \pi r ^ { 2 } ( g _ { i j } - 1 ) \frac { \sin ( Q r ) } { Q r } \frac { \sin ( \pi r / R ) } { \pi r / R } d r, \ ( 3 )$$

where Q is the wave vector, ρ 0 is the average atomic density, R is the real space maximum value of integration, gij is the partial radial distribution function and sin ( π r / R ) / ( π r / R ) is a Lorch-type function used to reduce effects from having finite box sizes. The neutron total structure factor was computed from the partial structure factors,

$$S _ { N } ( Q ) = \left ( \sum _ { i = 1 } ^ { n } c _ { i } b _ { i } \right ) ^ { - 2 } \sum _ { i, j = 1 } ^ { n } c _ { i } c _ { j } b _ { i } b _ { j } S _ { i j } ( Q ), \quad ( 4 )$$

where ci is the fraction of atom species i and bi is neutron scattering length of atomic species i . The neutron scattering lengths are 3 63 . · 10 -15 m, 5 803 . · 10 -15 m and 4 1491 . · 10 -15 mfor sodium, oxygen, and silicon, respectively 28 . Simulated structure factors were evaluated by comparing with experimental structure factors from experimental neutron scattering experiments by employing the R χ -factor as proposed by Wright 29 ,

$$R _ { \chi } = 1 0 0 \% \, \sqrt { \frac { \sum _ { r } [ S _ { N } ( Q ) ^ { \exp } - S _ { N } ( Q ) ^ { \sin } ] ^ { 2 } } { \sum _ { r } [ S _ { N } ( Q ) ^ { \exp } ] ^ { 2 } } }, \quad ( 5 )$$

providing a quantitative measure of the agreement between simulated and experimental structure. Generally, a simulated structure with R χ &lt;10% is considered to have a good agreement with experiment.

## C. Mechanical properties

To further validate the ability of the Teter potential 21,22 to predict structure and properties, the elastic properties of simulated sodium silicate glasses were estimated by displacing the simulation boxes with a finite difference in 6 perturbations, namely in in the xx , yy , zz , xy , xz and yz directions, while measuring the 6 components of the pressure tensor, i.e. the directional pressure components. The simulation boxes were displaced in steps of 0.1 % of the current box length for the tensile strain and in steps of 0.1 Å for the shear strain. Performing linear regression in the elastic region yielded the elastic constants C 11 and C 44 , which are used to compute the elastic constant C 12 , bulk ( B ), Young's ( E ) and shear ( G ) moduli and Poisson's ratio ( v ),

$$C _ { 1 2 } = C _ { 1 1 } - 2 C _ { 4 4 }, \quad \ \ ( 6 ) \, \stackrel { \dots } { \ g l a } \\ \ u s i$$

$$B = \frac { C _ { 1 1 } + 2 C _ { 1 2 } } { 3 }, \quad \quad \quad ( 7 ) \quad \text{pe} \\ \text{th} \text{ }$$

$$E = \frac { ( C _ { 1 1 } - C _ { 1 2 } ) ( C _ { 1 1 } + 2 C _ { 1 2 } ) } { C _ { 1 1 } + C _ { 1 2 } }, \quad \quad ( 8 ) \quad \text{$\mathfrak E$.}$$

$$G = \frac { C _ { 1 1 } - C _ { 1 2 } } { 2 } = C _ { 4 4 }, \quad \ \ ( 9 ) \quad \ w a \\ \text{for}$$

$$v = \frac { C _ { 1 2 } } { C _ { 1 1 } + C _ { 1 2 } }. \quad \quad \ \ ( 1 0 )$$

## D. Thermal conductivity

To investigate the ability of the Teter potential to predict vibrational and thermal properties in sodium silicate glasses, we compute the thermal conductivity of simulated glasses by employing the recently developed Quasi Harmonic GreenKubo (QHGK) method 19,20 as implemented in the κ ALD o software 30 . The method is a unification of the phonon gas model 1 and the theory of heat transport in glasses provided by Allen and Feldman 4 , allowing for a generalized description of both crystalline and amorphous materials 20 , where κ is given as,

$$\kappa _ { \alpha \beta } = \frac { 1 } { V } \sum _ { n m } C _ { n m } v _ { n m } ^ { \alpha } v _ { n m } ^ { \beta } \tau _ { n m } ^ { \circ }, \quad \quad ( 1 1 ) \quad \text{max} \quad ( 1 1 ) \quad \text{max}$$

where V is the simulation cell size, n and m are mode indexes, α and β are cartesian components, Cnm is the generalized heat capacity, v α nm is the generalized phonon group velocities and τ ◦ nm is the generalized phonon relaxation times. Then, Eq. 11 conveniently provides the per mode thermal conductivity of simulated glasses. In order to estimate thermal conductivity with the QHGK method, the second and third-order IFC's were computed with the LAMMPS package (using the built-in finite displacement method using displacements of 10 -6 Å), and subsequently imported into κ ALD o along with a meltquenched energy minimized glass structure file. Computation of the third-order IFC's takes a high computational effort since the number of matrix elements scales like (3 N ) 3 where N is the number of atoms. Hence, box sizes of 600 atoms were employed for the QHGK calculations to significantly reduce the computational expenses. Simulation cell information of 600 atom systems can be found in Table S2 in the Supplementary Material. Moreover, evaluating the long-range Coulombic forces in reciprocal space with the PPPM algorithm was found to yield unphysically low thermal conductivity values and a significant increase in computational cost. Instead, similarly to a previous study 17 , the Coulombic damped shift force model (DSF) was employed for simulations of the 600 atom glasses and computation of the second and third-order IFC's, using a damping factor of 0.3 and a cutoff of 11 Å, which conveniently reduces the third-order IFC file sizes and computational demand. Thermal conductivity estimations were then performed at 300 K, while also taking quantum effects (i.e., the Debye scaling of the heat capacity) into account.

## E. Vibrational Characteristics

Further investigation of vibrational and thermal properties was conducted by evaluating the second order interatomic force constants (IFC). The second order IFC's ( Φ ′′ i α β j ) were then used to obtain the dynamical matrix ( D ) which is the second IFCs rescaled by the atomic masses, that is,

$$\mathbf D = \frac { \Phi _ { i \alpha j \beta } ^ { ^ { \prime \prime } } } { \sqrt { m _ { i } m _ { j } } },$$

where mi is the atomic mass of atom i . Then, from the dynamical matrix, the mode eigenfrequencies ( Ω ) and eigenvectors ( e ) were obtained by solving the following eigenvalue equation,

$$\mathbf e \cdot \Omega = D \cdot \mathbf e,$$

where e contains 3 N column vectors ( e i ), where N is the number of atoms, which provides information about the relative motion of all atoms in each eigenmode, and Ω is a square matrix with diagonal elements of the squared angular modal eigenfrequencies,

$$\Omega = \begin{bmatrix} \omega _ { 1 } ^ { 2 } & & & & & & \text{ulat} \\ & \omega _ { 2 } ^ { 2 } & & & \\ & & \omega _ { 3 } ^ { 2 } & & \\ & & & \cdot & \\ & & & & \omega _ { 3 N } ^ { 2 } \end{bmatrix} ^ { \quad \begin{pmatrix} \ u l a t \\ ( A P \end{pmatrix} } ( 1 4 ) & & & \\ \dots \dots \vdots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots \vdots \dots$$

which are readily converted into ordinary frequencies by ω 2 / 2 π . Wethen bin the 3 N eigenfrequencies into a histogram to obtain the vibrational density of states ( g ( ω )),

$$g ( \omega ) = \sum _ { n } ^ { 3 N } ( \omega _ { n } - \omega ) \delta, \text{ \quad \ \ } ( 1 5 ) \text{ \ \ } \text{co}$$

where δ is the Dirac delta. To further investigate the modal vibrational characteristics, we exploit the modal eigenvector normalization,

$$\sum _ { i } ^ { N } \vec { \mathbf e } _ { i } ( n ) \cdot \vec { \mathbf e } _ { i } ( n ) = \sum _ { i } ^ { N _ { N a } } \vec { \mathbf e } _ { i } ( n ) \cdot \vec { \mathbf e } _ { i } ( n ) + \int _ { i } ^ { N _ { O } } \vec { \mathbf e } _ { i } ( n ) \cdot \vec { \mathbf e } _ { i } ( n ) = 1, \int _ { i } ^ { N _ { O } } \vec { \mathbf e } _ { i } ( n ) \cdot \vec { \mathbf e } _ { i } ( n ) = 1, \int _ { i } ^ { N _ { O } } \vec { \mathbf e } _ { i } ( n ) \cdot \vec { \mathbf e } _ { i } ( n ) = 1,$$

where i denotes the atomic index number, n is the modal index number and N Na/O/Si is the number of sodium, oxygen, and silicon atoms, respectively. The inner product of each eigenvector essentially yield weighting factors for each atom and atomic species, providing a measure of the vibrational participation of atomic species in the collective modal motion. This is utilized to obtain the partial g ( ω ) 31

$$g _ { \alpha } ( \omega ) = g ( \omega ) \sum _ { i \in \alpha } | \vec { \mathbf e } _ { i } ( n ) | ^ { 2 }, \quad \quad ( 1 7 ) \quad \begin{pmatrix} \text{of} \, \mathfrak { i } \\ \text{its} \, \mathfrak { \ }pha } \end{pmatrix}$$

where α denotes atomic species with α = (Si, Na, O, BO, NBO, Si Q 1 , Si Q 2 , Si Q 3 , Si Q 4 ), where BO and NBO is bridging and non-bridging oxygen, respectively, and Si Q n is silicon atoms in Q n species of different degrees of polymerization ( n related to the number of bridging-oxygens attached to the Si atom).

We further utilize the eigenvector normalization from Eq. 16 to compute the per mode participation ratio ( PR ) to estimate the per type participation in each mode 32 ,

$$P R ( n ) = \left ( N \sum _ { i } ^ { N } | \vec { \mathrm e } _ { i } ( n ) | ^ { 4 } \right ) ^ { - 1 }, \quad \quad ( 1 8 ) \quad \mathbf A.$$

where all individual atomic participation ratio values for atoms i are summed over the entire n -th mode. Practically, the limit of PR = 1 shows an equal contribution from all atoms to the vibrations while the limit of 1 / N indicates a mode of only 1 atom contributing to the modal vibration.

To gain further insight into the participation ratio in simulated glasses we employ the atomic participation ratio ( APR ) 33 ,

$$A P R _ { \alpha } ( n ) = \frac { | \vec { \mathbf e } _ { i } ( n ) | ^ { 2 } } { \left ( N \sum _ { i \in \alpha } ^ { N _ { i } } | \vec { \mathbf e } _ { i } ( n ) | ^ { 4 } \right ) ^ { \frac { 1 } { 2 } } }, \quad \ \ ( 1 9 ) \\. \quad \. \quad. \quad. \quad. \quad. \quad.$$

which is an atomic decomposition of the participation ratio, providing the participation ratio of individual atomic species ( APR α ) that satisfies PR n ( ) = ∑ i ∈ α ( APR α ( n )) 2 .

Next, to distinguish localized modes from delocalized modes we compute the inverse participation ratio, simply by computing the inverse of the PR,

$$I P R ( n ) = ( P R ( n ) ) ^ { - 1 }, \text{ \quad \ \ } ( 2 0 )$$

which provides a measure of the frequency range at which the mobility edge appears and localized modes are prevalent.

To further characterize the modal motion, the phase quotient ( PQ ), originally proposed by Bell and Hibbins-Butler 34 , was used to distinguish between acoustic ( PQ = 1) and optical ( PQ = -1) modes, which provides a measure of whether nearest neighbor atoms vibrate mostly in phase or out of phase in each mode. Specifically, PQ is given as,

$$P Q ( n ) = \frac { \sum _ { m } \vec { \mathbf e } _ { i } ( n ) \cdot \vec { \mathbf e } _ { j } ( n ) } { \sum _ { m } | \mathbf e _ { i } ( n ) \cdot \vec { \mathbf e } _ { j } ( n ) | }. \quad \quad ( 2 1 ) \\.. \dots \dots \dots \dots, \dots \dots, \dots...............$$

We note that the terms 'acoustic' and 'optical' are commonly used to define fully in-phase and out of phase modes from the dispersion relation of crystalline phonons, yet, the terminology still somewhat applies and conveniently provides a quantitative measure of the average degree of phase motion of modes in glasses. We further study the atomic motion and its relation to nearest neighbors by investigating the atomic phase quotient ( APQ ) of different atomic species,

$$\ A P Q ( n ) & = \frac { \sum _ { m } \mathbf e _ { i, \alpha } ^ { \dagger } ( n ) \cdot \vec { \mathbf e } _ { j } ( n ) } { \sum _ { m } | \mathbf e _ { i, \alpha } ^ { \dagger } ( n ) \cdot \vec { \mathbf e } _ { j } ( n ) | }, \quad \quad ( 2 2 ) \\ \cdot \quad. \quad \cdot \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

which in nature is similar to the PQ but describes the direction of eigenvectors of each nearest neighbor atom adjacent to specific atomic species.

## III. RESULTS AND DISCUSSION

## A. Structural and mechanical characterization

In this work we employ the force field of Teter 21,22 which has previously been shown to provide a good description of the structure and mechanics in a number of alkali silicates 22,35-40 . We first prepare x Na2O- 100-( x ) SiO2 glasses with x = {0, 10, 15, 20, 25, 30, 33.3, 35, 40} to make an initial assesment of properties in a very broad compositional range.

First, to probe the structural agreement with experiments we present the calculated and measured neutron structure factor (experimental S Q ( ) from Ref. 41 ) in Figure 1a for a 30Na2O70SiO2 glass as well as a similar comparison for pure SiO2 in Figure S3 in the Supplementary Material.

FIG. 1. (a) Neutron-weighted structure factor of the studied glass with a composition of 30Na2O-70SiO2 from molecular dynamics (MD) simulations (green) as well as from experiments (black) showing a very good agreement and (b) Young's (E), Bulk (B) and Shear (G) moduli of simulated glasses (filled circles) compared to experimental values (open circles). The experimental structure factor data is from Ref. 41 and moduli data are from Ref. 42

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-2b041207-1683-4241-b94f-2191c798f1de.jpg)

Notably, by comparing the experimental and simulated S Q ( ) using the R χ factor introduced by Wright 29 , both 30Na2O-70SiO2 and SiO2 glasses see values well below 10%, indicating how the simulations provide a very good structural depiction of the experimental glass in good agreement with previous reports of the present pair potential 21,22 . We supplement the structure factor comparisons with calculated Q n -distributions (where n indicates the fraction of bridging oxygens per Si tetrahedra), bond-angle distributions (Si-Bridging oxygen-Si, O-Si-O as well as various Na-O related bond angles), total and partial structure factors as well as estimations of the fraction of bridging and non-bridging oxygen species for all studied compositions. These data are presented in Figures S4-11 in the Supplementary Material. While not all parameters are accessible experimentally, for those available we find good agreement between the simulated glasses and experiments. For example, we identify O-Si-O bond centered around angles of ∼ 109 ◦ (Figure S6) as expected by a tetrahedral configuration, overall good qualitative agreement with Q n distributions (Figure S11a), and an expected decrease in bridging oxygen with increasing Na2O content (Figure S11b).

Next, we probe the mechanical properties of the system by performing displacements in the tensile and shear directions to record the elastic constants of the system and ultimately to estimate the elastic moduli 43 (see Methods for details). The results are shown in Figure 1b. From this analysis we find decreasing Young's ( E ), bulk ( B ) and shear ( G ) moduli with increasing Na2O content. This is in good agreement with experiments for E and G , but disagrees with the somewhat constant experimental B (open green symbols in Figure 1b as obtained from Ref. 42 ) in the studied modifier concentration range. This disagreement is likely also the cause of the deviation between simulated and experimental Poisson's ratio although the general increase with increasing Na2O trend is reproduced (see Figure S12 in the Supplementary Material). In summary, we generally find good agreement between experimental and simulated mechanical properties, especially in the 20-35 mol% Na2O range, highlighting how the force field performs best in this concentration range. Given how the present study will focus on thermal transport, a property somewhat related to mechanical properties, the remaining part of this work will focus on x Na2O- 100-( x ) SiO2 glasses with Na2O contents of x = {20, 25, 30, 33.3, 35}.

## B. Vibrational features

While a number of properties have previously been tested using the Teter potential 21,22 , the force field has mainly been employed due to its superior modelling of structural and mechanical characteristics of modified alkali silicate glasses 25,35-37,44,45 . A few studies have probed the total and partial vibrational density of states (VDOS) 31,45 however with emphasis on the density/pressure response of glasses at or near the Na2O-2SiO2 composition. As such, to our knowledge the potential have never been used to test the thermal conductivity nor any detailed vibrational characteristics of glasses. In the present work we perform a detailed vibrational analysis of the compositions of interest as well as provide the first use of the potential to simulate the thermal conductivity of these glass types. First, we compute the total vibrational density of states ( g ( ω ) , see Methods) for a range of compositions as presented in Figure 2.

Comparing the simulated data of the glasses produced using the Teter potential to ab initio simulations (from Ref. 46 , See Figures S13-14) we find good agreement. For more, our calculations agree well with previous estimations of the VDOS from the same potential 31 . This generally suggests that our simulations are capable of decently describing the vibrational character of the studied glasses and thus supports the usage of

FIG. 2. (Top) Vibrational density of states ( g ( ω ) ) of simulated sodium silicate glasses with increasing sodium content from 20 to 35 mol% Na2O as well as the g ( ω ) from ab initio simulations as performed in Ref. 46 , and (bottom) illustrations of eigenvectors of modes in three different frequency ranges of g ( ω ) , i.e. (a) a low-frequency mode at ∼ 7 THz, (b) a mid-frequency mode at ∼ 20 THz and (c) a high-frequency mode at ∼ 32 THz, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-b6f8bed4-d579-4ec9-aa23-c21c04e40d04.jpg)

the potential to predict other vibrationally-related properties, e.g. thermal conductivity. In addition to the total VDOS we also calculate the partial vibrational density of states for O, Si, and Na (including subdivisions based on Q n distribution and oxygen coordination, see Figure S15 in the Supplementary Material as well as an example of the decomposed VDOS in Figure S16). From this analysis we find that the shift of the low-frequency band (0-10 THz) in Figure 2 to higher frequencies and intensity is caused by sodium contributions while the shifting in the 10-25 THz range are found to be caused by oxygen contributions upon increasing Na2O content.

We further show the eigenvectors (i.e., the solutions to Eq. 13) of three modes in Figure 2 to graphically illustrate the vast changes of modal characteristics in different frequency ranges. For example, in Figure 2a a low frequency mode at ∼ 7 THz is illustrated with a modal phase quotient (PQ) value (PQ=0.77), i.e. a mode where eigenvector displacements are mostly in-phase (vibrating in the same direction). Interestingly, this agrees with previous findings that modes in phase are common in the low-frequency region in both glasses and crystalline solids 47 , which could suggest that modes with high a PQ resemble that of a classical phonon. Observing the midfrequency range ( ∼ 10-30 THz) where modes tend to shift towards diffusive heat transfer processes (often denoted as diffusons), we illustrate a mode at 20 THz in Figure 2b. Com- pared to Figure 2a we find that displacements seem directionally random yet that vibrations are still extended throughout the solid in the majority of atoms. Finally, we inspect a highfrequency mode (32 THz) in Figure 2c which shows highly localised eigendisplacements with vibrations centered around only few atoms. These modes are often referred to as locons which are commonly found at higher frequencies above the so-called 'mobility edge' (generally, modes of participation ratios lower than 0.15) 4 . To probe the effect of the depolymerization of the network with increasing Na2O content we provide a more elaborate illustration of eigenvectors in multiple compositions at a series of frequencies in Figure S17 in the Supplementary Material. Notably, it is found that eigenvectors in modes with a frequency of ∼ 1 THz become less periodic and more diffuse when increasing the sodium content - a feature often related to decreasing heat transfer 48,49 . Furthermore, eigenvectors in modes with a frequency of 15 THz in the mid-frequency range seems to become more localized in certain spatial volumes when increasing the sodium concentration, which is an indication of lowered participation ratio of eigenvectors when disrupting the silica network. Finally, in the high-frequency range at 35 THz, modes are generally very localized in all simulated glasses, yet upon increasing Na2O content, eigenvectors seem to gradually shift towards higher localization.

While the visual inspection of eigenvectors is useful to provide an overview, to quantitatively assess the qualitative observations of mode localization provided above, we compute the participation ratio as presented in Figure 3a. Here, we observe that mid-frequency (10-20 THz) modes possess the highest PR ( ∼ 0.4), suggesting that the atomic eigenvectors in modes similar to the one illustrated in Fig 2b are somewhat the same magnitude and therefore contribute evenly to the vibrational motion in diffusive processes. Upon increasing the sodium concentration we observe how the PR decreases in the mid-frequency range and shifted towards an increase in the low-frequency range. This observation coincides with the shift of the partial VDOS of sodium and bridging oxygen species (See Figure S15 in the Supplementary Material), where the contribution from sodium related vibrations increase in the low frequency range while the contribution from bridging oxygen species decrease in the 10-23 THz range. The partial VDOS of silicon show that the atomic contribution from Q 4 and Q 3 silicon species also decreases at ∼ 20 THz, likely due to the induced network depolymerization, suggesting that sodium atoms in the interstitial sites are responsible for localization changes of atomic vibrations in the silica backbone.

Additionally, from Figure 2 we find that high-frequency modes (&gt;27 THz) are slightly shifted towards lower frequencies upon increasing Na2O content. From Figure 3a and Figure S18 in the Supplementary Material we find that this shift is accompanied by a slight decrease in PR . Generally, modes &lt;27 THz have a PR &gt; . 0 15 and are thus non-localized while while modes &gt;27 THz are localized ( PR &lt; 0 15), albeit we . note that localization is rather a gradual than a discrete change.

To obtain further insight, we compute the atomic participation ratio ( APR α where α denotes an atomic type) in Figure S19 in the Supplementary Material, which provides a neat way of describing contributions of different contributions to PR , in a way, such that PR n ( ) = ∑ i ∈ α ( APR α ( n )) 2 .

Here we find that the increased participation ratio in the low-frequency range ( ∼ 0-8 THz) are in fact a consequence of decreased localization of sodium vibrations. These localization changes quite impactfully increase the localization of bridging-oxygen vibrations in the range ∼ 0-20 THz, suggesting that sodium has much higher impact on oxygen vibrations than silicon-related vibrations, which seems to be reasonable observation since sodium is coordinating only directly to oxygens in the present systems and therefore has only weak interactions with silicon.

Due to their amorphous nature, dispersion relations are not well-described in glasses 10 . Instead, to access the acoustical/optical character of modes we compute the so-called phase quotient ( PQ ) per mode as introduced by Bell and HibbinsButler 34 . We present the PQ in Fig 3b, illustrating three distinct regions; i.e. modes vibrating in-phase (i.e., eigenvectors of neighboring atoms have angles of &lt;90 , and thus ◦ PQ &gt; 0, see Methods for definition); modes where nearby atoms vibrate perpendicular to each nearest neighbor eigenvector ( PQ = 0); and finally modes out of phase (angle between eigenvectors &gt;90 ◦ and PQ &lt; 0). We find that the majority of low-frequency modes are in-phase vibrations,

FIG. 3. (a) Participation ratio (PR) of all studied sodium silicate glasses (right side of a) along with three modes illustrating the connection between modal character and PR; and (b) phase quotient (PQ) as well as an illustration of the relation between nearest neighbor eigenvectors at PQ ={-1, 0, 1} for all studied x Na2O- 100-( x ) SiO2 glasses with increasing sodium contents.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-9dcaac98-1451-4e51-bc86-7103826b18ea.jpg)

which is a common feature in both crystalline and amorphous solids 34,47,50 .

Here, the PQ of the simulated glasses transitions from inphase to out of phase at ∼ 13 THz and then decreases until about ∼ 25 THz. Modes above this range is found to increase in PQ , but appears to level off at PQ ∼ 0. A plot of the PQ of a 3000 atom system is shown in Figure S20 in the Supplemental Material, showing that high-frequency out-of-phase modes vary considerably in PQ

We already established that depolymerization of the silica network with Na2O addition promotes considerable changes in modal characteristics of low and mid-frequency modes ( ∼ 0-25 THz), e.g. disruption of somewhat periodic eigenvectors at 1 THz (based on qualitative observations, see Figure S17 in the Supplementary Material) and decreased localization of modes in the mid-frequency range ( ∼ 8-25 THz) as a consequence of increased modal density of sodium related vibrations and extended localization of modes in the low-frequency range ( ∼ 0-8 THz). When observing the PQ in the mentioned frequency ranges, we find that the PQ of modes converge towards more perpendicular (or random) vibrations, i.e. PQ → 0, besides a small range of modes with frequencies of ∼ 14-18 THz. The most notable changes are found in

the low-frequency range, where the PQ of modes on average decrease, which correlates with the observations found in the modal illustrations in Figure S17 where eigenvectors seem to become more diffusive. In addition to PQ , we compute the atomic phase quotient ( APQ , see Methods) of Si, O, and Na as a function of modal frequency (see Figure S21). We note that the total PQ is not the sum of the APQs . Rather, the APQ provides a description of the phase-characteristics on a atomtype basis which is an interesting way of distinguishing the local vibrational environments.

Generally, we observe that the majority of changes are related to silicon and oxygen vibrations, where silicon seems to become more in-phase with adjacent atoms while oxygen vibrations become less in phase with adjacent atoms when increasing the sodium content. This might be due to the increased order in the tetrahedral environments of silicon as observed in the O-Si-O angles provided in Figure S6 in the Supplementary Material. Another interesting observation is how the Na APQ is near 0, indicating a significant degree of randomicity in the vibrational directions compared to the nearest neighbors, likely due to the weak bonding environment around Na atoms - a feature which is becoming more pronounced with a higher amount of Na2O in the glass.

## C. Heat transfer

Despite some interest in the fundamental vibrational characteristics of a few of the sodium silicate glasses 31,46 , derived properties, e.g. thermal conduction, are practically unexplored in simulations. This lack of study is notable given the well-understood structure and mechanical properties of these glasses as well as their usefulness as an archetypical model system for modified alkali silicate glasses, as well as its similarity to many glasses used industrially (e.g., windows and screens) 8 . In this work we employ the recently developed method of estimating the thermal conductivity ( κ ) in a unified model 19,20 coined as the Quasi Harmonic Green Kubo (QHGK) 20 . This model allows the estimation of κ beyond the pure phonon, diffuson (i.e. non-localized modes which see short relaxation times) and locon (i.e. localized modes which see only vibration of few atoms and hence negligible contribution to κ ) regimes. Practically, the model requires the input of force constants to the third order. We computed these through a finite displacement method implemented in LAMMPS and performed the following QHGK calculations as implemented in kALDo 30 . We present the total thermal conductivity estimations in Figure 4a where we also compare with experimental values of thermal conductivity 51,52 .

Generally, we find very good agreement between simulated and experimental data in the observed range, which seems to be consistent with the well reproduced VDOS of simulated glasses (Figure 2). Notably, the glasses with molar compositions x = 20 and x = 33.3 deviate from experimental data with 2.6 and 1.3%, respectively, showcasing how the Teter potential is able to predict vibrational features and thermal conductivity nicely. This agreement is a significant finding given how MD simulations have been notoriously known

FIG. 4. Effect of sodium concentration on (a) total thermal conductivity (b), cumulative thermal conductivity and (c) modal thermal conductivity of simulated x Na2O- 100-( x ) SiO2 glasses. Experimental values are from Ref. 51,52 . Error bars for MD glasses are smaller than the symbols.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-51fa7169-4879-4601-ae5e-170086448866.jpg)

for overestimating κ of many glass systems 16,53 . Practically, the estimations show a linear decrease in κ from ∼ 1.2 to ∼ 0.9 W m -1 K 1 -when increasing the sodium concentration from 20 to 35 mol%, evidently showing that formation of non-bridging oxygen from network depolymerization results in linearly decreasing κ in the simulated glasses.

A clear strength of QHGK calculations is that it does not only provide the total κ but also an individual contribution of each mode to the total κ , a so-called 'modal' thermal conductivity. From this, we present the cumulative thermal conductivity in Figure 4b to investigate the contribution to the thermal conductivity ( κ cum) in different frequency ranges upon increasing the sodium concentration. Upon increasing Na2O contents we find that the contribution to κ from low-frequency modes (0-8 THz) remain somewhat unchanged, while the slope of the curves decrease in the mid-frequency range (8-

24 THz), i.e. the contribution to κ of modes in the midfrequency range decreases.

We furthermore present the modal contribution to the thermal conductivity ( κ mode) in Figure 4c to investigate the contribution to κ from each of the 3 N modes individually. From this plot, it is evident that lower frequency modes contribute the most to κ , and the modal contribution seems to decrease with increasing frequency. Notably, above 25 THz, modes exhibit very low contribution to κ in line with their very low PR (Figure 3a) 4 . For more, we find a significant decrease in κ mode values in the frequency range 0-8 THz when increasing the sodium concentration, while the number of modes increases as seen from the increased modal density in the VDOS (Figure 2), which explains the somewhat static κ cum in the lowfrequency range (Figure 4b). Notably, the increased modal density of the VDOS in the range 0-8 THz originates from increased vibrational contributions from sodium related vibrations as seen in the partial VDOS of sodium (See Figure S15 in the Supplementary Material), suggesting that sodium related vibrations are responsible for a decreased κ contribution in the lower frequency modes. In addition to values of κ , the modal diffusivities Di from the QHGK calculations are shown in Figure S22 in the Supplementary Material illustrating that the diffusivity of modes significantly decrease in the 0-25 THz frequency range upon increasing Na2O. However, we find that the diffusivity of modes feature a rather constant contribution across the frequency range. This difference in comparison to thermal conductivity (which see a decreasing contribution, see Figure 4c) is caused by the quantum effect of heat capacity with decreasing mode occupancy levels at higher frequencies.

## D. Correlating modal characteristics and heat transfer

Next, to provide additional insight to the heat transfer mechanisms of the studied glasses, the modal thermal conductivities provided in Figure 4c are coupled with the computed vibrational properties of PR and APQ (Figure 3 as well as Figures S19-21). In Figure 5a we compare the modal participation ratio of glasses with the modal thermal conductivity and study the effect of localization changes upon increasing the Na2O concentration.

Wefirst notice the small band in the PR range ∼ 0-0.2 where modes are contributing insignificantly to κ (i.e., κ mode &lt; 0 1 . · 10 -3 W m -1 K 1 ). -These are high frequency modes (&gt;27 THz) which see a high spatial localization, which is consistent with previous findings of other amorphous systems in the literature 4,54 . For modes found in the frequency range 8-27 THz, we find that κ mode increases in a linear fashion with increasing PR . Increasing the PR of modes could rationally lead to the assumption that the modal contribution to κ increases as a consequence of more extended energy propagation. However, for modes found in the frequency range 0-8 THz, the modal contribution to κ seems to be inversely proportional to PR , which is in disagreement with the trend found for modes in the range of 8-27 THz. These observations show that the correlation between the modal participation ratio and the per-mode contribution to κ of extended modes is not straightforward in the studied sodium silicate glasses, thus suggesting that PR may not be the ideal descriptor of the underlying heat transfer.

In Figures 5b-f, we coupled κ mode with the APR of different atomic species (Si, Na, O, as well as the decompositions into bridging and non-bridging oxygens). We first turn our attention to the APR of sodium atoms coupled with κ mode in Figure 5b. As previously assessed, the modal participation ratio increases in the low-frequency range when increasing the sodium concentration (Figure 3), which results from a decreased vibrational localization of sodium atoms. Furthermore, κ mode decreases as the localization of sodium related vibrations decreases, yet, above APR ∼ 0.22 the κ mode gradually increases in each glass with more than 25 mol% of Na2O, which could be correlated to a greater heat transfer contributed from more evenly distributed sodium vibrations. Since higher degree of polymerization in the silica backbone (and thus less Na2O) has shown to be associated with greater heat transfer, we now study the localization of silion and oxygen to investigate how disruption of the silica network impacts localization of the network forming atoms and ultimately κ .

In Figure 5c the κ mode coupled with APR of silicon is presented. Similar to the total PR (Figure 5a), a clear correlation between localization of silicon and κ mode is difficult to establish. Although, we observe that an increasing Na2O content increases the localization of silicon in almost the entire spectral range, which seems to correlate with a decrease in κ mode. In general, the silicon vibrations contribute very little to the modes of high κ mode as opposed to sodium, suggesting that

FIG. 5. Effect of sodium concentration on (a) total participation ratio and atomic participation ratio of (b) sodium, (c) silicon, (d) oxygen, (e) bridging-oxygens and (f) non-bridging oxygens coupled with modal thermal conductivities estimated with the QHGK method.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-fffbea53-91ab-4164-822d-801e521b454e.jpg)

the localization changes of silicon may have less impact on overall thermal conduction. Lastly, we showcase the correlation between oxygen APR and κ mode in Figure 5c-f. We find that oxygen in general show similar trends to that of the total PR , especially in glasses with lower sodium contents. In fact, when observing the APR of bridging oxygens presented in Figure 5d-f, we find that oxygen in the silica backbone possess higher atomic participation ratio than NBOs and silicon atoms as seen in Figure 5c,e,f , suggesting BOs to be the main contributor to the heat transfer in the mid-low frequency range ( ∼ 8-20 THz) in the studied glasses. As might be noted, the APR is dependent on the number of atoms in each sample. As such, the changing composition also inherently changes the APR . To make a more unbiased comparison we also computed an APR normalized by the number of the atom type in each composition and correlated this to κ mode. We present this plot in Figure S23 in the Supplementary Material showcasing how the range of the normalized APR is largely similar for the different atomic types across the varying compositions. This indicates a simple scaling of the APR with the composition change. In summary, while there seems to be some correlations between thermal conductivity and PR , the correlation is nontrivial.

Instead, we investigate the impact of in-phase and out-ofphase modes, through the use of the phase-quotient ( PQ ), and its relation to the thermal conductivity. We provide these correlations for the total system as well as for all atom types in Figure 6. We generally observe higher PQ (i.e., more in-phase vibrations) at lower frequencies (Figure 3) which is somewhat similar to many crystalline solids 47 . For more, we generally find that higher PQ is associated with higher modal κ (Figure 6). This may be correlated to the inherent periodicity of strongly acoustic phonons found in crystalline solids which are also known as strong contributors to thermal conductivity. This suggests that the in-phase modes observed in the sodium silicate glasses possess some degree of periodicity, and hence explain the relatively high contribution to κ . With increasing Na2O content we generally observe small changes to total PQ (Figures 3 and 6), however a slight decrease in the associated κ mode is found. For more, a broad band can be found in the range of PQ = -0 5 to . PQ = 0 with very low contribution to κ corresponding to the localized modes observed in Figure 3a. This indicates that localized modes with low contribution to the thermal conductivity inherently consist of out-of-phase vibrations of nearest neighbour atoms.

For more, from Figure 3 we find that the PQ becomes less spread in the low-frequency range, while also shifting slightly towards PQ = 0 (i.e., they loose their character as either acoustic or optical) in the whole spectral range when the sodium concentration is increased, indicating a sort of disordering of the vibrations. To study these features and their impact on the modal contribution to κ , we investigate the atomic phase quotient (APQ), which describes the angle between specific atomic types (O, Si, Na, bridging-, and non-bridging O) and the eigenvector of their nearest neighbour atoms (see illustration in Figure 3b). The APQs are plotted in Figure S21 in the Supplementary Material as a function of modal frequencies, showing some of the same trends as the total PQ shown in Figure 3. We now couple the PQ and APQs to κ mode (Figure 6) to investigate the impact of increasing the Na2O on nearest neighbor atomic vibrations and the thermal conductivity.

Similarly , when observing the sodium APQ provided in Figure 6b, we observe a rather distinct APQ trend. Although the correlation between APQ and κ mode is not easily defined, we observe that modes generally lie closer to PQ=0 than the silicon and oxygen APQs. This may be linked to how sodium atoms are more loosely bound in the glass structure, effec-

The silicon and oxygen APQs in Figure 6(c,d) show many similarities to the trends observed in the total PQ. For example, there seems to be a very large degree of similarity in plots in Figure 6 when changing the molar concentration of Na2O. Effectively, this means that modes containing similar oxygen and silicon vibrational characteristics contribute largely similarly to heat conduction, independent of the glass composition. The small decrease in κ mode for vibrations of otherwise similar PQ/APQ is likely caused by the disruptive nature of the Na2O addition. Interestingly, for oxygen these points seem to be gradually different for BOs and NBOs (Figure 6(e,f)), namely that the BOs seem to provide modes of stronger acoustic and optical character than NBOs. This points towards a slightly higher "disorder" in the vibrations of the NBOs, a less distinct correlation with κ mode and hence likely a smaller direct correlation with heat transfer.

FIG. 6. κ mode plotted as a function of (a) the total phase quotient (PQ) of all nearest neighbor atomic interactions in each mode, as well as the atomic phase quotient (APQ) providing the angle between eigenvectors of nearest neighbor atoms of (b) sodium, (c) silicon, (d) oxygen (e) bridging-oxygens, and (e) non-bridging oxygens, respectively, for each of the studied sodium silicate glasses. Dashed lines provide the transition point (PQ=0) from in-phase to out of phase modes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-d2e689ea-05f6-4331-9d48-4b7ca7d731e7.jpg)

tively causing the vibrational direction in each mode to be less coupled to the direction of vibration of its nearest neighbor (which is always an oxygen, most commonly in an octahedral state, see Table S1 in the Supplementary Material). Effectively, this is believed to disrupt effective heat transfer throughout the network despite its large contribution to lower frequency vibrations, effectively reducing the per-mode thermal conductivity contribution (Figures 3-6). Naturally, upon increasing contents of Na2O, this effect is enhanced, explaining the overall reduction in κ (Figure 4) in modified oxide glass systems.

## IV. CONCLUSIONS

In summary, using classical molecular dynamics simulations we studied a series of sodium silicate ( x Na2O- 100-( x ) SiO2) glasses with varying Na2O content ( x = 20 -35 mol%) using the recently developed quasi-harmonic Green Kubo method. We find the used classical potential to be capable of predicting the observed experimental thermal conductivity of the systems and the decreasing thermal conductivity with increasing Na2O content. By modal deconvolution of the contributions to the total thermal conductivity, we find the per-mode contribution to conductivity to generally decrease in a wide frequency span (0-25 THz) upon increasing Na2O content, however due to a shifting of modes from 8-25 THz to 0-8 THz, the overall contribution to thermal conductivity is solely decreasing in the 8-22 THz frequency range.

From a detailed vibrational analysis we identify a clear trend of increasing thermal conductivity with increasing inphase motion of the underlying vibrations, especially for the network forming atoms (Si and O). However, upon increasing Na2O concentrations we find atomic vibrations to become less 'ordered', effectively causing a disruption of heat transfer throughout the glass network. Ultimately, our results provide new insight to the connection between vibrational character and mode-resolved thermal conductivity in the family of modified oxide glasses and may aid in the future design of oxide glasses of tailored thermal properties.

## SUPPLEMENTARY MATERIAL

See Supplemental Material for tabulated values of structural, mechanical, and thermal properties as well as detailed plots of pair-potential, density, pair-distribution functions, bond-angle distributions, structure factors, vibrational density of states, eigenmodes, and various vibrational characteristics of the studied glass systems.

## AUTHOR DECLARATIONS

## Conflicts of interest

The authors declare that they have no competing interests.

## ACKNOWLEDGMENTS

The authors acknowledge computational resources provided by CLAAUDIA at Aalborg University. The present work was initiated as a student thesis project and we acknowledge the thesis opponent Nicholas Bailey (Roskilde University) for fruitful discussions.

## AUTHOR CONTRIBUTIONS

Philip Rasmussen : Conceptualization (supporting); Methodology (equal); Investigation (lead); Formal analysis (lead); Visualization (lead); Writing -original draft (lead); Writing - review and editing (supporting). Søren S. Sørensen : Conceptualization (lead); Methodology (equal); Formal analysis (supporting); Supervision (lead); Writing - original draft (supporting); Writing review and editing (lead).

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

- 1 R. Peierls, 'Zur kinetischen theorie der wärmeleitung in kristallen,' Annalen der Physik 395 , 1055-1101 (1929).
- 3 R. C. Zeller and R. O. Pohl, 'Thermal conductivity and specific heat of non-crystalline solids,' Phys. Rev. B 4 , 2029-2041 (1971).
- 2 C. Kittel, 'Interpretation of the thermal conductivity of glasses,' Physical Review 75 , 972-974 (1949).
- 4 P. B. Allen and J. L. Feldman, 'Thermal conductivity of glasses: Theory and application to amorphous si,' Physical Review Letters 62 , 645-648 (1989).
- 6 E. D. Zanotto and J. C. Mauro, 'The glassy state of matter: Its definition
- 5 S. S. Sørensen, M. B. Østergaard, M. Stepniewska, H. Johra, Y. Yue, and M. M. Smedskjaer, 'Metal-organic framework glasses possess higher thermal conductivity than their crystalline counterparts,' ACS Applied Materials and Interfaces 12 , 18893-18903 (2020).
- and ultimate fate,' Journal of Non-Crystalline Solids 471 , 490-495 (2017). 7 W.-X. Zhou, Y. Cheng, K.-Q. Chen, G. Xie, T. Wang, and G. Zhang, 'Thermal conductivity of amorphous materials,' Advanced Functional Materials 30 , 1903829 (2020).
- 9 W. H. Zachariasen, 'The atomic arrangement in glass,' Journal of the American Chemical Society 54 , 3841-3851 (1932).
- 8 A. K. Varshneya, Fundamentals of inorganic glasses (Elsevier, 2013).
- 10 G. Greaves, 'Exafs and the structure of glass,' Journal of Non-Crystalline Solids 71 , 203-217 (1985).
- 12 Y. V. Vavilov, V. E. Komaroc, and N. A. Tabunova, 'Calculation of the thermal conductivity of glasses,' The Soviet Journal of Glass Physics and Chemistry 8 , 326-330 (1982).
- 11 E. H. Ratcliffe, 'A survey of most probable values for the thermal conductivies of glasses between about -150 and 100°c, including new data on twnety-two glasses and a working formula for the calculation of conductivity from composition,' Glass Technology 4 , 113-128 (1963).
- 13 D. Pye, I. Joseph, and A. Montenero, Properties of glass-forming melts (CrC Press, 2005).
- 15 E. Martin, G. Ori, T. Q. Duong, M. Boero, and C. Massobrio, 'Thermal conductivity of amorphous sio2 by first-principles molecular dynamics,' Journal of Non-Crystalline Solids 581 (2022), 10.1016/j.jnoncrysol.2022.121434.
- 14 P. Jund and R. Jullien, 'Molecular-dynamics calculation of the thermal conductivity of vitreous silica,' Phys. Rev. B 59 , 13707-13711 (1999).

- 16 S. S. Sørensen, H. Johra, J. C. Mauro, M. Bauchy, and M. M. Smedskjaer, 'Boron anomaly in the thermal conductivity of lithium borate glasses,' Physical Review Materials 3 , 075601 (2019).
- 18 A. J. H. McGaughey and M. Kaviany, 'Phonon transport in molecular dynamics simulations: Formulation and thermal conductivity prediction,' Advances in Heat Transfer 39 , 169-255 (2006).
- 17 S. S. Sørensen, P. P. Cielecki, H. Johra, M. Bockowski, E. Skovsen, Y. Yue, and M. M. Smedskjaer, 'Thermal conduction in a densified oxide glass: Insights from lattice dynamics,' Materials Today Communications 32 , 104160 (2022).
- 19 M. Simoncelli, N. Marzari, and F. Mauri, 'Unified theory of thermal transport in crystals and glasses,' Nature Physics 15 , 809-813 (2019).
- 21 A. N. Cormack, J. Du, and T. R. Zeitler, 'Alkali ion migration mechanisms in silicate glasses probed by molecular dynamics simulations,' Physical Chemistry Chemical Physics 4 , 3193-3197 (2002).
- 20 L. Isaeva, G. Barbalinardo, D. Donadio, and S. Baroni, 'Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach,' Nature Communications 10 , 3853 (2019).
- 22 J. Du, 'Challenges in molecular dynamics simulations of multicomponent oxide glasses,' in Springer Series in Materials Science (Springer International Publishing, 2015) p. 157-180.
- 24 S. Plimpton, 'Fast parallel algorithms for short-range molecular dynamics,' Journal of Computational Physics 117 , 1-19 (1995).
- 23 A. Carré, S. Ispas, J. Horbach, and W. Kob, 'Developing empirical potentials from ab initio simulations: The case of amorphous silica,' Computational Materials Science 124 , 323-334 (2016).
- 25 J. Du and A. Cormack, 'The medium range structure of sodium silicate glasses: a molecular dynamics simulation,' Journal of Non-Crystalline Solids 349 , 66-79 (2004).
- 27 J. Young, F. Glaze, C. Faick, and A. Finn, Density of some soda-potashsilica glasses as a function of the composition (National Bureau of Standards, 1939).
- 26 G. White, J. Birch, and M. H. Manghnani, 'Thermal properties of sodium silicate glasses at low temperatures,' Journal of Non-Crystalline Solids 23 , 99-110 (1977).
- 28 V. F. Sears, 'Neutron scattering lengths and cross sections,' Neutron news 3 , 26-37 (1992).
- 30 G. Barbalinardo, Z. Chen, N. W. Lundgren, and D. Donadio, 'Efficient anharmonic lattice dynamics calculations of thermal transport in crystalline and disordered solids,' Journal of Applied Physics 128 (2020), 10.1063/5.0020443.
- 29 A. C. Wright, 'The comparison of molecular dynamics simulations with diffraction experiments,' Journal of Non-Crystalline Solids 159 , 264-268 (1993).
- 31 M. Bauchy, 'Structural, vibrational, and thermal properties of densified silicates: Insights from molecular dynamics,' The Journal of chemical physics 137 (2012).
- 33 S. Pailhès, H. Euchner, V. Giordano, R. Debord, A. Assy, S. Gomès, A. Bosak, D. Machon, S. Paschen, and M. de Boissieu, 'Localization of propagative phonons in a perfectly crystalline solid,' Physical Review Letters 113 (2014), 10.1103/physrevlett.113.025506.
- 32 S. R. Elliott, 'Physics of amorphous materials,' (No Title) (1990).
- 34 R. J. Bell and D. C. Hibbins-Butler, 'Acoustic and optical modes in vitreous silica, germania and beryllium fluoride,' Journal of Physics C: Solid State Physics 8 , 787-792 (1975).
- 36 J. Du and A. Cormack, 'The structure of erbium doped sodium silicate glasses,' Journal of Non-Crystalline Solids 351 , 2263-2276 (2005).
- 35 J. Du and A. N. Cormack, 'Molecular dynamics simulation of the structure and hydroxylation of silica glass surfaces,' Journal of the American Ceramic Society 88 , 2532-2539 (2005).
- 37 J. Du and L. R. Corrales, 'Compositional dependence of the first sharp diffraction peaks in alkali silicate glasses: A molecular dynamics study,' Journal of Non-Crystalline Solids 352 , 3255-3269 (2006).
- 39 F. Lodesani, M. C. Menziani, H. Hijiya, Y. Takato, S. Urata, and A. Pedone, 'Structural origins of the mixed alkali effect in alkali aluminosilicate glasses: Molecular dynamics study and its assessment,' Scientific Reports 10 (2020), 10.1038/s41598-020-59875-7.
- 38 J. Du and L. Kokou, 'Europium environment and clustering in europium doped silica and sodium silicate glasses,' Journal of Non-Crystalline Solids 357 , 2235-2240 (2011).
- 40 Y. Yu, M. Wang, M. M. Smedskjaer, J. C. Mauro, G. Sant, and M. Bauchy, 'Thermometer effect: Origin of the mixed alkali effect in glass relaxation,' Physical Review Letters 119 (2017), 10.1103/physrevlett.119.095501.
- 42 Q. Zhao, M. Guerette, G. Scannell, and L. Huang, 'In-situ high temperature raman and brillouin light scattering studies of sodium silicate glasses,' Journal of Non-Crystalline Solids 358 , 3418-3426 (2012).
- 41 M. Fábián, P. Jóvári, E. Sváb, G. Mészáros, T. Proffen, and E. Veress, 'Network structure of 0.7sio2-0.3na2o glass from neutron and x-ray diffraction and rmc modelling,' Journal of Physics: Condensed Matter 19 , 335209 (2007).
- 43 J. I. Gersten and F. W. Smith, The Physics and Chemistry of Materials (Wiley, 2001).
- 45 B. Mantisi, M. Bauchy, and M. Micoulaut, 'Cycling through the glass transition: Evidence for reversibility windows and dynamic anomalies,' Physical Review B - Condensed Matter and Materials Physics 92 , 134201 (2015).
- 44 H. Liu, S. H. Hahn, M. Ren, M. Thiruvillamalai, T. M. Gross, J. Du, A. C. van Duin, and S. H. Kim, 'Searching for correlations between vibrational spectral features and structural parameters of silicate glass network,' Journal of the American Ceramic Society 103 , 3575-3589 (2020).
- 46 D. Kilymis, S. Ispas, B. Hehlen, S. Peuget, and J.-M. Delaye, 'Vibrational properties of sodosilicate glasses from first-principles calculations,' Physical Review B 99 (2019), 10.1103/physrevb.99.054209.
- 48 W. Zhou, Y. Cheng, K. Chen, G. Xie, T. Wang, and G. Zhang, 'Thermal conductivity of amorphous materials,' Advanced Functional Materials 30 , 1903829 (2019).
- 47 H. R. Seyf, W. Lv, A. Rohskopf, and A. Henry, 'The importance of phonons with negative phase quotient in disordered solids,' Scientific Reports 8 (2018), 10.1038/s41598-018-20704-7.
- 49 M. T. Agne, R. Hanus, and G. J. Snyder, 'Minimum thermal conductivity in the context of: Diffuson -mediated thermal transport,' Energy and Environmental Science 11 , 609-616 (2018).
- 51 Y. Hiroshima, Y. Hamamoto, S. Yoshida, and J. Matsuoka, 'Thermal conductivity of mixed alkali silicate glasses at low temperature,' Journal of Non-Crystalline Solids 354 , 341-344 (2008).
- 50 P. M. Voyles, N. Zotov, S. M. Nakhmanson, D. A. Drabold, J. M. Gibson, M. M. J. Treacy, and P. Keblinski, 'Structure and physical properties of paracrystalline atomistic models of amorphous silicon,' Journal of Applied Physics 90 , 4437-4451 (2001).
- 52 S. Sukenaga, T. Endo, T. Nishi, H. Yamada, K. Ohara, T. Wakihara, K. Inoue, S. Kawanishi, H. Ohta, and H. Shibata, 'Thermal conductivity of sodium silicate glasses and melts: Contribution of diffusive and propagative vibration modes,' Frontiers in Materials 8 (2021), 10.3389/fmats.2021.753746.
- 54 W. Lv and A. Henry, 'Non-negligible contributions to thermal conductivity from localized modes in amorphous silicon dioxide,' Scientific Reports 6 (2016), 10.1038/srep35720.
- 53 K. Chojin, M. Shimizu, Y. Shimotsuma, and K. Miura, 'Cooling-rate dependence of thermal conductivity in a sodium silicate glass: A molecular dynamics study,' Journal of the Ceramic Society of Japan 128 , 656-659 (2020).

## SUPPLEMENTARY MATERIAL Thermal conductivity in modified oxide glasses is

## governed by modal phase changes

P. Rasmussen 1 and S. S. Sørensen 1

1 Department of Chemistry and Bioscience, Aalborg University, Aalborg, Denmark, | soe@bio.aau.dk https://orcid.org/0000-0003-2230-7823

Table S 1: Molar concentration of Na2O ( x ) in the x Na2O-(100- )SiO2 x glasses, mass density ( ρ m ), atomic density ( ρ A ), average molar mass ( Mavg ), simulation cell size, l cell and number of atoms in each simulation cell used in simulations containing 3000 atom box sizes. Mass densities are adapted from Refs. [1], [2].

| x (mol%)   | ρ m (g cm - 3 )   | ρ A (Å - 3 · 10 - 2 )   | M avg (g mol - 1 )   | l cell (Å)   | # of atoms in simulation cell   | # of atoms in simulation cell   | # of atoms in simulation cell   |
|------------|-------------------|-------------------------|----------------------|--------------|---------------------------------|---------------------------------|---------------------------------|
|            |                   |                         |                      |              | Na                              | O                               | Si                              |
| 0          | 2.2               | 6.62                    | 20.028               | 35.66        | 0                               | 2000                            | 1000                            |
| 10         | 2.289             | 6.86                    | 22.091               | 35.23        | 200                             | 1900                            | 900                             |
| 15         | 2.336             | 6.99                    | 20.123               | 35.01        | 300                             | 1850                            | 850                             |
| 20         | 2.383             | 7.17                    | 20.154               | 34.80        | 400                             | 1800                            | 800                             |
| 25         | 2.431             | 7.25                    | 20.186               | 34.58        | 500                             | 1750                            | 750                             |
| 30         | 2.466             | 7.35                    | 20.217               | 34.44        | 600                             | 1700                            | 700                             |
| 33.3       | 2.491             | 7.41                    | 20.238               | 34.33        | 666                             | 1667                            | 667                             |
| 35         | 2.497             | 7.43                    | 20.249               | 34.31        | 700                             | 1650                            | 650                             |
| 40         | 2.532             | 7.52                    | 20.28                | 34.17        | 800                             | 1600                            | 600                             |
| 45         | 2.541             | 7.53                    | 20.312               | 34.15        | 900                             | 1550                            | 550                             |
| 50         | 2.560             | 7.58                    | 20.344               | 34.01        | 1000                            | 1500                            | 500                             |

Table S 2: Molar concentration of Na2O ( x ) in the x Na2O-(100- )SiO2 x glasses, mass density ( ρ m), atomic density ( ρ A), average molar mass ( M avg ), simulation cell size, l cell and number of atoms in simulations used for estimating thermal properties with the Quasi-harmonic Green-Kubo approach (600 atoms). Mass densities are adapted from Refs. [1], [2].

| x (mol%)   | ρ m (g cm - 3 )   | ρ A (Å - 3 · 10 - 2 )   | M avg (g mol - 1 )   | l cell (Å)   | # of atoms in simulation cell   | # of atoms in simulation cell   | # of atoms in simulation cell   |
|------------|-------------------|-------------------------|----------------------|--------------|---------------------------------|---------------------------------|---------------------------------|
|            |                   |                         |                      |              | Na                              | O                               | Si                              |
| 20         | 2.383             | 7.17                    | 20.154               | 20.35        | 80                              | 360                             | 160                             |
| 25         | 2.431             | 7.25                    | 20.186               | 20.22        | 100                             | 350                             | 150                             |
| 30         | 2.466             | 7.35                    | 20.217               | 20.14        | 120                             | 340                             | 140                             |
| 33.3       | 2.491             | 7.41                    | 20.238               | 20.08        | 132                             | 335                             | 133                             |
| 35         | 2.497             | 7.43                    | 20.249               | 20.07        | 140                             | 330                             | 130                             |

Table S 3: Average bond lengths ( ) of Si-O and Na-Ointeractions , and coordination l numbers (CN) of Na and Si of simulated sodium silicate glasses at various sodium concentrations ( x ). Experimental data is provided as extracted from neutron diffraction studies. References are provided with the experimental values.

|   x (mol%) |   l Si-O (Å) | l Na-O (Å)   | l Si-O,exp (Å)   | l Na-O,exp (Å)   | CN Na ( - )   |   CN Si ( - ) | CN Si-O,exp ( - )   | CN Na-O,exp ( - )   |
|------------|--------------|--------------|------------------|------------------|---------------|---------------|---------------------|---------------------|
|          0 |       1.6125 |              | 1.610            |                  |               |             4 | 4.0 [3]             |                     |
|         10 |       1.6125 | 2.3625       |                  |                  | 5.27          |             4 |                     |                     |
|         15 |       1.6125 | 2.3775       |                  |                  | 5.86          |             4 |                     |                     |
|         20 |       1.6125 | 2.3925       | 1.617            | 2.36             | 6.12          |             4 |                     |                     |
|         25 |       1.6125 | 2.3775       |                  |                  | 6.09          |             4 |                     | 5 [4]               |
|         30 |       1.6125 | 2.3775       | 1.622            | 2.36             | 5.89          |             4 |                     |                     |
|         33 |       1.6125 | 2.3775       |                  |                  | 5.79          |             4 | 4.0 [3]             | 5.5 [3]             |
|         35 |       1.6125 | 2.3775       |                  |                  | 6.08          |             4 |                     |                     |
|         40 |       1.6125 | 2.3775       |                  |                  | 6.14          |             4 |                     |                     |
|         45 |       1.5975 | 2.3775       |                  |                  | 6.06          |             4 |                     |                     |
|         50 |       1.5975 | 2.3775       |                  |                  | 5.98          |             4 |                     | 5 [4]               |

Table S 4: Elastic properties of glasses investigated, where x is the molar concentration of Na2O ( x ) in the x Na2O-(100- )SiO2 x glasses, E is the Young's modulus, B is the bulk modulus, G is the shear modulus and v is the Poisson's ratio. Data provided are estimated from 3000 atom system sizes.

|   x (mol%) |   E (GPa) |   B (GPa) |   G (GPa) |   v (-) |
|------------|-----------|-----------|-----------|---------|
|          0 |     101.3 |      65.5 |      40.8 |  0.2415 |
|         10 |      79.2 |      51.2 |      31.9 |  0.2422 |
|         15 |      72.2 |      47.1 |      29   |  0.2441 |
|         20 |      66.8 |      43.6 |      26.8 |  0.2441 |
|         25 |      63   |      40.6 |      25.4 |  0.2402 |
|         30 |      56.9 |      38.8 |      22.7 |  0.2553 |
|         33 |      56.2 |      39.9 |      22.3 |  0.2637 |
|         35 |      53.9 |      40.2 |      21.1 |  0.2756 |
|         40 |      52.3 |      39.7 |      20.4 |  0.2794 |
|         45 |      49.1 |      37.4 |      19.2 |  0.2811 |
|         50 |      45.7 |      39.4 |      17.5 |  0.3053 |

Table S 5: Thermal conductivities data computed with the Quasi-Harmonic Green-Kubo method of x Na2O-(100- )SiO2 x glasses where ( x ) is the molar concentraion of Na2 O, ( κ ) is the average thermal conductivity and (SD) is the standard deviations of estimated κ values.

|   x (mol%) |   κ (W m - 1 K - 1 ) |   SD (W m - 1 K - 1 ) |
|------------|----------------------|-----------------------|
|         20 |                1.185 |                 0.006 |
|         25 |                1.069 |                 0.029 |
|         30 |                0.986 |                 0.015 |
|         33 |                0.986 |                 0.012 |
|         35 |                0.915 |                 0.007 |

## Figure S1

Figure S 1: Potential energy curve for the Na-O, O-O and Si-O interactions of the Teter potential plotted without the r 24 short term correction (dashed lines) and with the short term correction (solid lines), illustrating the modification added to avoid the "Buckingham catastrophe".

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-30ae90d4-a7ae-450c-8e27-3683ecc4c791.jpg)

## Figure S2

Figure S 2: Mass densities ( ρ m) of x Na2O-(100- )SiO2 x glasses at room temperature used in simulations. Densities are adapted from Refs. [1], [2].

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-c217bd81-e3be-4478-8afe-4c6f6c3379ec.jpg)

## Figure S3

Figure S 3: Experimental structure factor from neutron diffraction experiments and simulated structure factor of SiO 2 glass using the potential of Teter. Based on the R χ factor (see Methods in the main text), the medium range order shows good agreement with experimental data (generally R &lt;10% is considered good), yet lower agreement than that of the simulated sodium silicate glasses χ (Figure 1a in the main text).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-bc226811-9d3d-46c4-9ec7-27b486cb7bf8.jpg)

## Figure S4

Figure S 4: Radial distribution functions (RDF's) of x Na2O-(100- )SiO2 x glasses with varying sodium concentrations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-aa1fe876-f056-4703-ae33-545a73a91649.jpg)

## Figure S5

Figure S 5: Partial distribution functions (PDF's) of x Na2O-(100- )SiO2 x glasses with varying sodium concentrations of pair-interactions (a) Na-O, (b) Na-Na, (c) O-O, (d) Na-Si, (e) Si-O and Si-Si. Insets are provided to show peak intensities.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-ee49a84b-0905-422d-9184-6cf7c078cf52.jpg)

## Figure S6

Figure S 6: Angular distribution functions (ADF's) of (a) inter-tetrahedral Silicon-Bridging oxygenSilicon bond angles and (b) O-Si-O intra-tetrahedral bond angle distributions with increasing sodium concentration.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-8a2f5af8-62cf-4655-bc5e-2f707e9d0169.jpg)

## Figure S7

Figure S 7: Angular distribution functions (ADFs) of (a) the Si-non-bridging oxygen-Na angle formed by non-bridging oxygens (NBOs) between silicon and sodium atoms and (b) the O-Na-O angles formed in the network modifier environments.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-1396eb4c-eaa4-4150-96b9-f4677a683f6b.jpg)

## Figure S8

Figure S 8: Partial Faber-Ziman structure factors of pair interactions (a) Na-O, (b) O-O and (c) Si-O

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-4b9f14a5-5247-44e5-be78-61804fa7c7d6.jpg)

## Figure S9

Figure S 9: Partial Faber-Ziman structure factors of cationic interactions (a) Na-Na, (b) Na-Si and (c) Si-Si.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-30c2d321-af7a-4c49-b5c0-32f136613082.jpg)

## Figure S10

Figure S 10: Neutron-weighted structure factor of simulated glasses where x is the molar concentration of Na2O ( x ) in the x Na2O-(100- )SiO2 x glasses. The dashed arrow is shown to emphasize the peak and positional change of the first sharp diffraction peak (FSDP), that is, the scattering peak at the lowest value of Q .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-4fd44432-d02a-4c0a-b9f5-eb5bd342f8eb.jpg)

## Figure S11

Figure S 11: Simulated and experimental data of (a) Q n -distributions in glasses of 3000 atom systems and (b) the distribution of bridging and non-bridging oxygen (BO and NBO, respectively) atoms. Lines in both plots are shown as guides for the eye.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-669934f4-1d34-4fe5-8590-a67dfe58d4ed.jpg)

## Figure S12

Figure S 12: Simulated and experimental Poisson's Ratio with increasing sodium content in the studied x Na2O-(100- )SiO2 x glasses. Simulations are somewhat able to predict the right trend, yet the increase in Poisson's ratio above 10 mol% Na2O in the experimental is less pronounced in simulated data and values generally deviate more than the other moduli data.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-0f326cd3-b957-4bb4-b1a3-c519178efe74.jpg)

## Figure S13

Figure S 13: Simulated vibrational density of states (VDOS) of (a) the 33 3Na2O-66 7SiO2 . . glass and (b) SiO2 glass. Simulated VDOS (red) is compared with ab initio calculations (dashed blue) from ref. [5] showing good agreement for the sodium silicate glass, but poor agreement for SiO2 glass when simulated with the potential of Teter.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-a85db548-e7a8-44dd-8d57-1906f73e2c0b.jpg)

## Figure S14

20Naz

0.05

0.04

0.03

0.01

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-8bcd8587-ef7d-4572-b461-d97fdd939e3b.jpg)

10

15

20

25

30

35

40

(THz)

Figure S 14: Simulated vibrational density of states (VDOS) of (a) the 20Na2O-80SiO2 glass and (b) 25Na2O-75SiO2 glass. Simulated VDOS (red) is compared with ab initio calculations (dashed blue) from ref. [5] showing fairly high agreement for both sodium silicate glass systems.

## Figure S15

Figure S 15: Partial vibrational density of states (pVDOS) of oxygen, sodium, silicon and further subdivisions of bridging and non-bridging oxygen (BO and NBO) as well as based on Q n value, showing the frequency ranges where different atomic species contribute to the vibrational density of states.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-5ba4329c-34a5-481e-8a7f-f88b77a6bb94.jpg)

## Figure S16

Figure S 16: (Partial) vibrational density of states (VDOS) of the 33 3Na2O-66 7SiO2 . . system illustrating the contribution to the total VDOS from atomic weighted vibrational density of states from various atomic species. We note that the contribution of O is also encompassed in the BO and NBO contributions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-3911cc20-a8e7-4457-8077-f9d4f27680ac.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-b9c93b4f-6d12-45e8-8556-38fbf369fc50.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-685ea575-c3ad-4712-b107-b8360482f5c3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-ff50c6ec-a943-496c-b9f7-b1baa262b4ff.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-67f01142-e0ac-4935-9943-5dd58be1b372.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-344cac39-451d-48b6-90cf-1aa4c2ab2716.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-4196193b-c1d7-4e77-a11c-a4dc0319dd36.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-578b9c8d-e2fd-4ac8-a087-c0f00efbd0c0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-2d8abd8f-f746-4ae4-b564-53cc518ae9b4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-cb300a98-a942-47ad-9d04-532ac4741e3d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-1fb2f26e-707f-49c8-b85e-a07bd0316de7.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-ecc40b4f-5634-41ed-94af-e41d14db2813.jpg)

Figure S17

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-b97daf65-ef65-4937-bea9-75c923d7d506.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-1c6936d4-1101-4c64-adda-4a3662237c89.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-ae46a6ac-b1ab-45fb-970c-93e59403ad33.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-a749b006-6b72-4100-a0fd-28ca0a35ed9c.jpg)

Figure S 17: Figure illustrating the changes of eigenvectors at different frequencies in glasses with increasing sodium content. Observing the low-frequency illustrations, i.e. 1 and 5 THz, eigenvectors seem to become more evenly distributed which correlates with an increasing participation ratio. The eigenvectors at 1 THz seem to become less periodic with increasing sodium contents, which could be a potential reason for the decreased per-mode heat transfer in this range, despite the increased participation ratio. In the mid-frequency range from 15-20 THz, eigenvector lengths becomes shorter in specific spatial regions of the glasses, meaning that modes become more localized with increasing sodium content. This correlates with the decreased participation ratio in the observed ranged. Lastly, high-frequency modes seem become even more localized when increasing the sodium content, suggesting that network depolymerization induces high-frequency mode localiztion.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-2719d142-e4dc-4260-a7da-33ef228033f5.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-914332e4-d944-4f15-be5d-585c0def8291.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-eb6fa08e-3105-4614-a580-0698df4e8d70.jpg)

## Figure S18

Figure S 18: Figure illustrating (a) participation ratio (PR) and (b) normalized inverse participation ratio (IPR) of simulated x Na2O-(100- )SiO2 x with increasing sodium content.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-228146bb-d0c5-48ef-9dad-3baae9096d95.jpg)

## Figure S19

Figure S 19: Figure illustrating the total participation ratio (top left), and the atomic participation ratio (APR) of various atomic species and further subdivisions of bridging and non-bridging oxygen (BO and NBO) as well as based on Q n value of the studied x Na2O-(100- )SiO2 x glasses with increasing sodium content.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-3cde0b4e-cc69-4ff8-8425-7ee633fe4905.jpg)

## Figure S20

Figure S 20: Phase quotient (PQ) of simulated x Na2O-(100- )SiO2 x glasses with increasing sodium contents ( x ). The dashed line indicates the transition point from in-phase to out of phase vibrations ( PQ = 0). The main plot is fitted to the data in the inset to provide easier comparison between the molar compositions of glasses.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-80afdd2c-73eb-4501-811f-b634d28e8542.jpg)

## Figure S21

Figure S 21: Figure illustrating (a) the total phase quotient (PQ) and the atomic phase quotient (APQ) of (b) sodium, (c) silicon and (d) oxygen, (e) bridging-oxygens, and (f) non-bridging oxygens in simulated x Na2O-(100- )SiO2 x glasses with increasing sodium content.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-fdad5e84-c334-49c9-aaa9-6e96f0da6fa1.jpg)

## Figure S22

Figure S 22: Modal diffusivity, Di , of simulated x Na2O- 100-( x ) SiO2 glasses with increasing sodium content.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-8f26beec-db97-49ba-98a2-0dbce24d56e8.jpg)

## Figure S23

Figure S 23: Effect of sodium concentration on (a) modal participation ratio and normalized atomic participation ratio of (b) oxygen, (c) silicon, (d) bridging-oxygens, (e) sodium and (f) non-bridging oxygens coupled with modal thermal conductivities estimated with the QHGK method.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00813v1-figure-853f8842-658e-4775-a086-207874df1eb3.jpg)

## Bibliography

- [1] Q. Zhao, M. Guerette, G. Scannell, and L. Huang, 'In-situ high temperature raman and brillouin light scattering studies of sodium silicate glasses,' Journal of Non-Crystalline Solids , vol. 358, no. 24, pp. 3418-3426, Dec. 2012, I S S N: 0022-3093.
- [2] J. Young, F. Glaze, C. Faick, and A. Finn, Density of some soda-potash-silica glasses as a function of the composition . National Bureau of Standards, 1939.
- [3] M. Misawa, D. Price, and K. Suzuki, 'The short-range structure of alkali disilicate glasses by pulsed neutron total scattering,' Journal of Non-Crystalline Solids , vol. 37, no. 1, pp. 85-97, 1980.
- [4] A. C. Hannon, S. Vaishnav, O. L. Alderman, and P. A. Bingham, 'The structure of sodium silicate glass from neutron diffraction and modeling of oxygen-oxygen correlations,' Journal of the American Ceramic Society , vol. 104, no. 12, pp. 61556171, 2021.
- [5] D. Kilymis, S. Ispas, B. Hehlen, S. Peuget, and J.-M. Delaye, 'Vibrational properties of sodosilicate glasses from first-principles calculations,' Phys. Rev. B , vol. 99, p. 054 209, 5 Feb. 2019. DOI: 10.1103/PhysRevB.99.054209. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.99.054209.
Philip Rasmussen Søren Strandskov Sørensen
Thermal conductivity in modified oxide glasses is governed by modal phase changes
The thermal conductivity of glasses is well-known to be significantly harder
to theoretically describe compared to crystalline materials. Because of this
fact, the fundamental understanding of thermal conductivity in glasses remain
extremely poor when moving beyond the case of simple glasses, e.g., glassy
SiO$_2$, and into so-called 'modified' oxide glasses, that is, glasses where
other oxides (e.g. alkali oxides) have been added to break up the network and
alter e.g. elastic and thermal properties. This lack of knowledge is apparent
despite how modified glasses comprise the far majority of known glasses. In the
present work we study an archetypical series of sodium silicate
($x\text{Na}_2\text{O}\text{-}(100\text{-}x)\text{SiO}_2$) glasses. Analyses of
modal contributions reveal how increasing Na$_2$O content induces increasing
vibrational localization with a change of vibrations to be less ordered, and a
related general decrease in modal contributions to thermal conductivity. We
find the vibrational phases (acoustic vs. optical) of sodium vibrations to be
relatively disordered compared to the network-forming silicon and oxygen
species, explaining how increasing Na$_2$O content decreases thermal
conductivity. Our work sheds new light on the fundamentals of glassy heat
transfer as well as the interplay between thermal conduction and modal
characteristics in glasses.