# Evidence and origin of anomalous diffusion of Ag + ion in amorphous silica: a molecular dynamics study with neural network interatomic potentials

Salomé Trillot, Nathalie Tarrat, Nicolas Combe, Patrizio Benzo, Caroline

Bonafos, Magali Benoit

► To cite this version:

Salomé Trillot, Nathalie Tarrat, Nicolas Combe, Patrizio Benzo, Caroline Bonafos, et al.. Evidence and origin of anomalous diffusion of Ag + ion in amorphous silica: a molecular dynamics study with neural network interatomic potentials. The Journal of Chemical Physics, inPress, Special Issue "Molecular Dynamics, Methods and Applications 60 Years after Rahman". hal-04973509

HAL Id: hal-04973509

https://hal.science/hal-04973509v1

Submitted on 3 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Evidence and origin of anomalous diffusion of Ag⁺ ion in amorphous silica: a molecular dynamics study with neural network interatomic potentials

Salomé Trillot,¹ Nathalie Tarrat,¹, a) Nicolas Combe,¹ Patrizio Benzo,¹ Caroline Bonafos,¹ and Magali Benoit¹, b)

CEMES-CNRS and Université de Toulouse, 29 rue Jeanne Marvig, 31055 Toulouse Cedex,
France

(Dated: March 3, 2025)

The release of Ag⁺ ions into the environment through silica layers is a promising strategy for the development of anti-microbial surface coating devices. The aim of the present study is to provide some insight into the elementary mechanisms of Ag⁺ ions diffusion through silica with the objective of proposing control strategies. Thanks to the development of interaction potentials based on neural networks, the diffusion processes were studied via molecular dynamics simulations. Silver diffusion was found to be anomalous and sub-diffusive, the origin of which could be attributed to deceleration and temporal anti-correlations. This sub-diffusion has been attributed primarily to the disordered nature of the silica matrix. Furthermore, it is magnified by the presence of coordination defects within the silica matrix. These defects, in particular the under-coordinated oxygen atoms, act as traps for Ag⁺ by forming O-Ag bonds, thereby limiting the jump length and retaining the ion for long duration. By comparison with existing diffusion models, the diffusion mechanism in the absence of defects appears to be of the Fractional Brownian motion type, substantially modified by the presence of defects. Two possible approaches have emerged to tune the Ag⁺ ions release through the silica layer: the monitoring of the number of defects and the opening/closing of diffusion paths via e.g. a modification of the silica density.

PACS numbers: 61.43.Fs; 66.30.hh; 66.30.J-; 02.70.Ns; 71.15.Mb; 61.43.Bn

## I. Introduction

Silica in its amorphous form is very often used for nanoparticles (NPs) coating in the field of nanotechnology, materials science, and biology¹, not only for its surface characteristics², but also for protection purposes. Indeed, it protects the NPs from the effects of the environment, thus keeping their physical properties intact, and conversely, protects the environment from the toxic effects of the NPs. Moreover, in the case of noble metal NPs, the silica layer prevents their aggregation³,⁴ and, as it is optically transparent, the optical properties of the coated NPs remain exploitable. It is particularly interesting in the case of NPs exhibiting localized surface plasmonic resonance, such as silver ones⁵,⁶.

Silica coating of silver NPs (AgNP) can be experimentally achieved in several ways, including chemical synthesis⁷, micro-emulsion routes⁸, flame aerosol⁹, heat treating of composite SiO₂/Ag powders¹⁰, consolidation of a Ag-attached silica colloid¹¹, combination of sputtering and plasma polymerization¹² and low energy ion beam synthesis¹³. Interestingly, despite the presence of the protective silica layer, these synthesized coated NPs can affect a biological medium in which they are placed¹³⁻¹⁵. Experimental evidences have proven that the observed effects do not come from the NP directly, but from metal ions which detach from the NP and diffuse through the silica layer to reach the surrounding medium¹³,¹⁵,¹⁶. Since silver ions are known for their ability to inhibit microbial

colonization¹⁷,¹⁸, such silver ion release properties open the way for the design of antimicrobial and antibiofilm surface coating strategies.

In previous works, we took advantage of these properties to develop bactericidal surface coating devices¹²,¹³ consisting of AgNPs embedded a few nanometers below the surface of a silica layer. The NPs act as a reservoir of Ag⁺ ions and the silica coating, in addition to its protective properties, can be used to modulate the release of Ag⁺ ions into the environment. Up to date, we have shown experimentally that, depending on the distance at which the NPs are embedded beneath the surface, the Ag⁺ ions release is activated or not¹³. A consequent step forward would be to control the release rate, which implies the need of understanding the mechanisms involved at the atomic scale. However, a number of questions remain unanswered, notably concerning the key factors driving the quantity of ions released from the AgNPs and their diffusion paths and rate in such devices. Answering these questions is of paramount importance for the development of effective devices.

To probe the nature of the interface between the AgNPs surface and the surrounding silica, we undertook theoretical studies using DFT-type numerical simulations, first in the form of thin films¹⁹ and then as NPs embedded in a silica matrix²⁰. These investigations enabled us to show that (i) Ag-O and Ag-Si bonds spontaneously form at the AgNP-silica interface, and (ii) AgNP loose electrons that are captured by the surrounding silica, and become positively charged. The formation of such bonds could be the first step in the detachment of silver ions (Ag⁺) from the surface of the NPs. This release could be assisted by species entering the silica, such as oxygen molecules, which are known to be one of the main causes of dissolution of AgNPs placed

<footnote>a) nathalie.tarrat@cemes.fr</footnote>

<footnote>b) magali.benoit@cemes.fr</footnote>
in aqueous and biological media²¹⁻²³. The fact that the Ag-NPs are positively charged when embedded in silica could also play in favor of cation detachment. Once detached from the NP, the diffusion nature and mechanisms of the Ag⁺ ions through the silica matrix remain completely unknown at present, as well as the parameters which could influence it.

The aim of the present study is to provide some insight into the elementary mechanisms of Ag⁺ ions diffusion in an amorphous silica matrix. To address this question from a theoretical standpoint, a major limitation is the size of the systems that can be simulated using DFT methods, as well as the simulation times achievable with such first-principles methods. To circumvent this limitation, we have developed a machine learning potential based on a neural network (NNP) trained with data from DFT calculations. This enabled us to perform molecular dynamics trajectories over relatively long times and distances. By analyzing the nature of diffusion and the underlying elementary mechanisms involved, suggestions can be made as to possible ways of controlling diffusion. In these numerical simulations, the diffusion of a single Ag⁺ ion is considered. This choice was motivated by experimental measurements which have shown that the amount per time unit of Ag⁺ ions detached from a layer of AgNPs embedded in silica is extremely low¹³, suggesting that Ag⁺-Ag⁺ interactions in the silica layer are negligible.

The paper is organized as follows. A first part is dedicated to the description of the methods. Details are given concerning the development of the NNP (database generation, NN training, and validity tests) and the molecular dynamics simulations, and the tools used to examine the diffusion processes are introduced. In a second part the nature of the Ag⁺ diffusion is analyzed and discussed. Focus is put on the effect of the presence of defects in the silica matrix and a comparison with anomalous diffusion models is presented. Finally a conclusion is drawn.

## II. Methods

In this section, we first present the methodology used to build the NNP for the Ag⁺/SiO₂ interactions (Sec. II A). Then we give the main parameters used to generate the glass samples and the diffusion trajectories of Ag⁺ inside the silica matrix (Sec. II B). The method and tools used to analyze the diffusion mechanisms are given in Sec. II C.

### A. Machine-learning potential

In order to generate a machine learning potential for the interactions between Si, O and Ag⁺, we used the N2P2 software²⁴,²⁵. It is based on the training, on a database made of DFT energies and forces of atomic configurations, of a NN with Behler-Parrinello (BP) atomic descriptors²⁶. This potential being intended to model the diffusion of a single silver ion, the Ag⁺-Ag⁺ interaction was not considered.

#### 1. Database generation: atomic configurations

To be as accurate and transferable as possible, the NNP must be trained on energies and forces computed on a broad sampling of atomic configurations. In our previous

work²⁷, we have developed a NNP for silica using configurations extracted from DFT-based *ab initio* molecular dynamics (AIMD) simulations. Especially, we have shown that a few hundred configurations from the high-temperature liquids were sufficient to train an accurate potential for amorphous silica.

Unfortunately, the same strategy can not be used to develop a NNP for Ag⁺-O and Ag⁺-Si interactions. Indeed, since there is only 1 Ag⁺ ion in the simulation box, a DFT-based *ab initio* molecular dynamics (AIMD) simulation of Ag⁺ diffusion in silica, can not produce enough configurations to feed the training of Ag⁺-(O/Si) interactions. To overcome this problem, we first made a preliminary Ag⁺-O and Ag⁺-Si simple model in terms of Lennard-Jones-type potentials (LJ) whose parameters were fitted on short AIMD trajectories of an Ag⁺ ion in a 384-atom silica glass at 500K (NVE, cubic box with 17.96 Å edge length). The parameters of these AIMD calculations are given in Sec. II A 2 and the time step is 1 fs. Here although the charge of the Ag⁺ ion is not explicitly taken into account, it is included effectively in the LJ interactions.

