Dynamics of ultrathin gold layers on vitreous silica probed by density functional theory

The structure and properties of Au ultrathin films on hydroxyl-free and hydroxylated silica glass surfaces are investigated using ab initio molecular dynamics simulations. Substantial surface structure dependence of Au agglomeration behavior (solid-state dewetting) is found. On hydroxyl-free surfaces, the Au film virtually undergoes instantaneous agglomeration accompanied by the formation of voids exposing a bare silica glass surface. In contrast, simulated annealing of the Au film on hydroxylated surface models leaves its structure unchanged within the simulation time. This points to a key role of reactive defect sites in the kinetics of solid-state dewetting processes of metals deposited on the glass surface. Such sites are important for initial void nucleation and formation of metal clusters. In addition, our calculations demonstrate the crucial role of the appropriate inclusion of dispersion interactions in density functional theory simulations of metals deposited on glass surfaces. For defective, hydroxyl-free glass surfaces the dispersion correction accounts for 35% of the total adhesion energy. The eﬀect is even more dramatic for hydroxylated glass surfaces, where adhesion energies are almost entirely due to dispersion interactions. The Au adhesion energies of 200 and 160 kJ (mol nm 2 ) (cid:2) 1 calculated for hydroxylated glass surfaces are in good agreement with the experimental data.


Introduction
Ultrathin metal layers and metal clusters deposited on silica and silicate glasses play an important role in a broad variety of technical applications ranging from architecture and automotive glazing to electronic devices, sensor technology and catalysis. [1][2][3][4][5][6][7][8] In particular, continuous layers, particulate deposits and clusters of noble metals have attained key importance due to their high imaginary part of the refractive index and, consequently, high optical reflectivity along with the specific position of their plasmonic resonance frequency. This enables a precise control of spectral reflection and transmission (e.g., for low-emissivity and solar protection windows) as well as other properties such as their catalytic or antimicrobial activity. 9,10 Tailoring of the metal layer structure and thickness (usually between a monolayer and a few nm), controlling the formation of atomic agglomerates, cluster deposits or nanoparticles, is the essential step in all these applications. It has been shown that film homogeneity and noble metal agglomeration behavior very strongly depend on subtle variations in the substrate conditions such as the degree of surface hydroxylation. 11 During the deposition process but also during off-line annealing or long-term operation thin metal films can disintegrate into particles, even upon annealing at temperatures well below the melting point, for example, through solid-state dewetting. [12][13][14][15][16][17][18] It is assumed that this process proceeds in three steps: (i) void initiation, (ii) void growth and (iii) void coalescence and particle formation. 17 It has further been postulated that the initial voids form at grain boundaries, vacancies, impurities or pores. [15][16][17] However, the role of surface defect sites of substrates in the solid-state dewetting process has been largely disregarded in this context. The present hypothesis is that voids, initiated by any nucleation center, grow towards the film surface. The time needed for this process is called incubation time, ranging typically from some minutes to hours, depending on film thickness and annealing temperature. 17 After breaking through the entire film, the voids grow and coalesce to form metal particles. However, Seguini et al. 18 have found that ultrathin Au films (0.5 to 6 nm thickness) deposited at room temperature on amorphous SiO 2 by means of e-beam evaporation spontaneously aggregate, forming nanoparticles without any additional thermal treatment.
Computational studies of noble metal agglomeration on silica surfaces are so far scarce. 19 To our best knowledge ab initio investigations of Au on SiO 2 surfaces were limited to a single Au atom or an Au dimer interacting with various cluster or crystalline models of silica. [20][21][22][23] In the present study the structure and properties of Au ultrathin films on hydroxyl-free and hydroxylated silica glass surfaces are investigated for the first time using ab initio molecular dynamics (MD) simulations.

