Dissociative Chemisorption of Water onto Silica Surfaces and

J. Phys. Chem. C 2008, 112, 1507-1515

1507

Dissociative Chemisorption of Water onto Silica Surfaces and Formation of Hydronium Ions T. S. Mahadevan and S. H. Garofalini* Interfacial Molecular Science Laboratory, Department of Materials Science and Engineering, Rutgers UniVersity, Piscataway, New Jersey 08855 ReceiVed: August 29, 2007; In Final Form: October 24, 2007

Molecular dynamics (MD) computer simulation of the adsorption of water molecules onto the vitreous silica surface was performed using a new dissociative water potential.58 The simulations showed dissociative chemisorption of water molecules onto the silica surface, forming silanol (SiOH) groups at a concentration consistent with experimental data. Water penetration and silanol formation ∼7-8 Å below the outermost oxygen are observed. Because of the dissociative nature of the water potential, formation of hydronium ions is allowed, and, whereas seldom observed in the simulations of bulk water, hydronium ions are formed during the reactions causing the formation of the silanols. The formation of hydronium ions has also been observed in ab initio calculations of water adsorption onto silica surfaces. The time evolution of the reactions involving hydronium ions in our MD simulations is similar to that observed in first-principles MD calculations. Hydronium ions offer a mechanism by which initially singly coordinated terminal oxygen (Si-O-) receives a H+ ion from a relatively distant chemisorbed H2O molecule via multiple H+ ion transfer, creating two SiOH sites.

Introduction There have been numerous studies of water interactions with silica because of the ubiquitous nature of each of these in both man-made and environmental systems. For instance, silica and silicate systems are important in technologically relevant materials such as optical fibers,1 glass,2 MEMS devices, microelectronics and wafer bonding,3 separation systems,4 and catalysis,5,6 as well as being an abundant mineral with significant geological implications.7 Studies have shown that water can rupture siloxane (Si-O-Si) bonds via dissociative chemisorption, forming silanol (Si-O-H) sites, which can have significantly deleterious effects on materials such as optical fibers under moderate strain, where stress corrosion cracking can occur.8-13 The presence of water in vitreous silica has a profound influence in many of its thermal, mechanical, and optical properties, and experimental techniques for studying these changes can provide information into the macroscopic aspects of the interactions.14 However, the molecular interactions are more difficult to explain with experimental techniques alone and computational methods have become a useful tool to study these interactions. Molecular orbital and electronic density functional theory calculations at various levels of approximation have been applied to calculate the molecular interaction between silicic acid and water15-20 or water interactions with oxide surfaces.21-25 Of course, the main problem with ab initio calculations is the very limited system sizes (hundreds of atoms, depending upon the level of approximation). In a recent study using a combination of quantum mechanical/classical mechanics calculations, Du et al. attest to the need for larger system sizes (what they call an “extended surface”) to obtain more accurate results for the energy barriers to the dissociation of water on strained siloxane bonds.23 They also observed that the hydration energy is highly dependent upon system size.23 Classical molecular dynamics (MD) simulations have been used to study changes in structure

and energy in SiO2-H2O interactions using non-dissociative water potentials or artificially inducing silanol formation.26-28 The use of a dissociative water potential was applied in MD simulations of water adsorption onto vitreous silica surfaces29,30 and sol-gel polymerization.31-34 Whereas that water potential matched only a few of the bulk water properties, the low-energy structures and energies of the H2O-Si(OH)4 clusters were found to be comparable to ab initio calculated data.31 The commonly shown reaction mechanism for the breaking of siloxane bonds and the formation of silanols was similar to that previously proposed35 and shown schematically below:

The siloxane (Si-O-Si) bond in (a) is strained (either by induced stress in the system or the bond is located at a small, two- or three-membered ring, which can exist at the silica surface). In (b), one silicon becomes pentacoordinated, forming a distorted trigonal bipyramid, which has been discussed much in earlier literature regarding SN2 reactions between silica and water based on ab initio calculations36-38 or molecular dynamics simulations of water adsorption at strained siloxane bonds.29,30 The SN2 reaction and the formation of pentacoordinated silicon as a reaction intermediate is also seen in computational studies of the polymerization of silicic acid (H4SiO4) molecules.31,34 Further interaction between the adsorbed oxygen from water with the silicon and the hydrogen with the bridging oxygen causes the rupture of one of the Si-O bonds in the pentacoordinated silicon (returning the silicon to the tetrahedral coordination) and the formation of silanols (SiOH), as drawn in (c) and (d) above. However, the previous MD simulations also showed that more complex behavior arises when other water molecules are present.30

10.1021/jp076936c CCC: $40.75 © 2008 American Chemical Society Published on Web 01/17/2008