Mixing the NNP for SiO₂²⁷ and this LJ potential for Ag-O and Ag-Si, molecular dynamics (MD) simulations were performed using the LAMMPS code²⁸ on a system made of 128 SiO₂ and a single Ag⁺ ion (NVT ensemble, time step of 1 fs, Nosé-Hoover thermostat with damping parameter set to 100 fs). MD simulations at 2 different temperatures have been used to generate the database configurations. At 3600K, a single trajectory of 500 ps was simulated in the liquid state from which 125 configurations were extracted (1 configuration every 4 ps). Then, in order to generate glass samples, the system was cooled down to 300K using two different quench rates (1.6 10¹¹ K/s and 3.2 10¹¹ K/s). In each glass sample, 5 different trajectories of Ag⁺ at 500K were generated using 5 different initial positions, making a total of 10 trajectories of 500 ps each. These simulations evidenced that the Ag⁺ ion moves by jumps: it oscillates around a given position for a while, before jumping to another position, and so on throughout the trajectory. To avoid having redundant configurations in the database, and to better sample the energy barriers, snapshots were not extracted at regular intervals along the trajectory but mainly at the Ag⁺ ion jumps. Thus, the jumps were identified and configurations were extracted every 1.2 ps during each jump, and every 30 ps elsewhere, leading to 1128 configurations. In total, the database contains 1253 configurations.

#### 2. Database generation: energies and forces

The energies and forces of the 1253 extracted configurations were then computed at the DFT level. The DFT calculations were carried out using VASP²⁹⁻³¹. For these calculations, the wavefunctions were expanded into plane-waves up to a cutoff energy of 450 eV at the Γ-point of the Brillouin zone and the PBE exchange and correlation functional was used together with PAW pseudopotentials for Si, O and Ag. The positive charge of the Ag⁺ ion was compensated by an uniform charged background.
### 3. Training

The parameters of the N2P2 NN have been optimized in order to obtain the best root mean square error (RMSE) for both forces and energies. The NN consists of two layers of 20 neurons each, the activation function is a hyperbolic tangent, and we used 71 BP radial functions and 101 BP angular functions as atomic descriptors. Among this set of descriptors, the radial and angular functions involving only Si and O atoms are the same as those used for the NN developed for silica in our previous work²⁷, with a cutoff of 6 Å. The other descriptors, involving the Ag⁺ ion, are defined with a cutoff of 9 Å. For the radial functions, the values of η and $r_s$ were set to {0.001, 0.01, 0.03, 0.05, 0.06, 0.1, 0.15, 0.3, 0.6, 1.5, 2} Å⁻² and {0, 2, 4, 5, 6} Å, respectively, and for the angular functions, the values of η, λ and ζ were set to {0.001, 0.003, 0.007, 0.01, 0.015, 0.05, 0.1, 2} Å⁻², {-1, 1} and {1, 2}, respectively²⁶. Each combination of these parameters can form a descriptor: a large number are tested, and those that have a negligible effect on the training performance are removed.

Training was performed on all energies and on 4.33% of the forces. 90% of database configurations were used for the training, and tests were done on the remaining 10%. After an initial training, the NNP was tested by running a MD simulation on silica liquid at 3000K. A hundred configurations, whose forces and energies were poorly reproduced by the NNP, were extracted and recalculated using DFT, before being added to the training set. 5 training runs were performed with different initial seeds, yielding highly similar results. After the final training, the values of the lowest RMSE for the energies and the forces are respectively 1.61 meV/atom and 0.13 eV/Å for the training and 1.28 meV/atom and 0.12 eV/Å for the test.

### 4. Validity tests

In order to test the accuracy of the NNP, two trajectories of a single Ag⁺ ion in SiO₂ have been produced at 500K (1 ns, amorphous silica) and 3000K (500 ps, liquid silica). 50 configurations were extracted from the trajectory at 500K (mainly at the Ag⁺ ion jumps, as described in Sec. II A 1) and 10 configurations were extracted from the liquid trajectory at 3000K. For each extracted configuration, forces and energies were recalculated using DFT and compared to the ones obtained using the NNP (Fig. 1). The agreement between the NNP and the DFT is very good, both for energies and forces. For energies (Fig. 1(a)), both liquid (high energies) and glass (low energies) configurations are well reproduced.

Finally, in order to assess the reliability of the structural properties of the resulting system, the Ag-Si and Ag-O pair correlation functions (PCF) were computed using configurations extracted (every 0.01 ps) from the MD trajectory (1 ns) at 500K. They are compared to the ones calculated using configurations extracted from the AIMD trajectories used to develop the LJ Ag-Si and Ag-O potential in Sec. II A 1 (see Fig. S1 in Supp. Mat.). The two PCFs agree well, especially for the first peaks.

<figure><img src="image_2.png" /><figcaption>Figure 1. Comparison of energies (a) and forces (b) between DFT and NNP for selected snapshots along 500K and 3000K trajectories.</figcaption></figure>

Therefore, the NNP developed for Ag/SiO₂ reproduces very well this system with a degree of precision comparable to the DFT level. It is worth pointing out that this agreement is very good, even though the ion charge is not explicitly accounted for in the NNP. The NN great flexibility enables us to accurately reproduce Ag⁺-Si and Ag⁺-O interactions and even chemical reactivity in the system, as we will show hereafter.

### B. Molecular dynamics

MD simulations were carried out using the LAMMPS code in the NVT ensemble (Nosé-Hoover thermostat with damping parameter of 100 fs) with a time step of 1 fs.

#### 1. Glass samples generation

In order to study the diffusion of Ag⁺ ion inside the silica matrix, silica glass samples were generated using the NNP by cooling down a liquid silica sufficiently rapidly to avoid crystallization. Fifteen silica glass samples of 6000 atoms each in a cubic box of 44.70 Å edge length have been generated using three different quench rates (5 samples per quench rate): 3.8 10¹¹ K/s (Q1), 3.8 10¹³ K/s (Q2), 3.8 10¹⁴ K/s (Q3). These samples differ by the number of coordination defects. Indeed, it has been shown that the faster the quench rate, the higher the concentration of defects in the silica glass³²,³³. Defects are defined from the coordination of Si (resp. O) atoms (using a cutoff radius $r = 2$ Å). Tab. I shows the average number of over- or under-coordinated Si (Si₅, Si₃) and O (O₃, O₁) atoms, as well as
<figure><img src="image_3.png" /><figcaption>Figure 2. Ring size distribution for glass samples generated with different quench rates.</figcaption></figure>

tetrahedra connected by their edge instead of their corner (edge-sharing tetrahedra, EST), in the glass at 300K, for the different quench rates.

<table><caption>Table I. Average number of structural defects in the 6000-atom glass samples at 300K for each quench rate.</caption><thead><tr><th>Name</th><th>Quench rate [K/s]</th><th>Si<sub>3</sub></th><th>Si<sub>5</sub></th><th>O<sub>1</sub></th><th>O<sub>3</sub></th><th>EST</th></tr></thead><tbody><tr><td>Q1</td><td>3.8 &times; 10<sup>11</sup></td><td>0</td><td>1.6</td><td>0</td><td>1.6</td><td>9.8</td></tr><tr><td>Q2</td><td>3.8 &times; 10<sup>13</sup></td><td>0</td><td>10.4</td><td>0</td><td>10.4</td><td>39.8</td></tr><tr><td>Q3</td><td>3.8 &times; 10<sup>14</sup></td><td>0.8</td><td>18.2</td><td>0.2</td><td>17.6</td><td>59.8</td></tr></tbody></table>

The structure of the different glass samples has been analyzed in terms of pair correlation functions and bond angle distributions (Figs. S2 and S3 in Supp. Mat) and are in good agreement with the ones from the literature, validating the NNP for silica²⁷.

The rings size distributions were computed using the OVITO package³⁴ and are shown in Fig. 2 for the different quench rates. The distributions are centered around sizes of 6- and 7-membered rings for all glasses, as expected in amorphous silica. When the quench rate is increased, we can observe a slight increase of the small (2 and 3) and large (9-11) rings at the expense of the medium-size rings (5-7). The 2-membered rings correspond to EST. This evolution of the ring size distribution with the quench rate is well referenced in the literature²⁷.

## 2. Diffusion of Ag⁺ ions

To study the diffusion properties of Ag⁺ in silica, we have performed 10 MD simulations of 3 ns at 1200K for each glass sample, each simulation being started using different initial position and velocity for the ion.

### C. Diffusion analysis

In this section, the definitions of the variables used to analyze diffusion data are presented. Informed readers can skip this section if they wish, and go directly to the Results (Sec. III).

It is generally admitted that anomalous diffusion can occur when the transport (or diffusion) process no longer obeys the central limit theorem (CLT)³⁵,³⁶.

In the case of normal diffusion, the CLT states that the probability density function P(x,t) of a series of random variables x(t) is a Gaussian if the function P(x,t) (i) is stationary, (ii) has finite variance and (iii) if random variables are independent. Violation of one or more of these conditions can lead to anomalous diffusion. In this case, the mean-square displacement of the random variables (MSD) is no longer proportional to t but:

$$ \langle [x(t) - x(0)]^2 \rangle = K t^{2H} $$

where, for anomalous diffusion, $H \neq 1/2$ and K is a coefficient in units of $[distance^2]/[time^2H]$. In the case of Brownian diffusion, $K = 2dD$ where d is the dimension and D the diffusion coefficient.

In the present case, the random variable x(t) is given by the position $\mathbf{R}(t)$ of the Ag⁺ ion. The square displacement of the ion as a function of time is given by:

$$ \mathrm{MSD}(t) = |\mathbf{R}(t) - \mathbf{R}(t_0)|^2 \quad (1) $$

where $\mathbf{R}(t_0)$ is the initial position of Ag⁺. Even if this quantity is not a "mean" value since only a single ion is considered, the "MSD" designation has been kept hereafter to make reading easier.

The ensemble average of the square displacement is computed by averaging the MSD(t) over the different trajectories:

$$ \begin{aligned} \langle \text{MSD}(t) \rangle &= \left\langle \left| \mathbf{R}(t) - \mathbf{R}(t_0) \right|^2 \right\rangle \\ &= \frac{1}{N_{\text{traj}}} \sum_{i=1}^{N_{\text{traj}}} \left| \mathbf{R}^{(i)}(t) - \mathbf{R}^{(i)}(t_0) \right|^2 \end{aligned} \quad (2) $$

where $\mathbf{R}^{(i)}(t)$ is the position of Ag⁺ at time t in the trajectory (i), and $N_{\text{traj}}$ is the total number of independent trajectories.

It can be also useful to compute the Mean Maximum Excursion (MME), as shown by Tejedor et al.³⁷, in order to better estimate the exponent H. Indeed these authors have shown that the MSD and the MME scale similarly with time in most of the cases. The MME can be computed using:

$$ \mathrm{MME}(t) = \max(|\mathbf{R}(t) - \mathbf{R}(t_0)|^2) \quad (3) $$

As for the MSD, it is possible to average this quantity over the trajectories:

$$ \langle \mathrm{MME}(t) \rangle = \left\langle \max\left(\left|\mathbf{R}(t) - \mathbf{R}(t_0)\right|^2\right) \right\rangle \quad (4) $$

The origin of the anomalous diffusion can be traced back to the violation of one or more of the CLT conditions³⁶. Indeed, it is possible to establish the stationarity of a trajectory by analyzing whether or not the increments of $\mathbf{R}(t)$ (or jumps) are identically distributed along time. This means that the dynamics is statistically invariant at every stage of the diffusion process. If this is not the case, the first CLT condition is violated, resulting in anomalous diffusion due
<table><caption>Table II. Different possible effects leading to anomalous diffusion and the values of the corresponding exponents.</caption><thead><tr><th>Effect</th><th>&lt; 1/2</th><th>1/2</th><th>&gt; 1/2</th></tr></thead><tbody><tr><td>Moses <i>M</i></td><td>slowing down/aging</td><td>no effect</td><td>accelerating</td></tr><tr><td>Noah <i>L</i></td><td>-</td><td>no effect</td><td>rare extreme events</td></tr><tr><td>Joseph <i>J</i></td><td>anti-correlations</td><td>no effect</td><td>positive, long-ranged correlations</td></tr></tbody></table>

to the so-called "Moses" effect, characterized by the exponent *M*. The time average of the velocity vector norm is used to quantify the non-stationarity of the diffusion:

$$ \langle |\mathbf{v}|(t) \rangle = \left\langle \frac{\Delta}{t} \sum_{j=1}^{t/\Delta} |\mathbf{v}(j)| \right\rangle \propto t^{M-1/2} \quad (5) $$

where the velocity vector is defined with respect to the distance traveled by the Ag⁺ ion during an arbitrary time increment Δ, i.e. $\mathbf{v}(j) = \delta\mathbf{R}(j)/\Delta$ with $\delta\mathbf{R}(j) = \mathbf{R}(j\Delta) - \mathbf{R}([j-1]\Delta)$. For normal diffusion, $M = 1/2$ and the process has a stationary distribution of jumps. However, if $M < 1/2$ or $M > 1/2$, the dynamics is decelerating or accelerating, respectively, leading to a non-stationary distribution in time. This can be observed for instance during aging, and results in weak-ergodicity breaking.

The second condition of the CLT can be violated if the jumps variance is not finite. This can happen for instance if the jumps distribution has an infinite tail or grows in time and can lead to anomalous diffusion. It is called the "Noah" effect and is related to the occurrence of extreme events. It can be quantified by the scaling of the square velocity as:

$$ \langle \overline{\mathbf{v}^2(t)} \rangle = \left\langle \frac{\Delta}{t} \sum_{j=1}^{t/\Delta} \mathbf{v}^2(j) \right\rangle \propto t^{2L+2M-2} \quad (6) $$

By definition $L \ge 1/2$. If $L > 1/2$, its value will grow in time even if $M = 1/2$. This happens for instance in the case of Lévy flights$^{38,39}$.

Long range temporal correlations and anti-correlations can also lead to a nonlinear scaling of the mean square displacement. This is called the "Joseph" effect and it can be formally quantified by the scaling of the velocity-velocity auto-correlation function. However it is very difficult to reliably compute this quantity but Vilk at al.$^{36}$ have shown that one can equivalently use the time-averaged MSD instead:

$$ \overline{\text{MSD}(\tau, t)} = \frac{1}{t-\tau} \sum_{j=0}^{t-\tau} |\mathbf{R}(j+\tau) - \mathbf{R}(j)|^2 \propto t^{2L+2M-2} \tau^{2J} \quad (7) $$

where the overline denotes the time average, $\mathbf{R}(j+\tau)$ and $\mathbf{R}(j)$ are the positions of the Ag⁺ ion at times $j+\tau$ and $j$, respectively. It is also possible to average this quantity over independent trajectories, $\langle \overline{\text{MSD}(\tau, t)} \rangle$. For long-range temporal correlations, $J \ne 1/2$, thus violating the third condition of the CLT. Temporally anti-correlated diffusion was reported in the case of fluorescent rhodamine molecules on wet silica surfaces$^{40}$.

Tab. II presents a summary of the different effects that can lead to anomalous diffusion and the values of the cor-

responding exponents. It has been shown$^{36,39,41,42}$ that the above relations yield a fundamental summation relation between the different exponents:

$$ H = J + L + M - 1 $$

This summation relation will be used hereafter to analyze the diffusion process of Ag⁺ in silica.

### III. Results

In this section, the trajectories of the Ag⁺ ion in silica at 1200K are first analyzed using the tools presented in Sec. II C. Then the focus is put on the search for the origin of the observed anomalous diffusion. In a third part, the effect of the concentration of coordination defects in the silica matrix is investigated. Finally, the Ag⁺ diffusion is compared to existing models.

#### A. Nature of the Ag⁺ ion diffusion

Figure 3 presents the square displacement ($\text{MSD}(t)$, Eq. 1) of 4 independent trajectories of 15 ns each, of a Ag⁺ ion in a 384-atom silica glass generated with a quench rate of $1.6 \times 10^{11}$ K/s at 500K, which were obtained using the same initial position but different initial velocities. The evolution of the square displacement is typical of a stochastic process, in which the Ag⁺ ion diffuses by jumps that take place randomly in time. Due to the amorphous structure of the silica matrix, the Ag⁺ ion can be trapped in an residence zone of the silica (cavity or pore) which is surrounded by multiple possible escape paths. These paths are characterized by potential barriers that the ion must cross to escape the residence zone and diffuse.

The jumps can thus be considered as rare events and, in order to obtain sufficient statistics for these latter, we decided to work at higher temperature (1200K) which enables to increase the jump frequency without breaking the silica network. Moreover, the statistics was improved by averaging over (i) 5 different 6000-atom silica glass samples (Q1 samples) and (ii) 10 different initial positions and velocities of the Ag⁺ ion in each glass sample. Overall, 50 independent trajectories of 3 ns were produced.

An illustration of ten independent trajectories of Ag⁺ in one of the 6000-atom silica glass is depicted in Figure 4. For each trajectory, the silver ion moves significantly in the simulation box, validating the choice of a high temperature to increase the jumps statistics. Moreover, the ion seems to move by following a large number of distinct diffusion paths.

An example of square displacement ($\text{MSD}(t)$, Eq. 1) computed on one of these trajectories is shown in Fig. 5(a). As observed previously, the square displacement comprises alternating segments corresponding to periods
<figure><img src="image_4.png" /><figcaption>Figure 3. Square displacements computed on 4 independent trajectories of a Ag⁺ ion as a function of time in a 384-atom silica glass at 500K. The Ag⁺ ion was placed at the same initial position with different initial velocities (v̇<sub>0</sub><sup>1</sup>, v̇<sub>0</sub><sup>2</sup>, v̇<sub>0</sub><sup>3</sup>, v̇<sub>0</sub><sup>4</sup>).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Atomic positions of 10 independent trajectories of Ag⁺ in a 6000-atom silica glass sample, starting from a different initial position. The Si and O atoms are not shown for clarity.</figcaption></figure>

during which the silver ion remains localized in residence zones of the silica matrix, and jumps during which the ion moves through the matrix.

In order to count the number and the size of the jumps and residence times, the Ag⁺ ion position coordinates have been used: a jump is defined when at least one of the Ag⁺ coordinates varies monotonically over a distance greater than a threshold judiciously chosen to distinguish a diffusion from a position fluctuation (see Supp. Mat. for details). Once the jumps were determined, the remaining of the trajectory was attributed to residence times. In Fig. 5(a), the jumps have been colored in magenta on the MSD(t) curve and the residence times have been identified

by cyan horizontal segments. The duration of the jumps is always quite short, whereas their extension (i.e. the distance travelled) exhibits a broad distribution. The residence times are also quite widespread. The corresponding silver ion positions in Fig. 5(b) and (c) are colored following the same color code for jumps (magenta) and residence times (cyan). It can be seen that the residence times are associated to spaces (the residence zones) of different sizes and shapes, ranging from a small spheroid to a much more elongated form with occasional protrusions. Most of the jumps connect the residence spaces, but a few are located within these spaces, indicating that some of the large displacements may take place in the same place as the residence zones.

