In a liquid, particles are free to flow. In contrast, a solid forms because the energy of the particles decreases, allowing them to take on a relatively ordered, three-dimensional, and more rigid structure. Even so, particles can flow within a solid, given the right circumstances. For example, in glasses, monovalent cations have been found to diffuse throughout their matrix relatively quickly. Some glass compositions have already been studied for different applications because of their remarkably high alkali ion diffusion, especially under an external electric field (Daiko et al., 2022). The interest in this particular situation goes beyond scientific curiosity. These materials may be applied in solid-state batteries and sensors, and the phenomenon is closely related to the mixed-alkali effect (MAE) (Calahoo et al., 2020).

In the case of a single charge carrier, the diffusion under the influence of an external electric field at a temperature T can be given by the ionic conductivity ( σ ) using the Nernst–Einstein equation (Varshneya and Mauro, 2019): σ = Z · e 2 · n · D k B · T , (1) where Z is the valence of the charge carrier (in the case of a monovalent cation, Z = 1 ) , e is the elementary charge, n is the density of effective charge carriers, D is the diffusion coefficient, k B is Boltzmann’s constant, and T is the absolute temperature.

For a given temperature and chemical composition, σ is entirely dependent on n and D . However, determining which parameters control the phenomenon can be complex, as obtaining the terms empirically and independently from each other is challenging with the currently available technology. Consequently, no universally accepted model explains the mechanisms behind the ionic conductivity in all glass materials (Dyre et al., 2009). Generally, the models found in the literature can be divided into strong electrolytes or weak electrolytes, depending on whether the ions are entirely mobile in the glass matrix or partially immobilized (Bragatto, 2020). This distinction implies that n might be equal to the total density of the species responsible for the ionic conductivity or just a fraction of the total density. Therefore, an essential step in developing such a universal model is to gain a better understanding of the nature of a charge carrier.

One way to look at this problem is by analyzing Haven’s ratio ( H R ) (Varshneya and Mauro, 2019). H R is defined as the ratio of the self-diffusion coefficient of ions ( D S e l f ) to the diffusion coefficient calculated from conductivity ( D σ ) using Eq. 1: H R = D S e l f D σ . (2)

According to Haven and Verkerk (1965), H R can be interpreted as follows: for an ionic crystal in which the diffusion is interstitial, the values for D S e l f and D σ are the same, and H R = 1. On the other hand, if the mobile ion follows a different random walk, in which some of the jumps are not allowed during D S e l f but are allowed during D σ , the diffusion coefficient will be different, with D σ > D S e l f . A more in-depth analysis of H R along with its meaning and impact on the mechanisms of ionic conduction in glasses can be found in the research by Isard (1999), Kahnt (1996), and Murch (1982).

In their work, Haven and Verkek found the values of H R were between 0.4 and 0.6 for different sodium silicate glasses. The authors also observed that H R decreased sharply with increasing composition up to 10 mol% and kept dropping more slowly at higher concentrations. Since their work, many authors have investigated H R , obtaining similar values for other ionic-conductive glasses.

D S e l f can be measured in a laboratory using radioactive tracers. However, calculating D σ from Eq. 1 might be problematic as the true nature of n must be known. In general, n is taken as the total ion concentration per volume, implying that the charge carriers act as a strong electrolyte. Molecular dynamics (MD) is a novel and powerful tool that can shed some light on this problem. Previous studies (Welch et al., 2019; Atila et al., 2020; Zhao et al., 2020) were able to measure the diffusion of different species of glassy materials using MD without making assumptions about the nature of the charge carriers. In other words, D can be understood from the mean squared displacement ( r 2 ) of mobile ions by D = r 2 2 · n d · ∆ t , (3) where n d is the dimensionality within which the process occurs and ∆ t is the observational time. In turn, r 2 is given by r 2 = ∆ x → t 2 = 1 N ∑ i = 1 N x → i t − x → i 0 2 , (4) where N is the number of averaged particles, x → i 0 is an initial reference position, and x → i t is the position with respect to time.

Although experiments on the diffusion of species in glasses are common, the study of the same diffusion with an external electric field is not. In this work, glasses of the composition xA₂O⋅(1-x)SiO₂, with x = 0.1, 0.2, and 0.3, and A = Li, Na, and K were prepared using MD, and the r 2 values of the alkali ions with and without an external field were calculated. These values were used to obtain D S e l f , D σ , and H R . We hypothesize that studying these parameters with MD and comparing them to laboratory experimental results and interpretations in the literature will validate the use of this technique and expand the current application of MD in the glass science field to study the material under the effects of an external electric field.

All simulations were run using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), an open-source molecular dynamics software application (Du and Cormack, 2022). Alkali silicate glasses of the form xA₂O∙(1-x)SiO₂ were simulated for x = 0.1, 0.2, and 0.3 and A = Na, Li, and K. A visual representation can be seen in Figure 1.

Pairwise potentials developed by Pedone et al. (2006) were used to describe the interactions between the species according to U r = Z i Z j e 2 r + D i j [ { 1 − e − a i j r − r 0 } 2 − 1 ] + C i j r 12 , (5) where D _ij is the bond dissociation energy, a _ij is a function of the slope of the potential energy well, and r ₀ is the equilibrium bond distance. The parameters derived by Pedone et al. from binary oxide crystals pertinent to the systems in this work are given in Table 1. These potentials were chosen due to their widely reported accuracy for silicates (Pedone et al., 2008; Pedone, 2009; Jabraoui et al., 2016).