1508 J. Phys. Chem. C, Vol. 112, No. 5, 2008 More recent computational studies have shown additional mechanisms for the adsorption of water onto silica surfaces and the formation of silanol sites.21-25,39 For instance, Del Bene et al.39 showed the benefit of a second water molecule participating in the dissociation of the siloxane bond in an F2H4Si2O molecule. The reaction occurs via transfer of a hydrogen ion between the neighboring water molecules, with the hydrogen of one water molecule attaching to the bridging oxygen, while the second water molecule attaches to the silicon and transfers one of its hydrogen ions to the first water molecule. A fivecoordinated silicon forms during the reaction, with a structure that appears to be similar to the trigonal bipyramid observed in previous computational studies,30,38 although the Si-O bond length to the new oxygen from water is longer than observed in previous ab initio calculations.38 The authors use the need for the second water molecule to discount the MichalskeFreiman model, which involves only one water molecule. However, the ab initio calculations used an F2H4Si2O molecule, with hydrogen and fluorine rather than oxygen or OH terminating the other three bonds to each silicon ion. Previous ab initio calculations by Kudo and Gordon showed that terminating silicon with hydrogen versus OH significantly raised the barriers in hydrolysis and condensation reactions.40 Interestingly, ab initio calculations of cluster models showed rupture of the siloxane bond via reaction with a single water molecule in a manner similar to the reaction mechanism (a-d) shown above, even though those calculations also used hydrogen termination of the silicon ions.41 Therefore, depending upon the type of ab initio calculation, both types of reactions are observed. Several recent computational studies of water interactions with silica surfaces have also shown the important role of additional water molecules adjacent to the reaction site on the hydroxylation process.21-24,42 In these studies, the formation of hydronium ions (H3O+) plays an integral role in the reactions at silica surfaces. Modern MD simulations should be able to capture the complexity of the reactions observed in the less empirical calculations. Although previous MD simulations30 predicted reactions and structures that were qualitatively similar to subsequent ab initio calculations, the earlier water potential was not designed to reproduce most bulk water properties. In order for MD simulations to accurately reproduce the details of the reaction mechanisms involving hydronium formation near reaction sites an accurate and dissociative water potential is required. Although there are a large number of water potentials available in the literature,43-56 most describe rigid, or at best flexible, water molecules. The few fully atomistic water potentials that allow for dissociation of the water molecule do not sufficiently match bulk water properties, especially the liquid equation of state (Figure 15 in ref 57). The authors have recently developed a dissociative water potential that matches bulk water structure, diffusion constant, dipole moment, vibrational density of states, enthalpy of formation, and the liquid equation of state from 273 to 373 K.58 Although hydronium ions do not usually form in simulations of bulk water with this potential, the insertion of a hydronium shows H+ ion migration mechanisms involving both Eigen and Zundel complexes, as well as more complicated configurations.58 Similar results were observed in high level ab initio calculations.59 In the current work, the adsorption and reactions of water onto a silica surface are investigated using molecular dynamics with the new dissociative water potential.

Mahadevan and Garofalini Computational Procedure. The dissociative water potential is based on the rigid water potential developed by Guillot and Guissani.49 However, in our multibody interatomic water potential, intramolecular interactions are added so as to allow for dissociation of the molecule. The multibody potential has both two- and three-body terms. The two body term is given by

U2-body ) Uqq + Uqdqd + Uqqd + Uqdq + Urep + Udisp

(1)

where

Uqq(rij) ) Uqdqd(rij) )

()

qiqj rij × erfc rij β

( ) ( ) ( ) ( )

(2)

()

qid qjd rij rij erf ij × erfc rij β 2ξ

(3)

() () ()

rij rij qiqjd erf × erfc ij rij β x2ξ

(4)

qid qj rij rij erf × erfc Uqdq(rij) ) ij rij β x2ξ

(5)

qid qj rij rij erf × erfc ij rij β x2ξ

(6)

Uqqd(rij) )

Uqdq(rij) )

erfc Urep(rij) ) Arepij ξrij Udisp(rij) )

( ) rij

2ξrij rij

(7)

-C6ij

(8)

rij6

The values of the parameters used for the different atom species are given in Table 1. In general qi ) -4qid and the charge (q) on silicon is four times the charge on hydrogen, and oxygen has twice the opposite charge on hydrogen. The Wolf sum is used to handle the inclusion of the long-range terms in the Madelung potential.60 One of the effects of using the Wolf summation is the presence of the complementary error function term in the charge-based parts of the potential function. This does not appear in the repulsion and dispersion terms. In summary, the final effect of using the Wolf summation is to induce a correction term for all of the charge-based terms of the potential giving the final equation as:

U(1)ele )

1

N

∑ j∑ *i

2 i)l

(

×

rij < Rc

qiqjerfc(rij/β) rij

- lim

{

}) )∑

qiqjerfc(rij/β)

rij erfc(Rc/β) 1 + 2Rc βπ1/2

rij f Rc

(

-

N

qi2 (9)

i)1

The two-body term in the potential does not distinguish between molecules and acts between all of the ion pairs. In addition to the two-body term, the multibody potential has a three-body term of the following form

Dissociative Chemisorption of Water onto Silica Surfaces

U3(rij,rjk,θjik) ) λjik exp

[

γij (rij - rij0)

+

γij rik - rik0

]

[cos(θjik) -

cos(θjik0)]2 (10)