When averaged over the 50 trajectories, the <MSD(t)> of the Ag⁺ ion appears to be proportional to t<sup>2H</sup> with H = 0.25 (see figure 6(a)), which is characteristic of an anomalous diffusion of the sub-diffusive type. Such anomalous diffusion processes exist in nature<sup>43</sup> and can be observed for instance in situations as diverse as the diffusion of tracer particles in the cytoplasm of mammalian cells, the displacement of motile amoeba, the movement paths of ants, black-winged kite, white stork or eurasian griffon vulture<sup>36</sup>, the diffusion of water molecules on cell membranes<sup>44</sup> or confined in soft environment<sup>45</sup>, or the diffusion of carriers in semiconductors<sup>46</sup>. As the <MSD(t)> curve is noisy, the H exponent was also evaluated on the basis of the evolution with time of the MME averaged over the trajectories <MME(t)> (see figure 6(b) and Eq. 4 in Sec. II C for its definition). The value obtained H = 0.26 is in very good agreement with that based on the <MSD(t)> evolution.

As explained in Sec. II C, it is possible to trace back the processes behind anomalous diffusion by calculating some quantities and analyzing their behavior as a function of time. Following the work of Vilk at al.<sup>36</sup>, we have analyzed the trajectories of the Ag⁺ ions by calculating the velocities (Eq. 5) and the squared velocities (Eq. 6) in order to extract *M* and *L* exponents, allowing to quantify the effects of non-stationarity (Moses effect) and the presence of extreme events (Noah effect) on diffusion, respectively. The results, presented in Figs. 6(c) and (d) for the velocity and the square velocity, respectively, were obtained using a Δ value of 0.005 ns. They show a decrease in time which can be fitted by a power law, in both cases. The fitted exponents allow to extract the *M* and *L* values, equal to 0.44 and 0.51, respectively. *L* close to 0.5 means that extreme events are not the cause of the anomalous diffusion observed here. The *M* value below 0.5 indicates a slight deceleration in diffusion and suggests that the system may be subject to weak-ergodicity breaking<sup>47</sup> due to non-stationarity of the dynamics.

As stated by the authors of Refs. <sup>47</sup> and <sup>48</sup>, in such a case, the ensemble-averaged <MSD(t)> and time-averaged <MSD(τ,t)> do not coincide and the <span style="text-decoration: overline">MSD(τ,t)</span> for each trajectory strongly differs from one trajectory to the other. It has also been shown that the diffusion coefficients obtained from the fits of each MSD(τ,t) are randomly distributed and that their mean value is time-dependent<sup>48</sup>.
<figure><img src="image_6.png" /><figcaption>Figure 5. (a) Square displacement (MSD) of an Ag⁺ ion as a function of time (blue) and its interpretation in terms of jumps (magenta) and residence times (cyan). The inset shows a zoom and the choice of the Δ value (0.005 ns) used to calculate velocities. (b) All positions of the Ag⁺ ion during a single trajectory with jumps in magenta and residence times in cyan. The Si and O atoms are not shown for clarity. (c) Particular zoom of (b) where the Si-O bonds are depicted as two-color lines.</figcaption></figure>

The $\overline{MSD}(\tau,t)$ were calculated using Eq. 7 for each in-dependent trajectory of Ag⁺ and the results are plotted against lag time $\tau$ in Figure 6(e) (green lines). One can observe that most of the $\overline{MSD}(\tau,t)$ align with each other at short $\tau$ but are scattered over a wide range of y-values, in-dicating a large distribution of diffusion coefficients. More-over, a small number of trajectories are flat, indicating that the silver ion has not diffused at all during the 3 ns of MD. Overall, by fitting the ensemble-averaged $<\overline{MSD}(\tau,t)>$ over the short lag times ($\tau < 0.2$ ns), we obtain a J exponent equal to 0.32 (red dashed line in Fig.6(f) and black line in Fig.6(e)). The J exponent is characteristic of temporal correlations, a value different from 0.5 being associated to the Joseph effect. In the present case, $J < 1/2$ indicates the presence of temporal anti-correlations in the dynamics. This type of effect has already been associated with biased movements, diffusion in confined space$^{49}$ or long-range memory.

To summarize, the diffusion of silver ions in amorphous silica exhibits an anomalous sub-diffusive character evidenced by a scaling exponent *H* lower than 0.5. This sub-diffusion originates mainly from temporal anti-correlations (Joseph effect) and a slight deceleration (Moses effect) but is not due to extreme events (Noah effect).

## B. Origin of the anomalous diffusion

In order to understand the structural origin of this anomalous diffusion, we have analyzed the characteristics of the silica matrix along the Ag⁺ ion paths. To this end, we attempted to identify possible correlations between the structure of the silica matrix and the diffusion path of the silver ion, which could explain the stochastic nature of the diffusion. The fact that the ion remains stuck in one place

for a long time is due to its entry into a potential well of varying depth, surrounded by very high energy barriers. Thermal fluctuations allow the ion to cross one of these energy barriers and jump to a nearby site and diffuse until it reaches another residence zone, i.e. another potential well. In the literature, it has been shown that the presence of disorder can induce abnormal diffusion resulting from a random distribution either of wells on the potential energy surface, or of barrier heights to be crossed in order to exit these wells, or of both$^{35,46,50}$. In the present case, due to the amorphous nature of silica, the potential energy surface is very complex and most certainly displays these characteristics.

The local environment of the Ag⁺ ion along the trajectory was characterized, depending on whether it is in a residence zone or in the process of jumping.

Firstly, correlations were sought between the lengths of the residence times and the corresponding silver ion environment, namely (i) the number of its neighbours and their type, (ii) the distances to these neighbours, and (iii) its Voronoï characteristics (cell volume, surface and number of facets) (see Supp. Mat.). The only correlation that was found concerned the type of neighbors around the Ag⁺ ions, in particular the presence of a defect close to the silver ion. In most cases, the defect involved is O1-type (non-bridging oxygen). Indeed, for 72 % of the long residence times ($> 0.2$ ns), a O1 defect was found at a distance lower than 4 Å of the Ag⁺ ion and in 3% of the cases, a EST defect was found in the vicinity of the silver ion.

Given the presence of these defects in the vicinity of the Ag⁺ ion in the majority of long residence times, the question arises as to whether these defects are responsible for the sub-diffusive nature of the diffusion. In order to test
<figure><img src="image_7.png" /></figure>

Figure 6. Ag⁺ diffusion at 1200K in 6000-atom silica samples generated with the Q1 quench rate: (a) <MSD(t)> as a function of time; (b) <MME(t)> as a function of time; (c) <|v|(t)> as a function of time; (d) <v²(t)> as a function of time; (e) Independent time-averaged MSD(τ,t) computed from the 50 Ag⁺ trajectories as a function of lag time τ; (f) <MSD(τ,t)> as a function of lag time τ. Power law fits were used to extract the scaling exponents J, L, M and H (red dashed lines), H was extracted from (a) and (b), M from (c), L from (d) and J from (f). The black line in (e) represents the fit of (f).

this hypothesis, we selected the trajectories of those ions which remain no more than 0.1 ns in the vicinity of a defect, and recalculated the characteristic exponents of diffusion on this subset, named hereafter Q1ND. The obtained values are given in Tab. III, together with the exponents found with the whole set of trajectories. The exponent values obtained on the Q1ND subset differ little from the values calculated for all the trajectories, indicating that the diffusion is sub-diffusive even without the presence of defects close to the ion. Nevertheless, one can notice that the M exponent is slightly higher for the defect-free trajectories and tends towards 0.5, indicating a less pronounced Moses effect.

Secondly, for the Q1ND subset, the size of the rings through which the silver ion jumps was compared with the overall ring size distribution in silica (Fig. 7). To do this, we selected the rings whose barycentre is the closest to the

ion position during the jump. The silver ion passes mainly through large rings (7- to 10-membered rings) while in silica the predominant rings are the 5- to 7-membered ones. This indicates that there may exist a critical ring size allowing the diffusion of the ion.

That being said, we may wonder whether diffusion is slowed down by the Ag-O1 interactions, as the presence of defects seems to be associated with a lower M exponent and long residence times.

### C. Effect of silica defects on the ion diffusion

In order to understand how the presence of defects in silica influences Ag⁺ ion diffusion, we set out to generate glass samples containing higher concentrations of defects (Q2 and Q3). As presented in Sec. II B 1, the obtained glass samples present a large number of defects at 300K, especially of the Si5, O3 and edge-sharing types (see Tab.
<figure><img src="image_8.png" /><figcaption>Figure 7. Probability at 1200K for the Ag⁺ ion to cross a silica ring of a given size during the jumps (for Q1), compared to the ring size distribution of the silica matrix.</figcaption></figure>

I). The defects of the O1 type, however, are not very numerous. By increasing the temperature to 1200K, the silica matrix structure evolves substantially with a modification of the proportion of the different defect types along the trajectory, i.e. the silica matrix is "healing". Indeed, during the temperature scaling from 300K to 1200K, we observe an increase of the numbers of Si3 and O1 defects and a decrease of the Si5 and O3 ones, especially for the Q2 and Q3 samples. Then, at 1200K, the average number of defects in the Q1 samples remains constant but it decreases with time in the Q2 and Q3 samples (Fig. 8). The number of edge-sharing tetrahedra also decreases (not shown). This "healing" of the matrix structure is due to the fact that the quenching rate used to generate the glass samples is very high and leads the system to a metastable state. By increasing the temperature, local reorganizations are possible which lead to the modification of the relative concentrations of the defect types and, eventually, to the decrease of the number of defects. The higher the quench rate, the more important the defect evolution and the matrix "healing". In the Q1 samples, the average number of defects is very low and the matrix healing already takes place during the temperature scaling to reach 1200K. As a consequence, there is no drastic evolution of the number of defects during the trajectories at 1200K (Fig. 8(a)).