Computational details
2.1. General simulation procedure Fig. 1 depicts the simulation procedure applied in this work. As the first step models of bulk silica glass were generated through simulated melting and quenching at the ReaxFF force field 24,25 level, starting from b-cristobalite supercells. The reliability of the models was assessed by comparison of the calculated structure and properties with the experimental data for silica glass. In the second step, slab models of hydroxyl-free glass surfaces were cut out of the bulk glass models, annealed using the ReaxFF force field and subsequently optimized at the density functional theory (DFT) level. From this, hydroxylated surface models were obtained by saturating all surface defects with oxygen and hydrogen atoms. Finally, a gold monolayer was deposited on each such surface model and MD simulations at the DFT level were carried out. The final equilibrated structures were optimized at the DFT level and used to evaluate gold-glass adhesion energies.

Methods
Structural optimizations and MD simulations at the ReaxFF force field 24 level were carried out using parameters reported in ref. 25. For all ReaxFF MD simulations, an NVT ensemble with a 0.5 fs time step and linear velocity scaling for cooling procedures was used, employing the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) program. 26 Structural optimizations and MD simulations at the DFT level were performed using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 27 and the projector-augmented plane-wave (PAW) 28,29 method as implemented in the VASP (Vienna ab initio Simulation Package) program. [30][31][32][33] Unless stated otherwise, a 250 eV cutoff for the plane wave basis set and the G point for the integration 34 of the first Brillouin zone were used. All MD simulations at the DFT level were performed with a 1 fs time step and linear velocity scaling applied for annealing procedures. For surface models the Grimme dispersion correction (DFT-D3) 35 was added. It has been shown that the DFT-D3 method achieves accuracy comparable with CCSD(T). For Au clusters DFT-D3 yields relative energies that are within 8 kJ mol À1 of the CCSD(T) reference values. 35

Bulk glass models
As noted above, bulk glass models were prepared by a simulated melting and quenching process at the ReaxFF force field level. First, supercells of b-cristobalite were constructed as summarized in Table 1. Using MD simulations, each supercell was equilibrated at 8000 K for 175 ps. At every 50 ps, 75 ps, 100 ps, 125 ps, 150 ps, and 175 ps a structure snapshot was taken. The snapshots were cooled down to 300 K within 2 ns, and optimized. The most stable structure for each supercell size, denoted as G2-G8 (cf. Table 1), was used for further calculations. G2 and G3 were subsequently optimized at the DFT level.
Harmonic vibrational frequencies for G2 were calculated using the density functional perturbation theory 36 as implemented in VASP. Infrared intensities were obtained from the matrix of Born effective charges employing the formula given by Baroni et al. 37 The vibrational frequencies were scaled 38 by 0.9825. This scaling factor was determined by fitting calculated infrared (IR) frequencies to the experimental data. 39 The stiffness tensor C for G2 was calculated at the DFT level using an energy cutoff of 400 eV for the plane wave basis set. Young's modulus E, shear modulus G and the Poisson ratio n were determined under the conditions of an isotropic material, which leads to the following expression for C: ; (1) with The values of E and n were then expressed accordingly as

Pristine glass surface models
Two slab models of pristine glass surfaces (denoted as S1 and S2) were constructed by cutting two-dimensional, about 9 Å thick slabs out of the G2 model. These slabs were optimized, annealed at 1500 K, cooled down to 300 K within 2 ns and subsequently optimized using the ReaxFF force field. Three-dimensional slab models for DFT calculations were prepared by adding a vacuum gap of 15 Å followed by structural optimization at the DFT level.

Hydroxylated glass surface models
Hydroxylated surface models were prepared by saturating all surface defects (threefold coordinated silicon (Si 3 ), non-bridging oxygen (NBO) and three-membered Si-O rings (3-ring)) of S1 and S2 with oxygen and hydrogen atoms. These hydroxylated slab models hS1 and hS2 were optimized, annealed at 300 K within 2 ns and subsequently optimized using the ReaxFF force field. A vacuum gap of 15 Å was added and the models were optimized at the DFT level.