The main purpose of having a three body term is to bias (but not fix) the angles in the water molecule at 104° and the SiO2 tetrahedral angle at 109.47°, and it acts only for rij < rij0 and rik < rik0. λ exists only for HOH, SiOH, SiOSi, and OSiO triplets and is zero for all of the other triplets. The effect of the threebody equation on a j-i-k triplet is to increase the potential energy of the system if the angle is different from the above values. The target angle θojik for the HOH triplet is set at 100° and not 104° to account for the combined effect of both the pair term and the three-body term. The values of the parameters for the three-body term are given in part c of Table 1. Whereas the three-body term biases the structure to the tetrahedral angle for silica, it does not prevent non-tetrahedral structures from forming, such as a trigonal (three-coordinated) silicon, or pentacoordinated silicon, which can exist as defects in silica surfaces or bulk, or three-coordinated oxygen (such as in hydronium ions in water). The parameters for the interactions between the ions in silica were designed to match the structure of amorphous silica, as well as the structure and energy of a silicic acid (H4SiO4) molecule interacting with a water molecule. The latter required the addition of the Si-H interactions (the O-O and Si-O potentials were similar for oxygen in silica or oxygen in water, although the presence of the complementary error function in eqs 2-6 in the Wolf sum meant that the effective interactions differed between these different pairs in silica or water). The structure of the silicic acid molecule was also designed to conform to quantum mechanical calculations by choosing an appropriate value for the Si-H repulsion. The H4SiO4-H2O pair was run for 20 000 moves at 10K, followed by a 20 000 move run at 1 K. Three low-energy configurations, shown in Figure 1, were observed. The interaction energy of the lowestenergy state, shown in part c of Figure 1, was found to be 10.7 kcal/mol, in which the water molecule was in a hydrogen bond single acceptor/double donor configuration. The full coulomb calculation is done on molecules or small clusters, not the Wolf summation, which is required for periodic systems. Using quantum mechanical calculations, Pelmenschikov et al. observed a similar single acceptor/double donor configuration as the lowest-energy configuration, with a minimum energy of 11.6 kcal/mol (varying from 8.8 to 11.6 kcal/mol for the same configuration, depending upon the basis set and electron correlation).20 Molecular dynamics simulations were carried out using a fifth-order Nordseick-Gear predictor corrector algorithm using the above potential to model the interactions between various ions. A 6.4 × 6.4 × 4.2 nm3 box of silica glass was first prepared by a melt quench procedure starting at 6000 K(30 000 steps) and cooled to 4000, 3000, 2000 (100 000 steps each), 1000 (40 000 steps), and 298 K (60 000 steps) under NVE (constant number, volume, energy) conditions, although the volume was changed as a function of temperature based on the thermal expansion coefficient of silica. Each time step equaled 1 femtosecond when simulating the silica system. This glass contained 11 664 atoms and was then equilibrated at 298 K in a 40 ps NPT (constant number, pressure, temperature) run with a hydrostatic pressure of 1 atm. The resultant structure was consistent with the pair distribution functions of bulk amorphous silica. From this NPT run, the configuration with a pressure closest to 1 atm was selected and used to make a surface in a

J. Phys. Chem. C, Vol. 112, No. 5, 2008 1509 manner similar to previous simulations.61,62 The surface was made by removing periodic boundaries in the Z dimension, freezing the bottom half of the glass ions and adding 25 Å of vacuum to the top of the glass. Freezing of the bottom of the system was done purely for computational purposes to reduce computational time. The system was relaxed at 298 K in an NVT simulation for 50 000 time steps. The removal of the periodic boundary conditions in the Z dimension created undercoordinated oxygen and silicon at the surface. Previous simulations showed that many of these defects would be removed by elevated temperature annealing with the formation of small (two- and three-membered) rings on the surface61-63 (an n-membered ring contains n silicon tetrahedra connected by bridging oxygen). Whereas three-membered rings are present in bulk silica,64,65 they increase in concentration at silica surfaces.12 The two-membered rings are not present in bulk silica, but do form at silica surfaces.66,67 These small rings preferentially reacted with adsorbed water in subsequent simulations,29 consistent with experimental findings.12,66 However, in the current work, the temperature was too low to allow for additional formation of small rings on the glass surface, so most defects involved undercoordinated oxygen and silicon. For the production runs of simulating the surface interacting with water, 169 water molecules were placed several angstroms above the surface in a square pattern. The 169 water molecules correspond to a density of 4 water molecules per nm2, which would allow for a maximum silanol concentration of 8 silanols per nm2 if all of the molecules reacted with the silica surface. This is in excess of the range of silanols observed experimentally for a variety of silicas, which is 2-6 SiOH/ per nm2.2,68 Therefore, it is anticipated that not all of the water molecules will react with the silica surface. With this as the starting configuration, a molecular dynamics simulation was run for 35 ps (350 000 timesteps, with a time step of 0.1 fs, which is required for the simulations with the hydrogen ions). Whereas most of the following will describe the results of this run, an additional simulation at 998 K was continued from the run at 298 K for 30 ps to observe reactions that might not occur at 298 K. Results Because the emphasis in this work is on observing the reactions between water molecules and the silica surface, the initial structure of the glass surface is depicted in Figure 2 in terms of the coordination numbers of silicon and oxygen in the glass as a function of the Z coordinate. The location of the oxygen indicates the outer surface of the glass. For the purpose of analysis, the entire glass was sectioned into x,y,z slabs of the dimensions 64 × 64 × 2 Å2 along the Z dimension and the number density of silicons and oxygens of specific coordination in the slabs are shown in the Figure 2. In all of the graphs showing the coordination, the units of the y axis are the number density of species in each slab and the x axis is the Z coordinate of the top of the slab. As expected, the defective species, threecoordinated silicon and one-coordinated oxygen, occur at the surface and are caused by the termination of glass via the removal of the periodic boundary in the Z dimension. Similar to previous simulations and experiment, oxygens are the outermost species and are seen as nonbridging oxygen (NBO) connected to only one silicon. Parts a-d of Figure 3 show one of the observed mechanisms where silanols are formed at the undercoordinated species without breakage of any of the siloxane bonds. A water molecule approaches a region of the surface that contains an NBO and a

1510 J. Phys. Chem. C, Vol. 112, No. 5, 2008

Mahadevan and Garofalini

TABLE 1. (a), Parameters of the Two Body Potential; (b), Charges on Species; (c), Three Body Parameters; (d), The A-matrix of ξrOH(T,P) Equation Given in Part a, Columns A(:,0)-A(:,2); (e), The A-Matrix of the ξrOH(T,P) Equation Given in Part a, Columns A(:,3)-A(:,6) a: ξ (Å)

species

Arep (J)

O-H O-O Si-O Si-Si Si-H H-H

2.283 × 4.250 × 10-17 2.67 × 10-16 7.00 × 10-17 5.00 × 10-16 10-16

C6 (J-Å6)

ξr (Å) P)a

24 24 24 24 24 24