If, as we suspect, the defects (especially of the O1 type) contribute to the slowing down of the Ag⁺ ion diffusion, the diffusion characteristics should be modified with increasing the quenching rate used to generate the silica matrix. From the evolution of the different quantities computed on the trajectories obtained with the two higher quench rates Q2 and Q3, the corresponding exponents *H*, *M*, *L* and *J* could be obtained (Figs. S5 and S6 of the Supp. Mat.). The results are presented in Tab. III and in Fig. 9 for the silica samples generated with the three quench rates. One can observe a significant decrease of the *H* exponent extracted from <MSD(*t*)> which decreases from 0.25 for Q1 samples to 0.15 for Q2 samples and to 0.11 for Q3 samples. The decrease is of the same order when the *H* exponent

<figure><img src="image_9.png" /><figcaption>Figure 8. Evolution of the average number of defects in silica as a function of time during the 1200K MD trajectories for the different quench rates. (a) Q1: 3.8 10¹¹ K/s, (b) Q2: 3.8 10¹³ K/s and (c) Q3: 3.8 10¹⁴ K/s.</figcaption></figure>

is extracted from <MME(t)> or obtained by summing the other exponents (J+L+M-1). Note that the three different estimates of the *H* exponent are similar, indicating that the statistics are sufficient to reliably determine the exponents. Only the *L* exponent remains unchanged when increasing the quench rate. The *J* exponent experiences a strong diminution (from 0.32 for Q1 to 0.24 and 0.23 for Q2 and Q3, respectively), as does *M* (from 0.44 for Q1 to 0.39 and 0.38 for Q2 and Q3, respectively). This decrease of the exponents is partly due to a significant increase in the number of trajectories in which the Ag⁺ ion did not diffuse, which is evidenced by a larger number of flat curves among the independent time-averaged MSD(τ,t) (Figs. S5(e) and S6(e)) indicative of a trapped silver ion. Thus, in comparison with the Q1 case, the increase of the average number of defects can be correlated with a decrease in the *M* exponent, indicating greater deceleration, as well as a decrease in the *J* exponent, indicating an increase in anti-correlation phenomena.
<figure><img src="image_10.png" /><figcaption>Figure 9. Evolution of the different exponents as a function of the quench rate used to generate the silica glass samples.</figcaption></figure>

<table><caption>Table III. Values of the different exponents and K coefficients extracted from the fits of the $<MSD(t)> (H (MSD))$, $<MME(t)> (H (MME))$, $<|v|(t)> (M)$, $<v^2(t)> (L)$ and $MSD(\tau,t)$ (J) for the diffusion of Ag⁺ in silica glass samples generated using different quench rates Q1, Q2 and Q3 and for the trajectories subset without defects Q1ND.</caption><thead><tr><th>Quench rate</th><th>Q1<sup>ND</sup></th><th>Q1</th><th>Q2</th><th>Q3</th></tr></thead><tbody><tr><td>M</td><td>0.48</td><td>0.44</td><td>0.39</td><td>0.38</td></tr><tr><td>L</td><td>0.50</td><td>0.51</td><td>0.50</td><td>0.50</td></tr><tr><td>J</td><td>0.31</td><td>0.32</td><td>0.24</td><td>0.23</td></tr><tr><td>H (MSD)</td><td>0.23</td><td>0.25</td><td>0.15</td><td>0.11</td></tr><tr><td>H (MME)</td><td>0.28</td><td>0.26</td><td>0.14</td><td>0.10</td></tr><tr><td>H (sum)</td><td>0.29</td><td>0.27</td><td>0.13</td><td>0.11</td></tr><tr><td>K [Å<sup>2</sup>/ns<sup>2H</sup>]</td><td>320</td><td>336</td><td>176</td><td>99</td></tr></tbody></table>

To analyze the effect of defects on Ag⁺ ion trajectories, the probabilities of the ion remaining for a given residence time and jumping a given distance were estimated for the three quench rates. The results are presented in Fig. 10(a) and (b) for the residence times and the jump distances probabilities, respectively. In the Q1, Q2, and Q3 samples, the residence times show an exponential decay with a long tail towards longer times. As the quench rate is increased, the residence times are longer and the tail is also increased, even reaching the time scale of the simulation. The presence of defects is thus directly correlated to the increase of the time spent by the ion in the residence zones. Correspondingly, the probability of jump lengths is modified as the quench rate is increased. For the Q1 case, the distribution is centered around 4 Å with a tail towards longer distances up to ≈ 12 Å. When the quench rate is increased (Q2 and Q3), this long tail disappears to the benefit of shorter jumps of ≈ 4 Å. To summarize, during diffusion in the silica matrix, the more the Ag⁺ ion encounters defects, the more likely it is not only to jump over shorter distances, but also to stay longer in one place (longer residence times). These two effects are accompanied by a reduction of the diffusion by increasing both the Moses effect (deceleration) and the

Joseph effect (anti-correlation).

In case Q1, we have shown that long residence times are 72 % of the cases associated with the presence of O1-type defects in the vicinity of the Ag⁺ ion. These percentages have also been estimated in the case of silica matrices with a greater number of defects (Q2 and Q3 samples) (Fig. S9 in Supp. Mat.). As shown in Tab. IV, for Q2 and Q3, the increase of the number of O1-type defects in silica correlates with the prevalence of this type of defects in the vicinity of the Ag⁺ ion during long residence times (> 0.2 ns). In around 12-13% of cases, we note the presence of Si5-type defects in the vicinity of the silver ion, as well as the punctual presence of EST-type defects. Furthermore, there are hardly any long residence times without associated defects.

Since the concentration of defects evolves during the dynamics for high quenching rates (matrix "healing", Fig. 8), it would be interesting to estimate the evolution of the percentage of long residence times associated with a defect with time. To this end, this percentage have been computed dividing the trajectories in three parts (0-1 ns, 1-2 ns and 2-3 ns). The results show that the percentage of residence times associated with a Si5 or EST defect does not appear to change significantly over time (Tab. IV). On the other hand, this percentage increases overall for O1-type defects over the course of the simulation, although the number of such defects decreases with time in cases Q2 and Q3. This can be explained by the fact that the probability of Ag⁺ having encountered a defect, and therefore becoming trapped, increases during the course of the simulation, even if the concentration of the average number of defects is lower at the end of the simulation. We can therefore conclude that, in cases Q2 and Q3, the Moses effect is the result of a deceleration that increases with time, as the probability of the Ag⁺ ion being trapped increases over time.

<table><caption>Table IV. Percentage of residence times > 0.2 ns for which a defect is found in the vicinity of the Ag⁺ ion (within a sphere of radius = 4 Å), for the three quench rates and for different parts of the trajectories at 1200K.</caption><thead><tr><th>Quench rate</th><th>Time [ns]</th><th>Defect type</th><th>No defect</th></tr><tr><th></th><th>Si5</th><th>O1</th><th>EST</th></tr></thead><tbody><tr><td rowspan="4">Q1</td><td>0-1</td><td>0</td><td>57</td><td>0</td><td>43</td></tr><tr><td>1-2</td><td>0</td><td>82</td><td>0</td><td>18</td></tr><tr><td>2-3</td><td>0</td><td>80</td><td>7</td><td>13</td></tr><tr><td><b>0-3</b></td><td><b>0</b></td><td><b>72</b></td><td><b>3</b></td><td><b>25</b></td></tr><tr><td rowspan="4">Q2</td><td>0-1</td><td>14</td><td>82</td><td>9</td><td>5</td></tr><tr><td>1-2</td><td>18</td><td>80</td><td>0</td><td>8</td></tr><tr><td>2-3</td><td>14</td><td>89</td><td>0</td><td>3</td></tr><tr><td><b>0-3</b></td><td><b>13</b></td><td><b>83</b></td><td><b>2</b></td><td><b>6</b></td></tr><tr><td rowspan="4">Q3</td><td>0-1</td><td>15</td><td>88</td><td>6</td><td>3</td></tr><tr><td>1-2</td><td>15</td><td>92</td><td>10</td><td>0</td></tr><tr><td>2-3</td><td>13</td><td>97</td><td>8</td><td>0</td></tr><tr><td><b>0-3</b></td><td><b>12</b></td><td><b>93</b></td><td><b>7</b></td><td><b>1</b></td></tr></tbody></table>

To summarize, the O1 defects seem to act as traps for the Ag⁺ ion, first by stopping it and then by keeping it for long times. The O1 atoms being under-coordinated oxygen
<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Residence time probability and (b) jump lengths probability for the Ag⁺ ion at 1200K in the silica matrix samples generated with different quench rates.</figcaption></figure>

atoms, they are negatively charged and are therefore keen to form a chemical bond with the Ag⁺ ion. The formation of these bonds is reflected in the Ag-O pair correlation function (see Fig. 11). Indeed, when moving from Q1 to Q2 samples, one can observe a decrease in amplitude of the first peak (centered around 2.6 Å) accompanied by an appearance of a peak centered around 2.2 Å, signature of the formation of Ag-O bonds. This new peak is even more pronounced in Q3. Once trapped in an Ag-O bond, the ion needs extra energy to escape from the residence zone and diffuse, as this bond needs to be broken. Note that at 1200K, coordination defects are likely to diffuse through the silica matrix.

<figure><img src="image_12.png" /><figcaption>Figure 11. Pair correlation functions between the Ag⁺ ion and the oxygen atoms of the silica samples generated with the different quench rates.</figcaption></figure>

In addition to the O1 defects that the silver ion may encounter on its path, the O1 can be created when an Ag⁺ meets a Si5-type defect. Indeed, in that case, the Si atom being over-coordinated, the Si-O bonds are longer and weaker, allowing the Ag⁺ ion to break one of them,