Simulation volumes were chosen to match experimentally obtained density values for each composition. The density values used to find the volume of the simulation cell for each glass composition from this work can be found in Table 2 ¹⁸.

A time step of 1 femtosecond (fs) was used for each composition. All simulations were split into two steps. The first step was the creation of the simulated glass, and the second step included testing each composition at various temperatures with and without an applied electric field. A benefit of splitting the simulations into steps was that an almost identical glass structure was present for each simulation. Each composition was simulated independently in phase one, with 5100 atoms correctly proportioned. Atoms were placed randomly inside the simulation cell, with densities, shown in Table 2, determining the volume.

In all simulations, a microcanonical ensemble that constrained the number of particles, volume, and internal energy (NVE) was used to facilitate energy distribution in the simulation cell. After initialization with an NVE ensemble for 100 ps, the simulations were introduced to a temperature of 3000 K in an NVT ensemble (constant number of particles, volume, and temperature) for 100 ps. Next, with energy equilibrium being reached, the pressure was allowed to equilibrate through the NPT ensemble (constant number of particles, pressure, and temperature) for 100 ps. The glasses were then cooled to 500 K at a rate of 1 K/ps. Here, each glass composition was duplicated among the various temperature and electric field tests.

Phase two started with another NVE initialization process for 100 ps. Then, the glass was raised to T at 1 K/ps in an NPT ensemble where T = 800 K, 1000 K, 1500 K, 2000 K, 2500 K, and 3000 K. Finally, an external electric field of 0.05 V/Å was applied across the x dimension of the simulation cell. This field strength was chosen as lower values had no detectable difference between exposed and unexposed samples. A schematic representation of the experiment is shown in Figure 2. In comparison, higher values resulted in unrealistic values for MSD and diffusion of the less mobile atoms, possibly resulting from dielectric breakdown. Similar electric field strengths were observed in other simulation works (Welch et al., 2019).

Another valid and more common concern in this kind of work is the time length of the diffusional process (Du and Cormack, 2022). Due to molecular dynamic limitations, the observational time is very short compared to laboratory experiments. Therefore, in our calculations of diffusion coefficient, we discarded the first 100 ps of every run to avoid any non-diffusive behavior by the ions.

For each glass composition and each temperature, six different experiments were completed, measuring r 2 as a function of time. Three were completed without any electric field applied, and three experiments applied an electric field of 0.05 V/Å on the x-axis. An example can be seen in Figure 3, with results for the glass 0.3 Na₂O 0.7 SiO₂ at 1000 K.

To obtain reliable values for the diffusion coefficient D , the diffusion must be at a steady-state regime (Varshneya and Mauro, 2019), meaning that the dependency of r 2 with time should be linear. To determine if the regime could be considered at a steady-state pace, all data were fitted using OriginPro 2019 software, and all results had R ² ≥ 0.997 after disregarding data points obtained before the 100 ps. Values of r 2 ∆ t were obtained from these fittings, and values for D S e l f and D σ were calculated using Eq. 3. Also, it was possible to obtain values of H R using Eq. 2. An example of these results for measurements determined at 1000 K can be found in Table 3.

Values for diffusion coefficients found in this work were comparable to laboratory diffusion coefficients for sodium in soda-lime silicate glasses (Mehrer et al., 2008), previous simulation results found by Welch et al. (2019), and thermal diffusion values found using MD in the literature (Du and Chen, 2012; Du and Cormack, 2022). This should be taken with caution because an artifact of the molecular dynamics experiment is that the glass’s temperature and the fictive temperature might be underestimated (Du and Cormack, 2022).

Values for activation energy ( E A ) can be obtained from the dependency of the diffusion coefficient with temperature, assuming an Arrhenius behavior (Varshneya and Mauro, 2019): D = D 0 + exp − E A k B T , (6) where D can be either D S e l f or D σ , and D 0 is a constant. Results for this work can be found in Figure 4, expressed as log ⁡ D as a function of 1000 / T . The activation energies were obtained by fitting the data in Figure 4 and can be found in Table 4.

Activation energy values for the self-diffusion were comparable to values found in the literature for lithium silicate glasses using MD (Du and Chen, 2012).

The data obtained to calculate Haven’s ratio for all glass compositions at the temperature interval range described earlier are shown in Figure 4. The results can be divided into two categories: composition dependence and temperature dependence.

Haven’s ratio and its origin are not yet fully understood, and obtaining good and reliable results can be challenging as it requires assumptions about the nature of ionic conductivity and expensive measurements for self-diffusion. However, as a relatively new technique, molecular dynamics can be a fast, cheap, and easy addition to research, given the proper considerations. In this work, the authors have examined simple alkali silicate glasses to verify the validity of MD techniques in this context.

Due to its nature, MD simulation experiments have limitations, especially with sample size and observational times, both relevant to the ionic conductivity phenomenon. Therefore, to observe a response from the system, the temperatures and voltage field applied had to be comparably higher than in laboratory experiments. The shift to higher temperatures is well-known in the literature. Still, the higher electric field used in this work was found by trial and error, and it was enough to create a significant diffusion of the alkali studied but not affect the diffusion of the glass-former units.

Values of r ² were obtained using the LAMMPS software, and from them, values for D_self and D_σ were calculated. Their dependency on temperature, alkali concentration, and alkali nature was also studied. Results obtained in this work agree with different points found in the literature, such as the compositional dependency of H_R, where the value sharply decreases at concentrations above x = 0.1 of alkali oxide, and temperature dependence, in which H_R approaches unity at higher temperatures. These results are encouraging and show the possibility of using MD techniques to further investigate the diffusional properties of oxide glasses, especially when under an external electric field.