f(T, 0.610 0.373 0.640 0.350

4.226 × 10-18 7.00 × 10-18 3.80 × 10-18

b: species/multiple O Si H

qd/e +0.226 -0.452 -0.113

q/e -0.904 +1.808 +0.452 c:

λ (ergs) γRβ (Å) γRδ (Å) rRβ0 (Å) rRδ0 (Å) θ°jik

O-Si-O ) β-R-β 1 × 10-11 2.8

Si-O-Si ) β-R-β 15 × 10-11 2.0

H-O-H ) β-R-β 30 × 10-11 1.3

3.0

2.6

1.6

109.47

109.47

100

Si-O-H ) β-R-δ 21.2 × 10-11 2.0 1.2 2.6 1.5 109.47

d: 0.655726502 3.403472 × 10-4 -4.057853 × 10-8 1.657262 × 10-12

-1.04442689 × 10-2 -3.986929 × 10-6 4.677537 × 10-10 -1.838785 × 10-14

8.31892416 × 10-5 1.742261 × 10-8 -2.007873 × 10-12 7.549619 × 10-17

e: -3.07929142 × -3.364186 × 10-11 3.800411 × 10-15 -1.355453 × 10-19

10-7

a

5.44770929 × 10-10 2.419996 × 10-14 -2.672717 × 10-18 8.939302 × 10-23

-3.73609493 × 10-13 0 0 0

) 3,n ) 5 m n f(T,P)) ξrOH(T,P) ) ∑mm ) 0,n ) 0 AmnP T .

Figure 1. Stable silicic acid-water configurations. Oxygen, gray; silicon, blue; hydrogen, red. (a) and (b) were observed at 10, 50, and 298 K whereas (c) was observed only at 1 and 10 K. The hydrogenbonded O-H distances in all of the above figures is below 2.4 Å.

nearby three-coordinated silicon (arrows in part a of Figure 3). (Because the bonds between the ions are only drawn if the pair of ions is in the image, some ions will look undercoordinated but are not. Additional hydrogen ions are seen, but again, the other ions in the water molecule to which these hydrogens are attached are not in the figure.) The single water molecule gets close enough to a dangling NBO and a three-coordinated silicon (part a of Figure 3) at the surface and gives up a proton to the oxygen (part b of Figure 3), shortly afterward attaching to the silicon (part c of Figure 3). This results in the formation of two silanols (part d of Figure 3). Parts a-c of Figure 4 show another mechanism for silanol formation in which adsorption of the water molecule occurs initially at the undercoordinated silicon site. The relevant atoms are labeled 1-4 in the images. Part a of Figure 4 shows the initial configuration of a water molecule (labeled 1) getting close

Figure 2. Density of the coordinations of the silicon and oxygen in the silica surface prior to exposure to water. Oxygens are the outermost species, consistent with previous experiments and simulations, with undercoordinated oxygen and silicon occurring at the surface. Dashed lines in the inset bound the silica-vacuum region of the system depicted in the density profiles.

to a three-coordinated silicon (labeled 2). (Note, the apparently small Si-O-Si bond angle at the arrow is really due to the 2D image of the 3D system.) In part b of Figure 4, the water molecule itself is captured by the silicon (indicated by the arrow), resulting in the formation of an SiOH2. The SiOH2 is


Figure 3. Silanol (SiOH) formation on the surface at undercoordinated oxygen and silicon (at arrows). Only a small section of the surface is shown here, with bonds to other ions not drawn in the image also not drawn. Light-blue color indicates oxygen in the water molecule; other ion colors given in Figure 1.

short-lived because of the over-coordination of the oxygen, and it soon dissociates a hydrogen ion to a neighboring nonbridging oxygen (labeled 3) on an originally three-coordinated silicon (labeled 4 in part a of Figure 4). The three-coordinated silicon (4) to which the hydroxyl (3) is attached also reacts with another water molecule to form a geminal site (two hydroxyls on the same silicon, creating an Si(OH)2). The mechanism described by Figure 3 occurs frequently in these simulations and may be the result of the high concentration of NBOs at the outer surface (Figure 2) to which incoming water molecules are initially attracted, whereas three-coordinated silicons (mechanism shown in Figure 4) are slightly below these NBO. Whereas the mechanisms described in Figures 3 and 4 depicted the formation of silanols without breaking the siloxane (Si-O-Si) bond, there were also silanols forming by the breaking of the siloxane bond. A common feature in the siloxane-breaking mechanisms was the formation of a fivecoordinated silicon intermediate that subsequently broke from one of the bridging oxygens, returning the silicon to fourcoordination. An example is shown in Figure 5, which was done at 998 K. Combined classical-quantum simulations at the silica surface and on two member chains of silica have also shown such reactions, and the requirement of the formation of a fivecoordinated silicon has been inferred from these studies.21,22 The dissociative chemisorption of water molecules on a silica surface observed in these MD simulations results in the distribution of species as a function of the Z distance (perpendicular to the silica surface) shown in Figure 6. The data from the initial distribution of four- and two-coordinated silicon and oxygen, respectively, are taken from the initial configuration. The data for the final silicon and oxygen coordination and additional species (SiOH and H2O) are averaged over the last 10 000 moves of the 350 000 move run. This difference in data collection causes a slight change in the distribution in the subsurface region because of the additional motion in the final