leading to the formation of the O1 defect and thus a Ag-O bond. At this point, we would like to draw attention to the fact that the machine learning potential we have developed for this study is capable of reproducing this reactivity, which would be impossible to do with a classic empirical potential.

### D. Comparison to anomalous diffusion models

Anomalous diffusion has been the subject of numerous theoretical studies, and several models have been developed to account for this phenomenon³⁵. Among these models, the three most popular are (i) the Continuous-Time Random Walk (CTRW)⁴⁶ in which the particle (walker) is allowed to move after having arrested for a waiting time $t_w$ which is drawn from a power-law distribution, (ii) the Fractional Brownian Motion (FBM)⁵¹ which is a generalization of the Brownian motion characterized by time correlations in the random walks, and (iii) the diffusion on fractals (Fractal Motion, FM)⁵⁰ which is a random walk on a fractal lattice, where the walker meets bottle-necks and dead-ends that are on all scales.

All these models may present a sub diffusive behavior that scales with time as a power law with an exponent lower than one. However, each of these models has its own particular characteristics, specific to the situations they were developed for⁵². For example, the FBM model is ergodic, while the CRTW and FM models are weakly non-ergodic. In other words, the FBM model is stationary ($M=0.5$) which implies that $H=J$ in the case $L=0.5$, whereas this is not the case for the two other models. The degree of deviation from ergodicity can be measured by estimating the quantity:

$$EB = \frac{\langle \overline{MSD(\tau, t)} \rangle}{\langle MSD(\tau) \rangle} \quad (8)$$

which is equal to 1 in the case of normal diffusion as well as in the case of anomalous diffusion when the system is ergodic⁵³.

Moreover, Tejedor et al.³⁷ have shown that the use of
MME can help distinguish between several of these models. For instance, for the FBM model, the H exponent extracted from the fit of the <MME(t)> should be different from the one obtained by fitting the <MSD(t)>, whereas they should be similar for the two other models. In addition, they have defined the following ratio:

$$R_{\text{MSD}} = \frac{\langle \text{MSD}(t)^2 \rangle}{\langle \text{MSD}(t) \rangle^2} \qquad (9)$$

$$R_{\text{MME}} = \frac{\langle \text{MME}(t)^2 \rangle}{\langle \text{MME}(t) \rangle^2} \qquad (10)$$

and determined the limiting values of these ratio for each model in 2 and 3 dimensions. In the case of the FBM model, these ratio are the same as the ones for the Brownian motion, i.e. $R_{\text{MSD}} = 5/3$ and $R_{\text{MME}} = 1.36$. For the CRTW case, these ratio depend on the value of H in dimension 3. Finally, it should be noted that in the case of the CRTW, the J exponent should be equal to 0.5.

By analyzing the diffusion characteristics of Ag⁺ in silica, we may be able to determine whether one of these diffusion models is able to describe the observed processes. For this purpose, the main characteristics of these 3 models are summarized in Tab. V and compared with those of the Ag⁺ diffusion behavior. In the first 3 columns, we reported the expected behavior of the <MSD(t)>, <MME(t)> and <MSD(τ,t)> as a function of time for the 3 models and the corresponding exponents that were extracted from our simulations. A quick look at these exponents for the different cases show that the H exponent fitted on MSD and on MME quantities differ significantly for the trajectories without defects Q1ND, whereas they are quite similar for the other cases, Q1, Q2, and Q3.

<figure><img src="image_13.png" /><figcaption>Figure 12. Evolution of the ergodicity parameter EB for the Ag⁺ diffusion as a function of the lag time τ in the Q1ND subset and in the Q1, Q2 and Q3 samples.</figcaption></figure>

In the 4th column, the ergodic behavior of the model is given. In order to evaluate it more quantitatively, the ergodicity parameter EB has been computed for the 4 considered cases, namely the Q1ND, Q1, Q2, and Q3 cases. The evolution of the EB parameter as a function of time is

presented in Fig. 12. For the Q1ND case, the EB parameter quickly reaches the value of 1 and oscillates around this value. When all trajectories are considered in the Q1 case, the EB parameter is close to 1, but remains slightly lower than 1 until the end of the simulation, when the value of 1 is reached. For cases Q2 and Q3, the EB parameter is very low (< 0.5) at the start of the simulation, then increases linearly to reach the value of 1 at the end of the simulation. With this result together with the values of the M exponent, everything seems to indicate that two different types of diffusion process are involved, depending on whether the silver is likely to encounter a defect on its path or not. Without defects, the ions trajectories are ergodic, whereas the more likely they are to encounter defects, the more they become non ergodic. In the case of the system without defects, the ergodic nature of the diffusion, along with the fact that different H exponents are obtained from the MSD and the MME, suggests that an FBM model could be used to describe the diffusion of silver in a defect-free silica matrix. Besides, for the other cases, the fact that the J exponent is always different from 0.5 rules out the CRTW model.

Let us now analyze the values of the $R_{\text{MSD}}$ and $R_{\text{MME}}$ ratios and compare them to the different models. The evolution of these ratio as a function of time is depicted in Fig. 13 and the values of these ratio at the end of the simulation are reported in the 5th and 6th columns of Tab. V. By comparison with the values expected from the different models, the ratio obtained for the Q1ND subset are quite compatible with those of the FBM model. In the opposite, none of the ratio $R_{\text{MSD}}$ and $R_{\text{MME}}$ computed for the other cases Q1, Q2, and Q3 are compatible with one of the considered models.

<figure><img src="image_14.png" /><figcaption>Figure 13. Evolution of the $R_{\text{MSD}}$ and $R_{\text{MME}}$ ratio for the Ag⁺ diffusion in the Q1ND subset and in the Q1, Q2 and Q3 samples. The expected values for the FBM model and the CRTW model in the case H=0.25 are shown on right side of the graph.</figcaption></figure>

Finally, in order to investigate an eventual fractal nature of the diffusion, the probability for the silver ion to be outside a growing sphere of radius $r_{0t}^H$ at time t was evaluated for the Q1, Q2 and Q3 cases. For this, the square root of
<table><caption>Table V. Main characteristics of the different diffusion models and comparison with the corresponding ones for the Ag⁺ ion in the samples Q1, Q2, and Q3, and in the Q1ND subset.</caption><thead><tr><th>Model &lt;MSD(t)&gt; &lt;MME(t)&gt; &lt;MSD(&tau;,t)&gt;</th><th>Ergodicity</th><th>R<sub>MSD</sub></th><th>R<sub>MME</sub></th><th>Growing sphere</th></tr></thead><tbody><tr><td>CTRW</td><td>&alpha; t<sup>2H</sup></td><td>&alpha; t<sup>2H</sup></td><td>&alpha; &tau;</td><td>no</td><td>&gt;2.62 (2H=0.5)<br/>&gt;3.00 (2H=0.3)<br/>&gt;3.15 (2H=0.22)</td><td>&gt;2.14 (2H=0.5)<br/>&gt;2.45 (2H=0.3)<br/>&gt;2.57 (2H=0.22)</td><td>constant</td></tr><tr><td>FBM</td><td>&alpha; t<sup>2H</sup></td><td>&alpha; t<sup>2H'</sup>, H' &ne; H</td><td>&alpha; &tau;<sup>2J</sup></td><td>yes</td><td>&lt;1.67</td><td>&lt; 1.36</td><td>constant</td></tr><tr><td>FM</td><td>&alpha; t<sup>2H</sup></td><td>&alpha; t<sup>2H</sup></td><td>&alpha; &tau;<sup>2J</sup></td><td>no</td><td>n.e.</td><td>n.e.</td><td>&alpha; t<sup>H(d-d<sub>f</sub>)</sup></td></tr><tr><td>Q1<sup>ND</sup></td><td>&alpha; t<sup>0.46</sup></td><td>&alpha; t<sup>0.56</sup></td><td>&alpha; &tau;<sup>0.62</sup></td><td>yes</td><td>&approx; 1.55</td><td>&approx; 1.33</td><td>constant</td></tr><tr><td>Q1</td><td>&alpha; t<sup>0.50</sup></td><td>&alpha; t<sup>0.52</sup></td><td>&alpha; &tau;<sup>0.64</sup></td><td>yes/no</td><td>&approx; 1.95</td><td>&approx; 1.52</td><td>constant</td></tr><tr><td>Q2</td><td>&alpha; t<sup>0.30</sup></td><td>&alpha; t<sup>0.28</sup></td><td>&alpha; &tau;<sup>0.48</sup></td><td>no</td><td>&approx; 2.5</td><td>&approx; 1.95</td><td>constant</td></tr><tr><td>Q3</td><td>&alpha; t<sup>0.22</sup></td><td>&alpha; t<sup>0.20</sup></td><td>&alpha; &tau;<sup>0.46</sup></td><td>no</td><td>&approx; 2.45</td><td>&approx; 1.89</td><td>constant</td></tr></tbody></table>

the K coefficient estimated from the fits of the $<MSD(t)>$ (Tab. III) were chosen for the $r_0$ parameter, so that the radius is consistent with a distance. As can be seen in Fig. S10 in Supp. Mat., in all cases, the probability oscillates around a constant or decreases slightly in the second half of the simulation but it never scales as a power law, as it should if we were dealing with diffusion on a fractal.