Gold monolayer models
An fcc gold unit cell was optimized at the DFT level using a 10 Â 10 Â 10 Monkhorst-Pack grid for the integrations of the first Brillouin zone and a 400 eV cutoff for the plane wave basis set.
To deposit an Au film a single (111) 14.5 Â 15 Å Au monolayer containing 30 atoms was cut out of the bulk structure. It was compressed to 14.2 Â 14.2 Å to fit the dimension of the glass surface models and placed at a distance of approximately 3.5 Å over the surface. The structures were subsequently optimized at the DFT level keeping the cell parameters constant. The glass surface atoms within the upper 5 Å were allowed to relax with the remaining atoms fixed to their original positions so as to mimic bulk glass behavior. Next, the systems were annealed from 1400 K to 300 K within 5 ps using MD calculations at the DFT level and were subsequently optimized.
The adhesion energy, E ad , of the gold film was calculated as where E Au/SiO 2 is the total energy of the annealed and optimized gold/glass interface model. E Au and E SiO 2 refer to the total energies of the gold layer and of the clean substrate, respectively, and x is the number of gold atoms per supercell (for E ad in kJ (mol Au) À1 ) or the surface area (for E ad in kJ (mol nm 2 ) À1 ). With this definition, adhesion is energetically favorable if the value of E ad is positive. The dispersion energy contribution E dis of the interface was calculated applying eqn (4) and substituting Grimme correction contributions for the total energies.

Glass models
Simulations of amorphous materials usually employ semiamorphous periodic models. Within this approximation a proper description of the amorphous nature requires sufficiently large simulation cells. Therefore, a number of differently sized cells were generated to identify the limit of percolation at which the ensemble size is sufficiently large to describe the properties of bulk silica glass. Fig. 2 shows two examples of the glass models summarized in Table 1. The differences between the energy per SiO 2 unit of models with different sizes are very small, within 1À4 kJ (mol SiO 2 ) À1 , for both ReaxFF and DFT methods. This indicates that already the smallest model is suitable to describe vitreous silica, and also validates the accuracy of the ReaxFF force field. For further assessment of the reliability of our model, selected properties were calculated and compared to the experimental data. Fig. 3 shows the comparison of the total correlation function of G2 with the experimental neutron diffraction data for silica glass. 40 The agreement between calculated and experimental peak positions is very good. The first peak at 1.61 Å corresponds to the average distance between the nearest Si and O atoms, the second peak at 2.64 Å to the average distance between the nearest O atoms and the peak at 3.08 Å to the average distance between the nearest Si atoms. Table 2 summarizes the average interatomic distances and average bond angles for both ReaxFF and DFT optimized models. Both methods reproduce the experimental average distances very well, with DFT slightly overestimating and ReaxFF slightly underestimating the distances. The absolute  Fig. 2 Bulk silica glass models G2 (a) and G5 (b) (see Table 1; Si: grey, O: red). The calculated Young's modulus, the shear modulus and the Poisson ratio of G2 are 95 AE 9 GPa, 35 AE 1 GPa and 0.176 AE 0.02, respectively. The calculated values agree well with the experimental data of 73.5 GPa, 31.5 GPa and 0.17, 43 respectively.
In summary, the properties calculated even for the smallest G2 model containing 64 SiO 2 units are in very good agreement with the experimental data available for bulk silica glass. This confirms the validity of our procedure for generating glass models.

Glass surfaces
Surface properties of real glasses vary widely depending on the environmental conditions. In a very moist environment, silicate glass surfaces are typically fully hydroxylated, whereas heattreated (B1200 1C) glass in a vacuum chamber displays a hydroxyl-free surface. 44 Therefore, two limiting cases of silica glass surfaces were considered here: fully hydroxylated and hydroxyl-free surfaces. The corresponding surface models are shown in Fig. 5. Due to the limited size of the models they do not include low concentration defects, such as two membered Si-O rings (concentration 0.2-0.4 nm À2 ). 45 Table 3 compares the calculated and experimental properties of hydroxylated surfaces. The concentration of surface silanol groups for both hS1 and hS2 models of 4.4 and 4.9 nm À2 , respectively, is in very good agreement with the experimental value 44 of 4.9 AE 0.05 nm À2 . The slight deviation of 0.5 nm À2 for hS1 is due to the relatively small surface area of the models, hence, the removal of one silanol group results in a concentration error of 0.5 nm À2 . However, the ratio of geminal silanols to the total number of silanol groups of 0.22 and 0.1 for both hS1 and hS2, respectively, is in very good agreement with the range of 0.08-0.24, which was observed in experiments. 46,47 Table 4 shows calculated concentrations of NBO and Si 3 surface sites as well as 3-membered Si-O rings for hydroxyl-free glass surface models S1 and S2. NBOs are only present in S1, whereas Si 3 and 3-membered Si-O rings are included in both models. There are only few qualitative 45,48,49 and quantitative [50][51][52] experimental data concerning the defect concentration and structure of hydroxyl-free silica surfaces. Bobyshev and Radzig 50,51 determined the concentration of Si 3 and NBO centers on activated silica surfaces as 0.02 and 0.02 to 0.04 nm À2 , respectively. McCrate et al. 52 found a slightly higher value of approximately 10 À1 nm À2 for Si 3 centers. The limited size of our surface models prevents achieving such low defect concentrations. Already one defect site per surface unit cell results in a concentration of 0.49 nm À2 and modeling even lower ones would require prohibitively large simulation cells.