J. Phys. Chem. C, Vol. 112, No. 5, 2008 1511 averages caused by normal vibrational motion of the ions that would not be present in the single configuration used for the initial distribution. The concentration of correctly coordinated oxygen and silicon has increased at the surface, which is due to the additional bonds formed with hydrogen and oxygen, respectively, from the water molecules. The curve showing twocoordinated oxygen final includes all of the two-coordinated oxygen attached to silicon, either as Si-O-Si or Si-O-H, thus showing the outermost location of the silica surface. The concentration of silanols (SiOH) indicates silanol formation even at a depth of Z ) ∼40 Å, which is ∼7 Å below the outermost oxygen in the glass. Part of this apparent penetration of reactions is due to the atomistic roughness of the vitreous silica surface, which has been previously discussed in the literature.69-71 Also, because of this atomistic roughness, molecular water penetrates to Z ) ∼40 Å, and again ∼7 Å below the outer glass oxygen. Figures 7-9 provide a descriptive view of this penetration of water into the subsurface. Figure 7 shows side views of two sections of the glass surface, in which hydrogen within 4 Å of a silicon are drawn large (and pink), clearly showing penetration of hydrogen in either H2O or SiOH into the subsurface. Part a of Figure 8 shows several H2O entering a relatively large siloxane ring at the outer surface. Part b of Figure 8 shows the same image but includes more of the subsurface silica in the image, showing how the open surface ring closes off to a smaller (five-membered) ring farther below the surface, forming a channel structure into the subsurface. The five-membered ring is much less reactive with H2O than the smaller two- to threemembered rings that can exist at the silica surface, so the penetration of water will be delayed until a siloxane bond rupture occurs in this ring, effectively extending the channel. Rupture of this ring could occur more readily if the silica were under stress, causing a strain in the siloxane bonds in this fivemembered ring, allowing for a more rapid attack by the water molecule. Such a reaction would open the channel farther into the silica, enabling farther penetration of other water molecules into the silica subsurface. This would provide a mechanism for the well-known enhanced diffusion of water into silica under tensile stress and stress corrosion cracking of amorphous silica under low stress. Because the ion sizes are small in Figure 8, the depth of the penetration of the H2O cannot be seen. However, in Figure 9, the oxygen ions are drawn very large for the graphics to provide for an indication of the third dimension that is not discernible in Figure 8. Figure 9 shows a series of images of three water molecules penetrating the glass. In part b of Figure 9, one of the H2O begins to enter the channel. In part c of Figure 9, the other H2O follow the first in the channel. In part d of Figure 9, the arrow points to the first H2O, showing that it is clearly well below the silica surface (also note that other H2O are attracted to this location). The time evolution of the formation of silanol (SiOH) sites is shown in Figure 10. The final silanol concentration is 3.8/ nm2 (this number includes the SiOH and the chemisorbed H2O as Si(OH2)), which is comparable to experimentally observed values68 as well as those obtained by MD simulation techniques.28 A few hydrogens adsorb onto bridging oxygen, forming a Si2OH complex, but are not considered in the count of silanols. In contrast to previous results obtained by this group, there was no significant change in the distribution of the rings in the glass sample. This could be attributed to the fact that there were no excess small two- or three-membered rings at the surface. It has been previously observed that the smallmembered rings involve strained siloxane bonds that are reactive sites and are readily attacked by water.11,12,29



Figure 4. Silanol formation through dissociative chemisorption of the H2O molecule on three-coordinated silicon. H2O (labeled 1 in (a)) attaches to three-coordinated silicon (2) at the arrow in (b). The H2O dissociates, giving a hydrogen ion to the NBO (3) in (c), which is attached to another three-coordinated silicon (4), which reacts with another H2O to give it a second SiOH (last reaction not shown).

Figure 6. Density profile of correctly coordinated silicon and oxygen in the glass and the concentration of the SiOH and H2O species for the initial glass configuration (prior to water exposure) and an average over the last 10 000 moves for the final configuration. The curve for the two-coordinated oxygen final includes all of the oxygen attached to the silicon in Si-O-Si or Si-O-H and shows the outermost location of the silica surface.

Figure 5. Portion of the surface showing siloxane (Si-O-Si) bond rupture. Light-blue sphere is reacting bridging oxygen, other colors as given in Figure 1. Distorted pentacoordinated trigonal bipyramidal silicon seen in (c), which also shows two elongated (reddish) Si-O bonds during the reaction. (d) Shows that the rupture of the original bridging oxygen (light blue) occurs, with the formation of two SiOHs.

Figure 7. Sideviews of two locations along the surface showing penetration of hydrogen (large, pink spheres) as SiOH or H2O into the silica subsurface. Other colors given in Figure 1.

Hydrogen Ion Transfer. It has been shown that vitreous silica surfaces exposed to water enhance proton conduction,72-76 with an increase in conduction as a function of the hydrogen and/or water concentration. It has been conjectured that the increase in conduction occurs by a hopping mechanism involving a silanol site and an adjacent water molecule.72,74,77 Several recent ab initio or combined QM/MD calculations have shown the role of hydronium ions during reactions between water and a silica surface.21-24 Hence, we looked at the details of our simulations to determine if hydronium ions form during the reactions between water molecules and the glass surface. Hydronium ions do not normally form in our simulations of

bulk water,58 so the formation of hydronium ions in the current simulations would be caused only by the presence of the silica surface. Several processes involving the formation of H3O+ ions and the transfer of H+ ions were observed. Parts a-d of Figure 11 show the process of hydroxylation of the silica surface at a threecoordinated silicon (labeled 1 in part a of Figure 11) and a nearby NBO (labeled 2 in part a of Figure 11) via the formation of a transient hydronium ion. In part a of Figure 11, a water molecule adsorbs onto a previously three-coordinated silicon (at 1), with a nearby water molecule present. In part b of Figure 11, a hydrogen ion is transferred from the adsorbed water molecule to the nearby water molecule (arrow). In part c of


J. Phys. Chem. C, Vol. 112, No. 5, 2008 1513

Figure 8. Top view of a small section of the silica surface with several water molecules (oxygen in water shown as light blue, some hydrogen shown as pink, some as red; other colors from Figure 1). (a), View showing the large ring into which the water migrates. (b) Same view as (a), but (b) includes more silica farther below surface, showing the closure of a large opening created by the large ring seen in (a).