In summary, by comparing the characteristics of the Ag⁺ diffusion in the different silica samples to existing diffusion models, we have identified two possible diffusion regimes. The first one corresponds to the ion diffusion in a defect-free silica matrix and can be related to a FBM-like diffusion. It is stationary ($M \approx 0.5$, $H \approx J$) and ergodic and its anomalous character comes from anti-correlations ($J < 0.5$). The stochastic behavior of the diffusion comes from the disordered nature of the amorphous silica, that creates a very complex potential energy surface for the ion, with a broad distribution of wells, of barriers and paths to escape those wells. The second regime is due to the presence of defects that act as traps for the Ag⁺ ion, implying long residence times when encountered. Consequently, the diffusion becomes non ergodic as non stationarity sets in ($M < 0.5$). As the concentration of defects increases, the more the silver ion is susceptible to be trapped and the more non ergodic and non stationary the diffusion becomes. This second type of diffusion is difficult to characterize, as it is added to the first type, thereby modifying exponents and long-term behaviors. Ultimately, silver ion diffusion can be described as diffusion on a potential energy surface consisting of a random distribution of wells and barriers, complemented by a number of extremely deep potential wells (Ag-O bond traps) randomly distributed in space and characterized by very long residence times.

## IV. Conclusion

The diffusion of an Ag⁺ ion in amorphous silica was simulated using an interaction potential based on a neural network. As the diffusion process is stochastic, statistically significant sampling was required to characterize it. In addition, as bonds break and form during diffusion, a method that takes chemical reactivity into account was essential. Machine learning potentials meet both of these needs, combining the low computational cost of empirical potentials

with the accuracy and chemical reactivity of *ab initio* methods.

Silver diffusion was found to be anomalous and subdiffusive, the origin of which could be attributed to deceleration and temporal anti-correlations. This sub-diffusion can be attributed primarily to the disordered nature of the silica matrix. Furthermore, it is magnified by the presence of coordination defects within the silica matrix. These defects, in particular the under-coordinated oxygen atoms (O1), act as traps for Ag⁺, by forming O-Ag bonds, thereby limiting the jump length and retaining the ion for long duration. By comparison with existing diffusion models, the diffusion mechanism in the absence of defects appears to be of the Fractional Brownian motion type. This mechanism is substantially modified by the presence of defects.

The diffusion of species in amorphous silica has been extensively studied⁵⁴ in particular the diffusion of water⁵⁵⁻⁵⁷ or the self-diffusion of network-modifier cations⁵⁸⁻⁶⁰. In the case of Ag⁺, the situation is significantly different than that observed with modifiers because (i) the ion is isolated, and (ii) the silica network is fully polymerized. The present mechanisms are therefore closer to those that might be observed, for example, in the case of gas diffusion such as Ne or He⁶¹, atomic hydrogen⁶² or other metallic species⁶³. Note that machine learning methods have recently been used to decode the diffusion mechanisms at work in many situations⁶⁴⁻⁶⁶. Notably, an open competition has been organized among the scientific teams to compare different approaches⁶⁷.

It is reasonable to question to what extent the results observed at 1200K are representative of what happens at room temperature. Since the silica network is not substantially modified at 1200K compared to 300K, we can assume that the main diffusion mechanism would be the same, yielding a similar H exponent. At low temperature, the silica matrix is less flexible, and the ion has less thermal energy. Consequently, a slower diffusion is expected, resulting in a lower diffusion coefficient. Regarding the diffusion mechanism involving defects, at lower temperatures, the Ag⁺ ion would have a lower probability of encountering a defect due to its reduced diffusion. However, once the Ag-O bond is formed, it would be more difficult to break due to insufficient energy. One point that remains
to be investigated is the influence of the presence of other chemical species ($O_2$, $H_2O$) or of a higher concentration of $Ag^+$ ions on the diffusion mechanisms.

Finally, to answer the initial question of this study regarding the development of anti-bacterial devices based on silica-coated AgNPs, two possible approaches emerge to tune the $Ag^+$ ions release through the silica layer. First, by controlling the amount of O1-type defects, the $Ag^+$ diffusion can be more or less slowed down. This could be done, for instance, by irradiation of the silica matrix$^{68,69}$. Secondly, diffusion paths could be opened or closed by changing the overall ring size distribution in the matrix, for example through a slight modification of its density. The effect of a change in density on He and Ne diffusion has already been studied and appears to be a promising avenue$^{61}$.

## Conflicts of interest

There are no conflicts to declare.

## Supplementary Material

The Supplementary Material contains (i) The comparison of the Ag-Si and Ag-O pair correlation functions at 500K for one $Ag^+$ ion in a 384-atom silica glass obtained from NNP and DFT MD; (ii) The pair correlation functions and bond angle distributions of a 6000-atom silica glass at 1200K for the three different quench rates Q1, Q2 and Q3; (iii) The protocol used to detect the jumps in the trajectories; (iv) The evolution of the diffusion quantities and their fits for the Q2 and Q3 quench rates; (v) The correlations between the residence times and the local environment of the $Ag^+$ ion. (vi) The probability for the ion to be outside a growing sphere of radius $K^{1/2}t^H$.

## Data availability statement

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

## Acknowledgements

This work was supported by the ANR "BENDIS" (grant no ANR-21-CE09-0008). This work was granted access to the HPC resources of IDRIS and TGCC under the allocations 2022-[A0110910409] and 2024-[AD010915309] made by GENCI and of CALMIP supercomputing center under the allocation p17024. We thank Simona Ispas, Joël Puibasset and Julien Lam for very fruitful discussions.

## Notes and references

<footnote n="1">¹C. I. Crucho, *Applied Materials Today*, 2024, **38**, 102179.</footnote>

<footnote n="2">²A. Rimola, D. Costa, M. Sodupe, J.-F. Lambert and P. Ugliengo, *Chemical Reviews*, 2013, **113**, 4216–4313.</footnote>

<footnote n="3">³X. Li and J. J. Lenhart, *Environmental Science & Technology*, 2012, **46**, 5378–5386.</footnote>

<footnote n="4">⁴J. Danischewski, D. Donelson, M. Farzansyed, E. Jacoski, H. Kato, Q. Lućin and M. Roca, *Langmuir*, 2023, **39**, 1786–1792.</footnote>

<footnote n="5">⁵P. Benzo, C. Bonafos, M. Bayle, R. Carles, L. Cattaneo, C. Farcau, G. Benassayag, B. Pécassou and D. Muller, *Journal of Applied Physics*, 2013, **113**, 193505.</footnote>

<footnote n="6">⁶Y. Hang, A. Wang and N. Wu, *Chem. Soc. Rev.*, 2024, **53**, 2932–2971.</footnote>

<footnote n="7">⁷S. Nuti, A. Fernández-Lodeiro, J. Galhano, E. Oliveira, M. P. Duarte, J. L. Capelo-Martínez, C. Lodeiro and J. Fernández-Lodeiro, *Nanomaterials*, 2024, **14**, 462.</footnote>

<footnote n="8">⁸M. A. Isaacs, L. J. Durndell, A. C. Hilton, L. Olivi, C. M. A. Parlett, K. Wilson and A. F. Lee, *RSC Adv.*, 2017, **7**, 23342–23347.</footnote>

<footnote n="9">⁹G. A. Sotiriou, T. Sannomiya, A. Teleki, F. Krumeich, J. Voros and S. E. Pratsinis, *Advanced Functional Materials*, 2010, **20**, 4250–4257.</footnote>

<footnote n="10">¹⁰H. Granbohm, J. Larismaa, S. Ali, L.-S. Johansson and S.-P. Hannula, *Materials*, 2018, **11**, 80.</footnote>

<footnote n="11">¹¹T.-G. Kim, Y. W. Kim, J. S. Kim and B. Park, *Journal of Materials Research*, 2004, **19**, 1400–1407.</footnote>

<footnote n="12">¹²A. Pugliara, C. Bonafos, R. Carles, B. Despax and K. Makasheva, *Materials Research Express*, 2015, **2**, 065005.</footnote>

<footnote n="13">¹³A. Pugliara, K. Makasheva, B. Despax, M. Bayle, R. Carles, P. Benzo, G. BenAssayag, B. Pécassou, M. C. Sancho, E. Navarro, Y. Echegoyen and C. Bonafos, *Science of The Total Environment*, 2016, **565**, 863–871.</footnote>

<footnote n="14">¹⁴M. Liong, B. France, K. A. Bradley and J. I. Zink, *Advanced Materials*, 2009, **21**, 1684–1689.</footnote>

<footnote n="15">¹⁵B. Despax, C. Saulou, P. Raynaud, L. Datas and M. Mercier-Bonin, *Nanotechnology*, 2011, **22**, 175101.</footnote>

<footnote n="16">¹⁶T. C. Nason, G. Yang, K. Park and T. Lu, *Journal of Applied Physics*, 1991, **70**, 1392–1396.</footnote>

<footnote n="17">¹⁷Z. Amin, N. A. Waly and S. E. Arshad, *Applied Sciences*, 2021, **11**, 5496.</footnote>

<footnote n="18">¹⁸M. M. Ghobashy, M. A. Elkodous, S. H. Shabaka, S. A. Younis, D. M. Alshangiti, M. Madani, S. A. Al-Gahtany, W. F. Elkhatib, A. M. Noreddin, N. Nady and G. S. El-Sayyad, *Nanotechnology Reviews*, 2021, **10**, 954–977.</footnote>

<footnote n="19">¹⁹H. Balout, N. Tarrat, J. Puibasset, S. Ispas, C. Bonafos and M. Benoit, *ACS Applied Nano Materials*, 2019, **2**, 5179–5189.</footnote>

<footnote n="20">²⁰M. Benoit, J. Puibasset, C. Bonafos and N. Tarrat, *Nanoscale*, 2022, **14**, 7280–7291.</footnote>