Gold monolayer on glass surfaces
Simulated annealing of gold monolayers on hydroxyl-free surface models S1 and S2 leads to its virtually instantaneous agglomeration and formation of voids exposing a bare glass surface as demonstrated for the final structures in Fig. 6a and b. The analysis of MD Fig. 3 Comparison of the total correlation function T(r) of G2 (solid line, see Table 1) with the experimental 40 data (dashed line).  trajectories with snapshots shown in Fig. 7 reveals that voids are initiated by opening the Au layer at the NBO sites with the Si 3 centers located beneath the forming Au clusters. The growing voids gradually expose defect-free parts of the silica surface with the NBO sites remaining in close contact with the edge of the voids and Si 3 remaining beneath Au clusters. In the final, optimized structures ( Fig. 6a and b) the corresponding bond distances are in average 2.17 Å and 2.33 Å for Au-NBO and Au-Si 3 , respectively. These findings indicate that defect sites such as NBO and Si 3 centers on hydroxyl-free glass surfaces play a key role in the solid-state dewetting processes by providing centers for initial void nucleation and cluster formation. In addition, the strong interaction with the defect sites stabilizes the forming Au clusters. Our results are supported by experimental 18 and theoretical 19 data. Seguini et al. 18 found that ultrathin Au films (0.5 to 6 nm thickness) deposited at room temperature on amorphous SiO 2 by means of e-beam evaporation spontaneously aggregate, forming nanoparticles without any additional thermal treatment. In addition, rapid Au agglomeration on hydroxyl-free (0001) a-quartz surfaces has also been found by Kuo et al. 19 using MD simulations employing interatomic potential functions.   Table 4 Calculated concentrations of non-bridging oxygen (NBO, nm À2 ), three-coordinated silicon (Si 3 , nm À2 ) and three-membered Si-O ring (3-ring, nm À2 ) surface defects in hydroxyl-free glass surface models S1 and S2 S1 S2 In contrast to hydroxyl-free surfaces, simulated annealing of a gold monolayer on hydroxylated surface models hS1 and hS2 leaves its structure virtually unchanged within the simulation time (see Fig. 6c and d). This finding is consistent with the observed prolonged incubation times for agglomeration of thin Au layers on hydroxylated silica glass surfaces. 14,17 In light of our results, it appears that this behavior can be attributed to the absence of nucleation centers in the form of surface defects, as they are saturated on hydroxylated surfaces. Table 5 shows calculated adhesion energies of the final Au layers on the hydroxyl-free and hydroxylated glass surface models along with the corresponding dispersion contributions. For hydroxyl-free models S1 and S2 the calculated adhesion energies are 40.2 and 33.3 kJ (mol Au) À1 , respectively. The dispersion energy contribution to the adhesion energies of 35% on average is quite significant. In the case of hydroxylated models hS1 and hS2 the adhesion energies of 13.5 and 10.8 kJ (mol Au) À1 , respectively, are substantially lower and consist almost entirely of the dispersion contribution. The calculated Au adhesion energies per area, 200 and 160 kJ (mol nm 2 ) À1 for hS1 and hS2, respectively, are in a good agreement with the experimental value of 216 kJ (mol nm 2 ) À1 measured for gold films on hydroxylated silica glass surfaces. 15 To our best knowledge, the present study represents the first investigation of the behavior of Au layers on realistic models of vitreous silica surfaces using ab initio methods. Previous theoretical investigations were limited to one Au atom or an Au dimer interacting with various models of the silica surface.
However, the reactivity of a single Au atom and small Au n clusters for n o 5 is significantly higher than that of larger clusters or two-dimensional Au layers. [53][54][55][56] In line with our findings for the hydroxylated glass surface, the DFT calculations of Del Vitto et al. 20 using cluster models of amorphous silica and periodic models of an all-silica edingtonite surface found very weak binding energies of less than 10 kJ (mol Au) À1 for an isolated Au atom and an Au dimer interacting with bridging Si-O-Si and terminal silanol groups. A similar result has been obtained for the adsorption of an Au atom on defectfree unsupported and supported crystalline silica films. 21,23 In contrast, for models of hydroxyl-free silica surfaces 20,22 and defective crystalline silica films 23 a very strong interaction of a single Au atom with surface defect sites has been predicted, 305-380 kJ (mol Au) À1 for Si 3 and 208-300 kJ (mol Au) À1 for NBO. The interaction energy decreases significantly with increasing number of Au atoms. For Au 2 the reported calculated binding energies are 93-124 kJ (mol Au) À1 for Si 3 and 70-140 kJ (mol Au) À1 for NBO centers. 20,22,23 Due to already mentioned high reactivity of individual Au atoms and of Au clusters these interaction energies are significantly higher than Au adhesion energies derived in the present work (cf. Table 5).