Figure 11. Silanol formation via formation of the hydronium ion, H3O+ (shown as yellow in (c)). (a) Non-dissociate chemisorption of H2O onto three-coordinated silicon at 1, followed by a hydrogen ion moving to an adjacent water molecule (at the arrow) in (b). (c) H3O+ ion shown with yellow oxygen. (d) The hydronium ion loses a hydrogen ion to the NBO labeled as 2 in (c) to form the second SiOH (arrow in (d)).

Figure 9. Snapshot of the region seen in Figure 8, but with ion sizes made large for image, providing a better view of the third dimension, the migration of H2O molecules into the channel (c) and the depth of penetration (d), with the arrow pointing to the lowest oxygen from water.

Figure 10. Time evolution of species formed during reactions of water molecules with the silica surface.

Figure 11, the intact hydronium ion moves to the NBO (labeled 2). In part d of Figure 11, the hydronium transfers one of its hydrogen ions to the NBO, forming the second silanol. From parts b to d of Figure 11, the time of the reaction was 130 fs. Using ab initio calculations, Du et al. showed the formation of a single hydronium in the reaction involving two water molecules near a two-membered siloxane ring.21,23 The first water molecule absorbed onto a silicon, forming a pentacoor-

diated silicon (similar to the structure previously discussed in the literature29,30,36-38). One of the hydrogens from this water molecule is transferred to a nearby H2O molecule, forming the hydronium, which rapidly transfers a hydrogen ion to the bridging oxygen in a barrier-free exchange,23 allowing for the rupture and formation of two SiOH. Cheng et al. observed the exchange of two hydrogen ions involving one hydronium ion to take 100 fs to form a silanol on a Si-O-Si molecule surrounded by four water molecules.42 Another multiple hydrogen-exchange process involving three water molecules and a three-coordinated silicon with a NBO took ∼170 fs from the onset of the first hydrogen exchange between water molecules to end with two additional silanols (SiOH) (see their Figure 542). Parts a-h of Figure 12 show the details of a more complicated hydrogen exchange reaction involving three water molecules in the formation of two SiOH. Part a of Figure 12 shows three water molecules over the silica surface with undercoordinated silicon and oxygen (shown by arrows). The relevant oxygens in the figure are colored light blue in the images (three in the water molecules and the fourth as the NBO). Numbers in some images relate to the oxygen and hydrogen labels discussed below in reference to Figure 13. In part b of Figure 12, a water molecule adsorbs onto the three-coordinated silicon (arrow). In part c of Figure 12, one hydrogen from this H2O molecule interacts with an adjacent H2O, seen at the arrow in mid-transit between the two oxygens. Part d of Figure 12 shows the formation of the first H3O ion (arrow), which exchanges a hydrogen to the next H2O molecule (arrow in part e of Figure 12), which forms the next H3O ion (arrow in part f of Figure 12). In part g of Figure 12, a hydrogen ion is split between the H3O ion and the NBO during the final hydrogen transfer. Part h of Figure 12 shows a later configuration, with the arrows pointing to the final two SiOH’s that form via this process. The total time from part b of Figure 12 to the final hydrogen transfer that creates part h of Figure 12 is ∼150 fs. The time evolution of the O-H bonds in the reaction is shown in Figure 13, in which the numbers on the oxygen and hydrogen



Figure 12. Silanol formation occurring with multiple hydrogen ion transfer via H3O+ ion formation. Light blue represents oxygen involved in this process. Other colors are given in Figure 1. (a) Arrows point to three-coordinated silicon and NBO that eventually become the SiOH sites. (b) Chemisorption of H2O onto three-coordinate silicon (arrow). (c) Hydrogen (labeled 1 and arrow) moves from chemisorbed H2O toward free H2O (labeled 2), forming a split O-H-O bond. (d) The arrow points to the H3O ion formed from the configuration in (c). (e) The hydrogen ion (labeled 2 and arrow) from H3O+ ion moves to next H2O (oxygen labeled 3), which forms the next H3O+ ion in (f) at the arrow. In (g), hydrogen (at arrow) split between oxygen in H3O+ ion and NBO. (h) The hydrogen ion attaches to NBO, forming a second SiOH (arrows point to both SiOH). This multiple H3O ion formation lasts ∼150 fs, similar to ab initio calculations of much smaller systems, showing similar transfer mechanisms (text).

of H3O ion formation similar to the ab initio calculations, with similar time evolutions. In most cases where opportunities for exchange of protons existed in our simulations, the lifetime of any individual hydronium ion was less than 100 fs, with an average around 55 fs. However, we did observe some stable hydroniums that lasted for greater than 300 fs, though none of these long-lasting hydroniums were observed after 6 ps. In these current simulations, the energetics of these reactions observed in our calculations are much more difficult to determine than in the ab initio calculations because of the very many atoms moving in our system and the number of other reactions occurring at any one time. Of course, smaller scale simulations could be performed using fewer water molecules at the sites shown in the images above to obtain energies, and this is planned for future work. Figure 13. Time evolution of specific Oi-Hj bond lengths during the reactions shown in Figure 12. i and j of specific oxygens and hydrogens are given in Figure 12. Arrows indicate a nominal start (when H1 attaches to O2, see Figure 12c) and the end of the reactions that form two SiOH’s, as shown in Figure 12.