<footnote n="21">²¹W. Zhang, Y. Yao, N. Sullivan and Y. Chen, *Environmental Science & Technology*, 2011, **45**, 4422–4428.</footnote>

<footnote n="22">²²K. Loza, J. Diendorf, C. Sengstock, L. Ruiz-Gonzalez, J. M. Gonzalez-Calbet, M. Vallet-Regi, M. Koller and M. Epple, *J. Mater. Chem. B*, 2014, **2**, 1634–1643.</footnote>

<footnote n="23">²³Z. Adamczyk, M. Oćwieja, H. Mrowiec, S. Walas and D. Lupa, *Journal of Colloid and Interface Science*, 2016, **469**, 355–364.</footnote>

<footnote n="24">²⁴A. Singraber, T. Morawietz, J. Behler and C. Dellago, *Journal of Chemical Theory and Computation*, 2019, **15**, 3075–3092.</footnote>

<footnote n="25">²⁵A. Singraber, J. Behler and C. Dellago, *J. Chem. Theory Comput.*, 2019, **15**, 1827–1840.</footnote>

<footnote n="26">²⁶J. Behler and M. Parrinello, *Phys. Rev. Lett.*, 2007, **98**, 146401.</footnote>

<footnote n="27">²⁷S. Trillot, J. Lam, S. Ispas, A. K. A. Kandy, M. E. Tuckerman, N. Tarrat and M. Benoit, *Computational Materials Science*, 2024, **236**, 112848.</footnote>

<footnote n="28">²⁸A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpeton.</footnote>

<footnote n="29">²⁹G. Kresse and J. Hafner, *Physical Review B*, 1993, **47**, 558–561.</footnote>

<footnote n="30">³⁰G. Kresse and J. Furthmüller, *Computational Materials Science*, 1996, **6**, 15–50.</footnote>

<footnote n="31">³¹G. Kresse and J. Furthmüller, *Physical Review B*, 1996, **54**, 11169–11186.</footnote>

<footnote n="32">³²K. Vollmayr, W. Kob and K. Binder, *Physical Review B*, 1996, **54**, 15808–15827.</footnote>

<footnote n="33">³³S. Halbert, S. Ispas, C. Raynaud and O. Eisenstein, *New J. Chem.*, 2018, **42**, 1356–1367.</footnote>

<footnote n="34">³⁴A. Stukowski, *Modelling and Simulation in Materials Science and Engineering*, 2009, **18**, 015012.</footnote>

<footnote n="35">³⁵J.-P. Bouchaud and A. Georges, *Physics Reports*, 1990, **195**, 127–293.</footnote>

<footnote n="36">³⁶O. Vilk, E. Aghion, T. Avgar, C. Beta, O. Nagel, A. Sabri, R. Sarfati, D. K. Schwartz, M. Weiss, D. Krapf, R. Nathan, R. Metzler and M. Assaf, *Phys. Rev. Research*, 2022, **4**, 033055.</footnote>

<footnote n="37">³⁷V. Tejedor, O. Bénichou, R. Voituriez, R. Jungmann, F. Simmel, C. Selhuber-Unkel, L. B. Oddershede and R. Metzler, *Biophysical Journal*, 2010, **98**, 1364–1372.</footnote>

<footnote n="38">³⁸V. Zaburdaev, S. Denisov and J. Klafter, *Rev. Mod. Phys.*, 2015, **87**, 483–530.</footnote>

<footnote n="39">³⁹E. Aghion, P. G. Meyer, V. Adlakha, H. Kantz and K. E. Bassler, *New J. Phys.*, 2021, **23**, 023002.</footnote>

<footnote n="40">⁴⁰R. Sarfati and D. K. Schwartz, *ACS Nano*, 2020, **14**, 3041–3047.</footnote>

<footnote n="41">⁴¹L. Chen, K. E. Bassler, J. L. McCauley and G. H. Gunaratne, *Phys. Rev. E*, 2017, **95**, 042141.</footnote>
<footnote n="42">⁴²P. G. Meyer, V. Adlakha, H. Kantz and K. E. Bassler, *New J. Phys.*, 2018, **20**, 113033.</footnote>

<footnote n="43">⁴³J. Klafter and I. M. Sokolov, *Phys. World*, 2005, **18**, 29–32.</footnote>

<footnote n="44">⁴⁴E. Yamamoto, T. Akimoto, M. Yasui and K. Yasuoka, *Sci Rep*, 2014, **4**, 4720.</footnote>

<footnote n="45">⁴⁵S. Hanot, S. Lyonnard and S. Mossa, *Nanoscale*, 2016, **8**, 3314–3325.</footnote>

<footnote n="46">⁴⁶H. Scher and E. W. Montroll, *Phys. Rev. B*, 1975, **12**, 2455–2477.</footnote>

<footnote n="47">⁴⁷S. Burov, J.-H. Jeon, R. Metzler and E. Barkai, *Phys. Chem. Chem. Phys.*, 2011, **13**, 1800.</footnote>

<footnote n="48">⁴⁸A. Lubelski, I. M. Sokolov and J. Klafter, *Phys. Rev. Lett.*, 2008, **100**, 250602.</footnote>

<footnote n="49">⁴⁹M. Magdziarz and J. Klafter, *Phys. Rev. E*, 2010, **82**, 011129.</footnote>

<footnote n="50">⁵⁰D. ben Avraham and S. Havlin, *Diffusion and Reactions in Fractals and Disordered Systems*, Cambridge University Press, 1st edn., 2000.</footnote>

<footnote n="51">⁵¹B. B. Mandelbrot and J. W. Van Ness, *SIAM Rev.*, 1968, **10**, 422–437.</footnote>

<footnote n="52">⁵²R. Metzler, J.-H. Jeon, A. G. Cherstvy and E. Barkai, *Phys. Chem. Chem. Phys.*, 2014, **16**, 24128–24164.</footnote>

<footnote n="53">⁵³A. c. v. Godec and R. Metzler, *Phys. Rev. Lett.*, 2013, **110**, 020603.</footnote>

<footnote n="54">⁵⁴G. N. Greaves, *Mineralogical Magazine*, 2000, **64**, 441–446.</footnote>

<footnote n="55">⁵⁵T. Bakos, S. N. Rashkeev and S. T. Pantelides, *Physical Review Letters*, 2002, **88**, year.</footnote>

<footnote n="56">⁵⁶T. Bakos, S. N. Rashkeev and S. T. Pantelides, *Physical Review B*, 2004, **69**, year.</footnote>

<footnote n="57">⁵⁷M. Tomozawa, K. M. Davis, J. H. Seaman and E. M. Aaldenberg, *J Am Ceram Soc*, 2017, **100**, 4548–4561.</footnote>

<footnote n="58">⁵⁸A. Meyer, J. Horbach, W. Kob, F. Kargl and H. Schober, *Phys. Rev. Lett.*, 2004, **93**, 027801.</footnote>

<footnote n="59">⁵⁹P. Hung, N. Yen, L. Vinh, F. Noritake, N. Ha and L. San, *Journal of Molecular Liquids*, 2020, **316**, 113834.</footnote>

<footnote n="60">⁶⁰N. T. Thao, N. V. Yen and P. T. Lien, *Eur. Phys. J. B*, 2023, **96**, 138.</footnote>

<footnote n="61">⁶¹S. V. Kukhtetskiy, E. V. Fomenko and A. G. Anshits, *Membranes*, 2023, **13**, 754.</footnote>

<footnote n="62">⁶²L. S. Morrissey, D. Pratt, W. Farrell, O. Tucker, S. Nakhla and R. Killen, *Icarus*, 2022, **379**, 114979.</footnote>

<footnote n="63">⁶³R. Büngener, W. Pamler and U. Gösele, *Materials Science in Semiconductor Processing*, 2003, **6**, 43–48.</footnote>

<footnote n="64">⁶⁴P. Kowalek, H. Loch-Olszewska and J. Szwabiński, *Phys. Rev. E*, 2019, **100**, 032410.</footnote>

<footnote n="65">⁶⁵H. Loch-Olszewska and J. Szwabiński, *Entropy*, 2020, **22**, 1436.</footnote>

<footnote n="66">⁶⁶X. Qu, Y. Hu, W. Cai, Y. Xu, H. Ke, G. Zhu and Z. Huang, *Phys. Rev. Research*, 2024, **6**, 013054.</footnote>

<footnote n="67">⁶⁷G. Muñoz-Gil, G. Volpe, M. A. Garcia-March, E. Aghion, A. Argun, C. B. Hong, T. Bland, S. Bo, J. A. Conejero, N. Firbas, Óscar Garibo i Orts, A. Gentili, Z. Huang, J.-H. Jeon, H. Kabbech, Y. Kim, P. Kowalek, D. Krapf, H. Loch-Olszewska, M. A. Lomholt, J.-B. Masson, P. G. Meyer, S. Park, B. Requena, I. Smal, T. Song, J. Szwabiński, S. Thapa, H. Verdier, G. Volpe, A. Widera, M. Lewenstein, R. Metzler and C. Manzo, *Nature Communications*, 2021, **12**, 6253.</footnote>

<footnote n="68">⁶⁸T. Bakos, S. Rashkeev and S. Pantelides, *IEEE Transactions on Nuclear Science*, 2002, **49**, 2713–2717.</footnote>

<footnote n="69">⁶⁹A. Alessi, S. Agnello, G. Buscarino, Y. Pan and R. I. Mashkovtsev, in *Applications of EPR in Radiation Research*, ed. A. Lund and M. Shiotani, Springer International Publishing, Cham, 2014, pp. 255–295.</footnote>