Conclusions
The results of our simulations demonstrate that properties of silica glass can be accurately described employing relatively small semi-amorphous periodic models. Simulation cells containing only 64 SiO 2 units properly reproduce experimental structural parameters, mechanical properties and IR vibrational spectra of bulk silica glass. The validity of the glass models is further confirmed by good agreement between the calculated and the experimental number of silanol and germinal silanol groups for hydroxylated silica glass surface models constructed from the semi-amorphous cells. Ab initio MD simulations of an Au monolayer deposited on hydroxyl-free and hydroxylated glass surface models point to a substantial  Table 5 Calculated adhesion energies E ad and their dispersion contributions E dis (kJ (mol Au) À1 and kJ (mol nm 2 ) À1 in brackets) of gold on hydroxylated (hS1 and hS2) and hydroxyl-free (S1 and S2) glass surfaces

Hydroxylated
Hydroxyl-free hS1 hS2 S1 S2  12.8 (190) surface dependence of Au agglomeration (solid-state dewetting) behavior. On hydroxyl-free surfaces, Au monolayers undergo virtually instantaneous agglomeration, which is accompanied by the formation of voids exposing a bare silica glass surface. In contrast, simulated annealing of Au monolayers on hydroxylated surface models reveals a significantly higher structural stability. This points to a key role of reactive defect sites such as Si 3 and NBO centers in the kinetics of solid-state dewetting processes of metals deposited on the glass surface. Such defect sites are important for initial void nucleation and formation of metal clusters. The same Au agglomeration behavior is observed for two different surface models with substantially different concentrations and locations of surface defects. This is a strong indication that our conclusion concerning the important role of defect sites in solid-state dewetting processes on the glass surface remains valid even in the case of lower concentrations of surface defects observed in experiments. In addition, the present calculations demonstrate the crucial role of appropriate inclusion of dispersion interactions in DFT simulations of metals deposited on the glass surface. For defective, hydroxyl-free glass surfaces the dispersion correction accounts for 35% of the total adhesion energy. This effect is even more dramatic for hydroxylated glass surfaces, where the adhesion energies are almost entirely due to the dispersion interactions. This finding is in line with a recent report by Ruiz Puigdollers et al. 57 for metal clusters on oxide surfaces. Finally, the gold adhesion energies of 200 and 160 kJ (mol nm 2 ) À1 calculated for hydroxylated glass surfaces are in good agreement with the experimental data.