relate to those in Figure 12. O2-H1 in Figure 13 involves the oxygen and hydrogen ions in the first hydrogen transfer as labeled in part c of Figure 12. The arrows in Figure 13 indicate the nominal start and end of the process. As soon as H1 bonds to O2, another hydrogen on O2 begins to oscillate away from this H3O ion (O2-H2 curve). H2 moves toward O3 as shown by the O3-H2 curve, oscillating in a manner complementary to the O2-H2 distance oscillations. When H2 finally bonds to O3, the final exchange occurs, as H3 leaves O3 (again with some oscillations). This time frame involving an adsorbed H2O molecule and hydrogen transfer and hydronium formation to two additional H2O molecules observed in our simulations is similar to the result obtained in ab initio calculations by Ma et al.24 In their figure 2, they also show a similar reaction, with a reaction time of 130 fs. However, this reaction involves the formation of a geminal site (two SiOH’s on the same silicon), meaning the transfer distance is quite small. Their Figure 7 shows a reaction more consistent with our Figure 12, again including an adsorbing H2O onto a three-coordinated silicon, formation of two intervening H3O ions, and a second SiOH. The time for this reaction was 146 fs. Therefore, our MD simulations show mechanisms

Conclusion Molecular dynamics simulations of the interactions between silica surfaces and water vapor using a new dissociative water potential revealed dissociative chemisorption of the water via mechanisms that included formation of pentacoordinated silicon as reaction intermediates as well as the formation of hydronium ions. The concentration of silanols observed in the simulations, 3.8SiOH/nm2, is consistent with experimental data. Penetration of water and the formation of silanols was observed ∼7 Å below the outer glass surface. This penetration is caused by the atomistic roughness of the glassy silica surface, which is caused by the normal siloxane-bonded ring structure of silica. This roughness and the ring structure allows for the formation of channels several angstroms deep into the subsurface into which water molecules can penetrate. Hydrogen ion transfer from an adsorption site to a nonbridging oxygen occurring through the formation of hydronium ions was also observed. This transfer process and its lifetime, and the lifetime of the hydronium ions, observed in these MD simulations are very similar to those observed in first principles MD calculations. References and Notes (1) Kurkjian, C. R.; Krause, J. T.; Paek, U. C. J. de Phys. 1982, 43, C9.

Dissociative Chemisorption of Water onto Silica Surfaces (2) Iler, R. K. The Chemistry of Silica; John Wiley & Sons: New York, 1979. (3) Tong, Q.-Y.; Gosele, U. Semiconductor Wafer Bonding: Science and Technology; John Wiley & Sons: New York, 1999. (4) Cabrera, K. J. Sep. Sci. 2004, 27, 843. (5) Roark, R. D.; Kohler, S. D.; Ekerdt, J. G.; Kim, D. S.; Wachs, I. E. Catal. Lett. 1992, 16, 77. (6) Chanthateyanonth, R.; Alper, H. J. Mol. Catal. 2003, 201, 23. (7) Chemical Weathering Rates of Silicate Minerals; White, A. F., Brantley, S. L., Eds.; Mineralogical Society of America: Washington, D.C., 1995; Vol. 31. (8) Brinker, C. J.; Bunker, B. C.; Tallant, D. R.; Ward, K. J. “Chemical Reactivity and the Structure of Gels”, 1986. (9) Michalske, T. A.; Bunker, B. C. J. Appl. Phys. 1984, 56, 2686. (10) Brinker, C. J.; Kirkpatrick, R. J.; Tallant, D. R.; Bunker, B. C.; Montez, B. J. Non-Cryst. Solids 1988, 99, 418. (11) Bunker, B. C.; Haaland, D. M.; Michalske, T. A.; Smith, W. L. Surf. Sci. 1989, 222, 95. (12) Brinker, C. J.; Brow, R. K.; Tallant, D. R.; Kirkpatrick, R. J. J. Non-Cryst. Solids 1990, 120, 26. (13) Chuang, I.-S.; Maciel, G. E. J. Am. Chem. Soc. 1996, 118, 401. (14) Doremus, R. H. J. Mater. Res. 1995, 10, 2379. (15) Chakoumakos, B. C.; Gibbs, G. V. J. Phys. Chem. 1986, 90, 996. (16) Sauer, J.; Morgeneyer, C.; Schroder, K.-P. J. Phys. Chem. 1984, 88, 6375. (17) Sauer, J. Chem. Phys. Lett. 1983, 97, 275. (18) Sauer, J. Chem. ReV. 1989, 89, 199. (19) Lasaga, A. C.; Gibbs, G. V. Am. J. Sci. 1990, 290, 263. (20) Pelmenschikov, A. G.; Morosi, G.; Gamba, A. J. Phys. Chem. 1992, 96, 7422. (21) Du, M.-H.; Kolchin, A.; Cheng, H.-P. J. Chem. Phys. 2004, 120, 1044. (22) Cao, C.; He, Y.; Torras, J.; Deumens, E.; Trickey, S. B.; Cheng, H.-P. J. Chem. Phys. 2007, 126, 211101. (23) Du, M.-H.; Kolchin, A.; Cheng, H.-P. J. Chem. Phys. 2003, 119, 6418. (24) Ma, Y.; Foster, A. S.; Nieminen, R. M. J. Chem. Phys. 2005, 122, 144709. (25) Mischler, C.; Horbach, J.; Kob, W.; Binder, K. J. Phys.: Condens. Matter 2005, 17, 4005. (26) Demontis, P.; Stara, G.; Suffritti, G. B. J. Phys. Chem. 2003, B 107, 4426. (27) Puibasset, J.; Pellenq, R. J. M. Phys. Chem. Chem. Phys. 2004, 6, 1933. (28) Du, J. C.; Cormack, A. N. J. Am. Ceram. Soc. 2005, 88. (29) Feuston, B. P.; Garofalini, S. H. J. Appl. Phys. 1990, 68, 4830. (30) Garofalini, S. H. J. Non-Cryst. Solids 1990, 120, 1. (31) Feuston, B. F.; Garofalini, S. H. J. Phys. Chem. 1990, 94, 5351. (32) Garofalini, S. H.; Martin, G. J. Phys. Chem. 1994, 98, 1311. (33) Martin, G.; Garofalini, S. H. J. Non-Cryst. Solids 1994, 171, 68. (34) Rao, N. Z.; Gelb, L. D. J. Phys. Chem. 2004, B 108, 12418. (35) Michalske, T.; Freiman, S. Nature 1982, 295, 511. (36) Lasaga, A. C. Atomic Treatment of Mineral-Water Surface Reactions. In Mineral-Water Interface Geochemistry; Hochella, M. F., White, A. F., Eds.; Mineralogical Society of America: Washington, D. C., 1990; Vol. 23; p 17. (37) Geisinger, K. L.; Gibbs, G. V.; Navrotsky, A. Phys. Chem. Minerals 1985, 11, 266. (38) Lasaga, A. C. Fundamental Approaches in Describing Mineral Dissolution and Precipitation Rates. In Chemical Weathering Rates of Silicate Minerals; White, A. F., Brantley, S. L., Eds.; Mineralogical Society of America: Washington, D. C., 1995; Vol. 31; p 23. (39) Del Bene, J. E.; Runge, K.; Bartlett, R. J. Comput. Mater. Sci. 2003, 27, 102.

J. Phys. Chem. C, Vol. 112, No. 5, 2008 1515 (40) Kudo, T.; Gordon, M. S. J. Am. Chem. Soc. 1998, 120, 11432. (41) Walsh, T. R.; Wilson, M.; Sutton, A. P. J. Chem. Phys. 2000, 113, 9191. (42) Cheng, H.-P.; Barnett, R. N.; Landman, U. J. Chem. Phys. 2002, 116, 9300. (43) Stillinger, F.; David, C. J. Chem. Phys. 1978, 64, 1473. (44) Rahman, A. J. Chem. Phys. 1972, 33, 3336. (45) Rahman, A.; Stillinger, F.; Lemberg, H. J. Chem. Phys. 1975, 63, 5223. (46) Jorgensen, W. J. Chem. Phys. 1982, 77, 4156. (47) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (48) Robinson, G. W.; Zhu, S.-B.; Singh, S.; Evans, M. W. Water in Biology, Chemistry, and Physics; World Scientific: Singapore, 1996. (49) Guillot, B.; Guissani, Y. J. Chem. Phys. 2001, 114, 6720. (50) Zhu, S.-B.; Singh, S.; Robinson, G. W. J. Chem. Phys. 1991, 95, 2791. (51) Halley, J. W.; Rustad, J. R.; Rahman, A. J. Chem. Phys. 1993, 98, 4110. (52) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (53) Jorgensen, W. L. J. Am. Chem. Soc. 1981, 103, 335. (54) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 83, 926. (55) Chialvo, A. A.; Cummings, P. T. Fluid Phase Equilib. 1998, 150151, 73. (56) Paricaud, P.; Predota, M.; Chialvo, A. A.; Cummings, P. T. J. Chem. Phys. 2005, 122, 244511. (57) Stern, H. A.; Rittner, F.; Berne, B. J.; Friesner, R. A. J. Chem. Phys. 2001, 115, 2237. (58) Mahadevan, T. S.; Garofalini, S. H. J. Phys. Chem. 2007, 111, 8919. (59) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601. (60) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. J. Chem. Phys. 1999, 110, 8254. (61) Garofalini, S. H. J. Chem. Phys. 1983, 78, 2069. (62) Feuston, B. F.; Garofalini, S. H. J. Chem. Phys. 1989, 91, 564. (63) Levine, S. M.; Garofalini, S. H. J. Chem. Phys. 1987, 86, 2997. (64) Galeener, F. L. J. Non-Cryst. Solids 1982, 49, 53. (65) Galeener, F. Vibrational Evidence for Intermediate Range Order in Glasses. In The Structure of Non-Crystalline Materials; Gaskell, P. H., Parker, J. M., Davis, E. A., Eds.; Taylor and Francis Inc.: New York, 1983; p 337. (66) Bunker, B.; Haaland, D.; Michalske, T.; Smith, W. Surf. Sci. 1989, 222, 95. (67) Bunker, B. C.; Haaland, D. M.; Ward, K. J.; Michalske, T. A.; Smith, W. L.; Balfe, C. A. Surf. Sci. 1989, 210, 406. (68) Zhuravlev, L. T. Langmuir 1987, 3, 316. (69) Garofalini, S. H.; Zirl, D. M. J. Vac. Sci. Technol., A 1988, 6, 975. (70) Garofalini, S. H. Atomistic Structure and Dynamic Behavior of Silica Surfaces. In Proceedings of the Second International Conference on Semiconductor Wafer Bonding; Schmidt, M. A., Abe, T., Hunt, C., Baumgart, H., Eds.; The Electrochemical Society: Pennington, NJ, 1993; Vol. 93-29; p 57. (71) Kohler, A. E.; Garofalini, S. H. Langmuir 1994, 10, 4664. (72) Nogami, M.; Abe, Y. Phys. ReV. B 1997, 55, 12108. (73) Nogami, M.; Nagao, R.; Wong, C.; Kasuga, T.; Hayakawa, T. J. Phys. Chem. 1999, 103, 9468. (74) Daiko, Y.; Kasuga, T.; Nogami, M. Microporous Mesoporous Mater. 2004, 69, 149. (75) Nogami, M.; Tarutani, Y.; Daiko, Y.; Izuhara, S.; Nakao, T.; Kasuga, T. J. Electrochem. Soc. 2004, 151, A2095. (76) Schober, T. Ionics 2006, 12, 131. (77) Nogami, M.; Nagao, R.; Wong, C. J. Phys. Chem. 1998, 102, 5772.

Dissociative Chemisorption of Water onto Silica Surfaces and

Recommend Documents