Oxygen Adsorption on β-Quartz Model Surfaces: Some Insights from

Feb 1, 2012 - Oxygen Adsorption on β-Quartz Model Surfaces: Some Insights from Density Functional Theory Calculations and Semiclassical ...
1 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCA

Oxygen Adsorption on β-Quartz Model Surfaces: Some Insights from Density Functional Theory Calculations and Semiclassical TimeDependent Dynamics Costantino Zazza,†,§ Maria Rutigliano,*,‡ Nico Sanna,† Vincenzo Barone,§ and Mario Cacciatore‡ †

CASPUR Consorzio Interuniversitario Applicazioni Supercalcolo per Università e Ricerca, Via dei Tizii 6b, 00185 Roma, Italy CNR-IMIP, Istituto di Metodologie Inorganiche e dei Plasmi, Via Amendola 122/D, 70126 Bari, Italy § Scuola Normale Superiore di Pisa, Piazza dei Cavalieri 7, 56126 Pisa, Italy ‡

ABSTRACT: The O/β-quartz interaction is described by combining our timedependent semiclassical approach to atom−molecule/surface scattering with firstprinciples electronic structure calculations at the DFT (PBE0) level of accuracy. In particular, the O, O2 interaction potentials with an on-top Si atom and its nearest O atom both localized over three different silica clusters have been calculated as a function of the oxygen−silica approaching distance. The calculated DFT potential energy surface has been used in semiclassical trajectory calculations to investigate the sticking and inelastic reflection of oxygen atoms from a model β-quartz surface. The collisional mechanism, including the role played by the phonon dynamics, is brought to light and accurate sticking probabilities are calculated at five impact energies in the range [0.05−0.8] eV and TS = 1000 K. The different catalytic response of β-quartz and β-cristobobalite to the atomic oxygen flux is also discussed and highlighted.

I. INTRODUCTION Recently, the interaction of atomic and molecular oxygen with silica and silica-based surfaces has been the object of various experimental and theoretical studies, due to the importance that this heterogeneous system has in various research fields and, specifically, in aerospace.1 Silica is, in fact, regarded as being among the less catalytic high-temperature resistant materials capable of inhibiting the exothermic surface caused by the interaction of atomic oxygen and nitrogen on the space shuttle’s walls during re-entry into the Earth’s atmosphere. As a consequence, the thermal protection system (TPS) of space shuttles used in planetary exploration missions is still made up mostly of silica-based materials. In this respect, the determination of the real heat flux, needed for realistic tailoring of the TPS system in exercise conditions, still represents a complicated task for the scientific community. From an experimental point of view, this is due to the difficulty of reproducing in the laboratory the extreme chemical−physical conditions typically encountered during the re-entry phase trajectory. On the other hand, kinetics or fluid-dynamics simulations are still approximate and far from being exhaustive; this is basically because there are many heterogeneous molecular processes that may occur at the boundary layer of the shuttle walls and, most importantly, they are still rather poorly understood.2 Among the many possible elementary surface processes at the silica−air interlayer, the oxygen−silica interactions can lead to either nonreactive or reactive processes. Nonreactive processes for instance include inelastic scattering phenomena in which the gas-phase oxygen atoms upon having hit the © 2012 American Chemical Society

surface are then scattered back into the gas phase. On the other hand, reactive adsorption channels Ogas + SiO2 →Oad *SiO2

(1)

can be observed when oxygen atoms approaching the silica surface from the gas phase are trapped in a silica chemi- or physisorption site and the available energy is not sufficient to enable them to escape from the potential energy well. The O adsorption processes are the precursors of another important class of reactive processes, i.e., the formation of energetically activated O2 molecules due to atom recombination on silica.3 In fact, this reaction, whether of Langmuir− Hinshelwood or of Eley−Rideal type, follows a two-step mechanism, the first of which involves the preadsorption of the oxygen atom at the silica surface. The recombination reactions can be highly exothermic and are therefore generally considered to be the main source responsible for a potentially relevant heat flux to the silica substrate that constitutes the space shuttle TPS covering the orbiter surface. A first study dealing with the reactive interaction of oxygen atoms recombining on a β-cristobalite surface from firstprinciples molecular dynamics calculations was undertaken by us.4 In this, and subsequent papers,5−7 the complexity of the O−silica interaction was highlighted, a complexity mainly caused by the need to include various factors controlling the dynamics of the atom−silica interaction in the MD simulation. Received: June 13, 2011 Revised: January 27, 2012 Published: February 1, 2012 1975

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

Further light is shed on other basic aspects concerning the O−silica system. Comparing the results with those obtained previously for the O−cristobalite system, the correlation between surface morphology and catalytic response is pointed out and discussed. Both aspects, i.e., the different electronic response and rearrangement of the two polymorph surfaces to the interacting oxygen atom and the consequently different adsorption dynamics, are discussed. In addition, we focus on a further aspect of the interaction, i.e., the different response of the phonons of the two different polymorphs to the interacting oxygen atom. Indeed, it is the phonons of the substrate that take part and, ultimately, control and determine the adsorption process. The inclusion of phonon dynamics in the dynamics of the molecule−surface interaction is one of the most interesting aspects of heterogeneous processes, generally overlooked in quantum calculations of an already complex nature. Another result of our analysis, relevant for the above-mentioned spaceshuttle thermal shield, is the great exothermicity of the adsorption processes. The oxygen chemisorption energy, which is different for the two polymorphs, is in fact mostly transferred to the phonons as a heat flow through the silica substrate, an effect not usually taken into consideration for simulations of the boundary layer thermofluid dynamics. Finally, it is worth noting that, although well-known in classical catalysis, the correlation between surface morphology and subsequent catalytic activity has been less thoroughly investigated theoretically using first-principles approaches.

One of the most influential factors in heterogeneous catalysis obviously concerns the atomistic structure and the related morphology of the surface exposed to the interacting gas-phase particles. As far as the O−silica interaction is concerned, the reactivity of two different polymorphs of silica, β-cristobalite and βquartz, was recently investigated5 across a large collisional energies range. The results found in ref 5 show that the crystalline structure of silica is determinant for the oxygen atom recombination process on silica. In fact, the recombination coefficient γ for O2 formation, at TS = 1000 K, was 3−4 times lower for the quartz surface than for cristobalite under the same recombination conditions. The phase stability diagram for silica shows that β-quartz, under standard pressure, is thermodynamically stable in the temperature range between 846 and 1145 K, whereas βcristobalite is stable up to the silica liquefaction temperature. The transitions between the different silica phases involve breaking up and re-forming the intra-atomic bonds with significant displacements of the atoms within the crystalline phase. Thus, at high temperatures, β-cristobalite has a facecentered cubic structure, whereas β-quartz belongs to the compact hexagonal system. The different rearrangement of the SiO4 tetrahedra within the crystal network, and the consequent diverse local electronic structure of the SiO2 units, is responsible for the different chemico-physical properties characterizing the two polymorphs. Despite the key role played by atom adsorption processes on substrate reactivity, and unlike the recombination processes the adsorption of atomic oxygen on silica surfaces has not yet been investigated experimentally, either directly or indirectly, due to obvious experimental difficulties. Hence, the need to explore in detail this key silica surface process on a theoretical ground. In the absence of experimental observations, results obtained using MD calculations, if obtained using reasonably accurate theoretical/computational methods, can be regarded as being of a highly predictive nature, as well as providing otherwise inaccessible collisional data on the dynamics and energetics of atom adsorption processes. In a previous work,7 we addressed this problem by performing accurate molecular dynamics simulations capable of describing in great detail the interaction of oxygen atoms with β-cristobalite. In the present work, which should be regarded as a further step in our investigation into O,O2−silica surface processes, we focus on the adsorption of atomic oxygen on a β-quartz surface, the other high-temperature crystallographic polymorph of silica. The methodology followed in this study is the same as that developed in another paper,7 so the O/β-quartz interaction is described by combining our time-dependent semiclassical approach to atom-molecule/surface scattering8 with firstprinciples electronic structure calculations at the DFT level of accuracy. The objectives of the work are many-fold. On the one hand, our global approach to O/β-quartz dynamics allows us to investigate in detail the atomistic steps underlying the chemisorption process, including the phonon excitation− annihilation processes of the silica substrate. The collisional mechanism, whether direct or indirect, is therefore brought to light and accurate adsorption probabilities are calculated at different impact kinetic energies of the O atom hitting the silica surface. The opening of inelastic processes other than adsorption is also investigated.

2. INTERACTION POTENTIALS OF O AND O2 OVER β-QUARTZ AND β-CRISTOBALITE SURFACES Empirical/classical molecular force fields are often further refined by introducing interaction potentials based on correlated electronic structure calculations accomplished using either conventional wave function based approaches9 or density functional theory (DFT) methods10,11 In this context, in a couple of recent applications we have recently and successfully applied DFT modeling for studying the O, N, N2, and O2 interaction patters over β-cristobalite polymorph.6,7 Furthermore, in the case of atomic species (i.e., O and N), the DFT minimum potential energy profiles computed considering sticking directions to the surface model have been explicitly introduced in a subsequent time-dependent semiclassical dynamics approach, with the scope of improving the computational description of the atom−surface local interactions during collisional events.6,7 Interestingly, in all the cases considered in our DFT computations, the basic picture that emerged reveals a complex chemical adsorption process both for the oxygen and for the nitrogen atoms. In fact, under our simulation conditions the chemical adsorption of O and N atoms over β-cristobalite was found to be strongly “spin-state-selected” with changes in the electronic configuration of the whole investigated model systems depending on the distances between the approaching atoms and the surface itself. In light of this observation, it therefore seems reasonable that any serious modeling must necessarily maintain the intrinsic electronic flexibility of these molecular systems during collisions. In turn, a fully classical description of the collisional pathways without any quantummechanics refinement of the N, O, N2, and O2 interactions with silica polymorphs should be considered, at least in principle, as affected by such a limitation during the attachment to the surface of the adsorbant.6,7 1976

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

In this contribution, inspired by the above computational evidence,6,7 we have addressed the O and O2 interaction potentials with an on top silicon atom and its nearest oxygen both localized over predefined β-quartz surface models cleaved from the unit cell as reported in ref 12 (Figure 1). To this end,

experimental coordinates; in particular, O and O2 gaseous species have been considered in interaction with Si3O4H6, Si8O16H14, and Si14O28H14 clusters (Figure 1). Furthermore, in line with previous investigations and with the aim to address the O and O2 interaction with an on-top Si atom over the βquartz, a potential energy surface (PES) scan has been made for each system model in the range RO,Si 1.00−4.00 Å; in this respect, it is worthwhile to mention that, as already applied in our recent work,7 the O−O interatomic distance used in all O2−surface PES scans was fixed to the experimental value reported in literature, that is, RO−O = 1.20 Å16 (Figure 1), along a perpendicular to the surface (i.e., sticking) direction of impact. Therefore, as a conclusive step we have also explicitly considered in our simulating scenario the oxygen(gas phase)-oxygen(β-quartz) interaction at the DFT level of theory. As a matter of fact, DFT-PES scans along the sticking direction of impact connecting the oxygen present in the gas phase with the oxygen atom onto the β-quartz surface have been computed (see the molecular structures depicted in Figure 2), and then

Figure 1. Cluster surface models of the β-quartz polymorph (i.e., Si3O4H6, Si8O16H14, and Si14O28H14) considered in our simulation conditions in interaction with O and O2 species both coming from the gas-phase condition and along a sticking direction of impact. Silicon atoms are reported in yellow, oxygen atoms are in red, and hydrogen atoms are in white.

Figure 2. Si3O4H6 (upper panel) and Si8O16H14 (lower panel) βquartz surface models simulated in interaction with the incoming gaseous O atom along a sticking R(Ogas‑phase−Osurface) reaction coordinate (see red dotted arrows). Colors as reported in Figure 1.

reported in section 3. It is worth noting that, for high precision PESs, extremely tight convergence criteria were used and an ultrafine grid has been chosen for DFT numerical integrations. PBE0/6-311+G* PES computations were run on computational facilities at Caspur (Rome) using the Gaussian03 package.17 As first step we started considering the gaseous oxygen atom in interaction with the on-top silicon of the smallest cluster of the β-quartz polymorph reported in Figure 1 (upper panel), that is, Si3O4H6, along the sticking direction and the results are then analyzed and compared with those previously obtained by taking into account the β-cristobalite surface model. As clearly reported in Figure 3 (upper panel), while the incoming atom

and for a direct comparison with previous data obtained for βcristobalite, we used the hybrid PBE0 functional13 based on Perdew−Burke−Erzenrhof exchange−correlation functionals14 in which a portion of Hartree−Fock exchange is added selfconsistently to the DFT (PBE) contribution. A triple split valence Pople style atomic basis set with additional diffuse spfunctions and polarization d-functions on heavy atoms, e.g., 6311+G*,15 has been selected for all our computations reported below. Then, three different cluster model systems of the βquartz silica surface have been selected for this purpose while keeping their geometries blocked according to the available 1977

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

theory considering free 1[Si3O4H6] and 3[Si3O4H6] systems. In this respect, according to these DFT potential energy curves the O atom from the gas phase, though it interacts with the substrate (i.e., Si3O4H6), does not alter its electronic configuration, remaining always confined in its own electronic ground state (i.e., 3P).16 The same trend has also been observed when we analyzed the potential energy curves calculated by increasing the size of the surface models to Si8O16H14 and Si14O28H14. The only remarkable difference observed in our computations is a slight increase of the binding energy from −5.51 eV (Si3O4H6) to −5.56 eV (Si8O16H14) and −5.58 eV (Si14O28H14), as can be observed from the data collected in Table 1. On the other hand, from this table we also note that Table 1. R(O−Si) Minimum Distances and Interaction Energies (Termed Be and Reported in eV) for the O− SixOyHz Model Systems of the Oxygen Sticking over βQuartz and β-Cristobalite Silica Polymorphs Computed with the PBE0/DFT Functional in Conjunction with the 6311+G* Atomic Basis Set β-quartz

Be/eV

Rmin (O−Si)/Å

O−Si3O4H6 O−Si8O16H14 O−Si14O28H14 β-cristobalite7

−5.51 −5.56 −5.58 Be/eV

1.53 1.52 1.52 (O−Si)/Å

O−Si3O4H6 O−Si7O14H14 O−Si17O34H18 ref 18a

−5.47 −5.60 −5.57 −5.60

Rmin

1.52 1.51 1.52 1.56

the minimum energy interatomic distance between the Ogas and the Si atom on top (termed as Rmin(O−Si)) remains almost constant as a function of the cluster size (Si3O4H6 → Si8O16H14 → Si14O28H14). When the O2 molecule is concerned, the potential energy minimum profile for the O2−Si3O4H6 system is found to involve, at a distance between the nearest atom to the surface oxygen of the O2 and the on top Si in the range 1.75−4.00 Å, two triplet states in which the incoming molecule can be found in its ground state, that is, X 3Σg−, and in interaction with the 1 [Si3O4H6] cluster (profile termed as tripleta in Figure 3, lower panel), or in a 1Σg+ state in combination with a 3[Si3O4H6] moiety (tripletc). At shorter distances, the presence of a spin crossing lying at almost 1.74 Å enables a second spin variation toward a singlet spin multiplicity curve (singlet profile in Figure 3, lower panel) with an O2 molecule featuring a 3Σ u+ occupation plus a 3[Si3O4H6] surface fragment. Intriguingly, the O2−β-cristobalite system showed, as explicitly reported in Figure 4 and in ref 7, virtually the same minimum energy profile computed for the O2−β-quartz one. Next, current results are more deeply compared with those recently obtained for a series of β-cristobalite cluster models having comparable dimensions (i.e., Si3O4H6, O−Si7O14H14 and Si17O34H18).7 This analysis clearly evidences that the two silica polymorphs (βquartz and β-cristobalite) show essentially the same PESs with respect to the oxygen sticking trajectory (Table 1 and Figure 4). This is due to the fact that the chemical interaction appears to be mainly confined within the silicon atom on top and, to a lesser extent, to its nearest oxygen atoms. In other words, our DFT results seem to suggest that along such a direction of impact the gaseous O atom and O2 molecule are not capable of discriminating between the two different silica substrates as a

Figure 3. Upper panel: oxygen atom interaction potential energy scan (PES) over the Si atom of Si3O4H6 cluster model of β-quartz surface; triplet (dotted line) and singlet (solid line) eigenstates of total electronic spin curves at PBE0/6-311+G* level of computation. Lower panel: minimum energy profiles, at PBE0/6-311+G* level of computation, for O (solid line) and O2 (dashed line) gas species while approaching the Si atom of the Si3O4H6 cluster model of the βquartz surface.

approaches the surface model, the entire investigated system undergoes a change in the electronic configuration as a response to the interaction between the O and Si3O4H6 moieties. More precisely, at R(O−Si) distances as long as 2.325 Å, the O−Si3O4H6 system shows a smooth crossing between triplet and singlet spin multiplicity profiles, which is ascribed to a 1[Si3O4H6] → 3[Si3O4H6] electronic rearrangement. Such a change in the electronic configuration of the substrate is known to be a key route for catalyzing the chemisorption processes in the presence of oxygen4,7 and nitrogen atoms.6 As a matter of fact, upon the 1[Si3O4H6] → 3 [Si3O4H6] change, the β-quartz polymorph shows two unpaired electrons that can be involved in the formation of two chemical bonds with the incoming oxygen: a polarized σ bond along the RO−Si reaction coordinate and a π Si−O double bond.18 This 1[Si3O4H6] → 3[Si3O4H6] local electronic response occurring in the β-quartz surface model has been further confirmed by the asymptotic energy difference between the two curves reported in Figure 3 (upper panel), which is found to be equal to 2.50 eV: this value perfectly agrees with the energy difference computed at the PBE0/6-311+G* level of 1978

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

while interacting with the two investigated silica polymorphs. In this context, we have seen that, at RO−O minimum distance of 2.18 Å for β-cristobalite and 2.17 Å for β-quartz the RO(gas phase)−Si(on top) interatomic distance is found to fall at 1.74 Å and at 1.95 Å, respectively. From this observation it is pertinent to underline one more time that, within our simulation conditions, the RO(gas‑phase)−O(surface) interaction profiles are seen to be strongly fine-tuned by the presence of an adjacent and unsatured silicon atom (that is the so-called on top Si atom). In the light of these findings, it has been therefore needed for a realistic description of the collisional pathways characterizing the dynamics of the oxygen interaction with βcristobalite and β-quartz, to take properly into account the two different RO(gas‑phase)−O(surface) PESs as reported in Figure 5. Before concluding this section, let us clarify that our choice to retain a constant O−O distance in the oxygen molecule as it is approaching the surface is justified simply because we are disregarding recombination processes catalyzed by the β-quartz substrate in the subsequent semiclassical dynamics (see section 3). Also, it is important to stress that this apparently unrealistic situation would seem to point out the potential involvement of several electronic eigenstates during recombination and/or dissociative chemisorption of the O2 species over β-quartz surface. A more complete and realistic DFT picture in the presence of a flexible oxygen molecule will represent the natural extension of the present work. Preliminary PBE0/6-311+G* data evidence that a fixed O−O distance represents a good approximation for values longer than 2.20 Å from the surface. At shorter distances the lowest energy profile has always a triplet character with the O−O distance changing from 1.20 to 1.40 Å before undergoing a dissociative chemisorption assisted by a singlet−triplet spin conversion of the surface model.

Figure 4. Minimum potential energy profiles, at PBE0/6-311+G* level of computation, for O and O2 incoming species over the Si atom (Figure 1) of the β-quartz (Si8O16H14, solid line) and β-cristobalite (Si7O14H14, dashed line) surfaces, respectively.

consequence of strongly localized interaction pattern. On the other hand, potential energy curves computed with respect to the RO(gas‑phase)−O(surface) collisional coordinate (see the starting conformation depicted in Figure 2) show that the interaction process along the sticking trajectory over β-quartz is fairly different from that observed for the β-cristobalite. In fact, as can be clearly seen in Figure 5, the chemisorption process of the O

3. OXYGEN ATOM ADSORPTION DYNAMICS OVER β-QUARTZ SURFACE The adsorption dynamics was carried out in the framework of the semiclassical collisional method8 according to which the O−silica interaction dynamics is assisted and controlled by the phonon excitation and de-excitation processes in the solid substrate. Therefore, as a first step, a 3D β-quartz model crystal was built up, consisting of 198 atoms (66 of Si) disposed on 8 layers according to the atom placement in the crystallographic unit cell.12 The phonon dynamics of the isolated cluster was then developed by numerical diagonalization of the 3Ddynamical matrix of the force constants obtained assuming the BKS potential19 as interaction potential of the Si and O lattice atoms. The phonons eigenvalues and eigenvectors that enter the scattering equations of motion are then obtained. The calculated phonon distribution function was reported in ref 5 and, as stated there, is in qualitative agreement with the phonon spectrum reported in ref 20. According to the semiclassical method, the O−silica interaction potential that governs the dynamics of the processes going on the surface is given as sum of pairwise atom−atom interactions between the gas-phase oxygen atom and the Si and O of the quartz surface. In the determination of complete potential energy surface (PES) we have considered the interaction potential obtained in DFT calculations for the O gas-phase atom impinging on the Si atom in the lattice as well

Figure 5. PBE0/6-311+G* calculated O(gas-phase)−O potential energy scan (PES) over the β-quartz (Si8O16H14, solid line) and βcristobalite (Si7O14H14, dashed line) surface models, respectively.

atom shows a higher binding energy when it is supposed to occur on a β-cristobalite cluster model; a binding energy of −2.08 eV was computed for the O(gas)/(O) β-cristobalite system, whereas the currently interacting O(gas)/(O) β-quartz one results in a binding energy of −1.05 eV. Notwithstanding such a different energetic aspect, it is, however, interesting to observe that the potential energy profiles are both characterized by the lack of any energetic barrier during the collision. The discrepancy in the energetic value (Be = −2.08 eV over βcristobalite vs Be = −1.05 eV over β-quartz), as well as the different interaction trend computed at PBE0/6-311+G* level of computation can be ascribed to a similar, but not equivalent, local atomistic environment sampled by the incoming atom 1979

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

perpendicular geometry is reported as a function of the O(gasphase)−Si distance R. The interaction potential for an O gasphase atom interacting on top of a lattice O atom is also reported as a function of the O(gas phase)−O distance R. The corresponding reference ab initio potentials reported in Figures 4 and 5 are also shown. The scattering equations are solved with respect to a Cartesian frame of reference centered on a Si surface atom having the Z-axis normal to the silica surface. The (X−Y) plane lies on the topmost surface layer. In the simulation the O atom is initially placed at a distance Z = 7 Å, far enough into the asymptotic free region. In this study we focus on the dynamics of the O atom approaching the Si site of the β-quartz surface with a specific initial geometry. In particular, we assume that the gas-phase O collides with the surface in the perpendicular geometry (polar angles (θ, φ) = (0, 0)) with a given kinetic energy, Ekin. The initial position coordinates (X, Y) are chosen randomly within an aiming area of 2.0 Å2 centered on the active silicon site. The dynamics was carried out at impact energies in the range [0.05−0.8]eV, including the same energies considered in our previous study on β-cristobalite. This collisional regime allows us to compare the results obtained in the two studies and, therefore, to point out the influence of the different surface polymorph on the adsorption dynamics. Also the temperature of the β-quartz was kept constant to TS = 1000 K, the same temperature assumed in ref 7. For each impact energy a batch of 5000 trajectories was integrated with a time step of 10−14 s. This ensures a numerical convergence of the calculated adsorption probabilities of about 5%. As expected, the trajectory analysis of the final states points toward the existence of two main scattering processes: (i) the O atom can be adsorbed at the quartz surface or it can be (ii) inelastically scattered back into the gas phase after hitting the surface. Process (i) occurs via a “direct” mechanism, according to which the oxygen atom is trapped at the surface after a single collision, or via an “indirect” multistep mechanism according to which the oxygen atom collides several times with the quartz surface before adsorption. Also the inelastic scattering (ii) can follow two different mechanisms: the O atom can be scattered directly into the gas phase or it can be temporarily adsorbed before desorption. We assume that O is adsorbed when the energy available to the O atom to escape from the chemisorption potential well, ΔEesc, is less than the effective potential, ΔVadd, experienced by the oxygen at a given interaction time. ΔEesc is defined as Etot − ΔEph − V0, where V0 is the zero point energy of the chemisorbed atom. Further, a second criterion can be followed according to which O is adsorbed whenever the Z component of the atom reaches the chemisorption distance, Zad = 1.52 Å, and the component of the momentum along the Z direction is negative for a sufficiently long time, of the order of few femto seconds. On the other hand, O desorbs when it is scattered from the silica surface far into the free space, i.e., when in the final condition of the trajectory the Z component is larger than 7.00 Å. In Table 3 the probabilities for inelastic reflection and adsorption processes are reported at the collisional energies explored in this study. From Table 3 we notice that, in the considered collisional energy range, the probability for adsorption decreases with the collisional energy, the contrary of what occurs for the reflection.

as for the O gas-phase atom impinging on the O atom in the lattice as reported in Figures 1 (upper panel) and 2. NSi

NO

i=1

i=1

VO−SiO2 = [ ∑ VO−Si]ft + [ ∑ VO−O](1 − ft ) ft = 1. − 1.17 tanh

(XO2 + YO2)

(2) (3)

where f t is a switching function and (XO, YO) are the coordinates of the gas-phase atom in the (X−Y) plane of the assumed frame of reference. The potential in eqs 2 and 3 features the same functional form to that used in ref 7 for the O−β-cristobalite system. Then we fitted these curves with Morse-like functions (eqs 4 and 5) and we combined them to obtain the complete PES that enables us to reproduce the corresponding DFT curve for the two different interaction geometry (Ogas−Sibulk and Ogas−Obulk) as a function of Ogas distance from the surface. The analytical expression for the Ogas−Sibulk and Ogas−Obulk interactions are, respectively, VO−Si = DO−Si e−pO−Si (R O−Si − r1)(e−pO−Si (R O−Si − r1) − 2)

(4)

VO−O = DO−Oe−pO−O (R O−O− r2)(e−pO−O (R O−O− r2) − 2)

(5)

where RO−Si and RO−O are the distances between the gas-phase oxygen atom and the Si and O lattice atoms, respectively. The Morse parameters, DO−Si, pO−Si, DO−O, pO−O, r1, and r2, are reported in Table 2. Table 2. Best-Fit Potential Parameters of Eqs 4 and 5 DO−Si (eV)

pO−Si (Å−1)

r1 (Å)

DO−O (eV)

pO−O (Å−1)

r2 (Å)

5.46

2.2

1.53

0.48

1.74

3.22

In Figure 6 the analytical interaction potential for atomic oxygen interacting with the on top Si lattice atom in the

Figure 6. Interaction potential for a gas-phase O atom interacting perpendicularly with the on top Si (black line) and O (gray line) atom of β-quartz (Figures 1 and 2) as a function of the O(gas-phase)−Si and O(gas-phase)−O distances: DFT-GGA calculations, full line; analytical fit, dashed line. 1980

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

Table 3. Semiclassical Probabilities for Inelastic Adsorption and Reflection Processes Due to O−Silica (β-Quartz) Interaction at Three Impact Energies (in Italics, βCristobalite) (TS = 1000 K) Ekin (eV)

adsorption

reflection

0.05

0.983 0.97 0.893 0.863 0.652 0.941 0.648 0.822 0.635 0.787

0.017 0.029 0.107 0.137 0.348 0.058 0.351 0.178 0.365 0.212

0.1 0.2 0.5 0.8

Figure 7. Time evolution of a typical trajectory for a multistep O−Si β-quartz adsorption process. The Z normal coordinate as a function of time for O atom impinging on the surface (left scale) and the energy transferred to the surface (right scale) are reported as a function of the collision time. Ekin = 0.2 eV and TS = 1000 K (1τ = 10−14 s).

This behavior is different from that observed in the case in which an oxygen atom impinges on a β-cristobalite surface (reported in italics). In that case, the most effective channel is the adsorption process for which the probability is nearly unitary for the lowest kinetic energy followed by a decreasing as the collisional energy increases. In the O−Si(β-quartz) trajectory analysis we distinguished the direct adsorption process, where the oxygen undergoes a single collision before adsorbing, an indirect adsorption process, where the oxygen atom rebounds several times on the surface before adsorbing. The obtained results are reported in Table 4. Table 4. Probabilities for Single and Multicollision Adsorption in O-Silica Interaction (in Italics, β-cristobalite) (TS = 1000 K)a Ekin (eV)

adsorption (single step)

adsorption (multi step)

0.05

0.807 0.754 0.582 0.397 0.355 0.725 0.602 0.687 0.625 0.721

0.176 0.217 0.311 0.466 0.297 0.216 0.046 0.135 0.010 0.066

0.1 0.2 0.5 0.8 a

Figure 8. Oxygen atom adsorbed (X−Y) coordinates (black dot) on the β-quartz surface (blue ball, Si atoms; red ball, O atoms) model assumed in the calculations. Ekin = 0.2 eV.

the geometry of the quartz model surface used in the simulation. From this figure we can infer that the adsorption occurs at the corners of a square centered on the Si atom. A very interesting aspect of the adsorption process concerns the energetics of the process, that is, the mechanism of energy exchange between the O atom and the phonons of the silica surface. Because the binding energy of O interacting with the Si active site is 5.58 eV, the energy delivered in the interaction process is substantial. The energy delivered to the quartz surface due to the phonon creation/annihilation processes is reported in Figure 9a, whereas the energy gained/lost by the oxygen atom as kinetic energy (Ekin) is reported in Figure 9d. The effective potential (Vadd), and the interaction potential (Vpot) are also reported as a function of the interaction time in Figure 9b,c. The energetics reported in Figure 9a−d corresponds to the multistep adsorption collision shown in Figure 7. At the first bounce on the surface, t ≈ 27τ, the O atom is initially attracted into the potential well due to attractive forces. The O−Si binding energy delivered in the adsorption process is in part transferred to the substrate (Figure 9a) and in part is transferred to the O atom as translational energy. As the collision evolves, the O comes closer to the surface, thus feeling the repulsive branch of the interaction potential. As a consequence, the energy transferred decreases in that the phonons are de-excited while the oxygen atom gains enough kinetic energy to temporarily escape from the potential well before being reattracted for the second time. The energy

The gas-phase O atom impacts perpendicularly the silica surface.

Comparing the results reported in Table 4 with those reported in Table 3 we can infer that, even though the behaviors of total adsorption probability as a function of the collisional energy are very similar for the two polymorphs, the dynamics of the adsorption process turns out to be much more complex when the distinction between the single and multistep processes is considered. The observed differences can be ascribed not only to the different phonon structures but also to a different reactivity of the substrates as reported in ref 5. In Figure 7 a typical trajectory for a multicollision adsorption is reported. Here the Z component of the gas-phase oxygen atom (left scale) together with the energy transferred to the surface (right scale) as a function of the collisional time (in units of τ, 1τ = 10−14 s) are reported. In Figure 8 the (X, Y) coordinates in the final state of the adsorbed O atom in single collisions are reported together with 1981

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

Figure 9. Time evolution for the energy exchanged with the surface (ΔEph), the effective potential (Vadd), the interaction potential (Vpot), and the oxygen kinetic energy (Ekin), for the trajectory in Figure 7. Ekin = 0.2 eV and TS = 1000 K.

exchange mechanism is evidently quite complex and the sequence just described is followed also in the successive bounces. After the third collision, the oxygen atom is finally adsorbed because the residual kinetic energy is not enough to escape from the potential well. It is interesting to note that, at each impact, nearly 1.0 eV is transferred to the silica surface so that at the end of the collision the total energy transferred to the surface is about 3.0 eV. The dynamics followed by the reflection collisions is the same to that observed in the case of O on β-cristobalite surface: a single impact with the surface after which the atom is scattered off into the gas phase. The most interesting result obtained in this comparative study stands out when we turn to consider the energy exchanged with the surface phonons. The average energy transferred is reported in Table 5.

Moreover, for the reflection processes the energy delivered to the quartz surface in reflection collisions is, at Ekin = 0.2 eV, lower than that transferred to the cristobalite surface, the reverse being true at the two higher collisional energies. This is a further manifestation of the large impact that different polymorphs have on the dynamics of atom−surface interactions. Finally, it is worth noting that the average energy transferred to the surface as heat flux is very high, so that a surface thermal damage can be expected. Indeed, it comes out that, very likely, the adsorption process of atomic oxygen could be very effective in the energy dissipation toward the substrate, more effective than the recombination reactions commonly assumed as the main ones responsible for the damage of the thermal protection systems of spatial vehicles re-entering the Earth’s atmosphere.

Table 5. Average Energy Transferred to the Silica Phonons, ΔEph, in Single-Step Adsorption and Reflection of Atomic Oxygen Interacting on Top of a Si Lattice Site (in Italics, βCristobalite) (TS = 1000 K)

4. CONCLUSIONS The adsorption dynamics of atomic oxygen on a model surface of β-quartz silica polymorph has been investigated by means of a time-dependent molecular dynamics semiclassical approach; to this end, the local force field between the incoming atom and the silicon and oxygen atoms of model clusters has been further refined using potential energy surfaces arising from density functional theory (DFT) calculations. As already emerged considering the β-cristobalite system interacting with O and O2 species, we observe that, also in the current case the O,O2-silica interaction is localized within the nearest neighbor atoms selected for the collisions. However, some differences between β-quartz and β-cristobalite silica polymorphs have emerged in our computations. Although, for both cases, while approaching the surface, the incoming O and O2 species may change their electronic configurations, the O/O2(gas phase)−oxygen(βquartz and β-cristobalite) energy profiles along the sticking direction appear to be depended on the individual local geometry characterizing the investigated polymorphs. Such

Ekin (eV)

single step adsorption

reflection

0.2

1.670 0.237 2.030 0.375 2.640 0.548

0.293 0.373 0.579 0.358 0.620 0.478

0.5 0.8

Comparing Table 5 with the results reported in Table 5 of ref 7, and here reported in italics, one can note that ΔEph transferred to the quartz surface in the adsorption processes is larger than that exchanged with the β-cristobalite surface by up to a factor of about 7. The phonon energy in the case of βquartz is systematically higher than that exchanged with the cristobalite surface. 1982

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983

The Journal of Physical Chemistry A

Article

evidence arising from a different surface configuration for βquartz and β-cristobalite coupled with a different phonon spectrum allows us to identify a rather different reactivity in the presence of incoming oxygen atoms at the considered kinetic energies. Interestingly, in this investigation we clearly report that the energy transferred to the β-quartz surface in the most statically relevant reactive routes (i.e., single step adsorption) is much higher than that calculated for the β-cristobalite. This finding let us suppose that a surface thermal damage should be presumably more effective in the β-quartz as a result of an efficient chemical conversion scheme (up to almost 80% at Ekin = 0.8 eV) capable of funneling the kinetic energy of the gaseous atomic oxygen into the polymorph structure itself.



(15) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650−654. (16) Herzberg, G. Spectra of Diatomic Molecules; Molecular Spectra and Molecular Structure Vol. 1; Van Nostrand Reinhold: New York, 1950. (17) Frisch, M. J.;et al. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (18) (a) Arasa, C.; Busnego, H. F.; Salin, A.; Sayòs, R. Surf. Sci. 2008, 602, 975−985. (b) Arasa, C.; Gamallo, P.; Sayòs, R. J. Phys. Chem. B 2005, 109, 14954−14964. (19) Van Beest, B. W. H.; Kramer, G. J.; Van Santen, R. A. Phys. Rev. Lett. 1990, 64, 1955−1958. (20) Hammonds, K. D.; Bosenick, A.; Dove, M. T.; Heine, V. Am. Mineral. 1998, 83, 476−479.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. A movie of trajectory in Figure 7 can be delivered upon request.



ACKNOWLEDGMENTS This research was supported by the Italian Ministry of Education, University and Research in the framework of the PRIN N. 2008KJX4SN Project. This work was inspired by the research made in the framework of the CIRA-IMIP project “CAST- Innovative Aerothermodynamic Configurations for Systems of Spatial Transport”



REFERENCES

(1) Jumper, E. J.; Seward, W. A. J. Therm. Heat Transf. 1994, 8, 460. Kurotaki, T. AIAA J. 2000, 2000−2366. (2) Cacciatore, M.; Rutigliano, M. Plasma Sources Sci. Technol. 2009, 18, 023002. Nasuti, F.; Barbato, M.; Bruno, C. J. Therm. Heat Transf. 1996, 10, 131−136. (3) Cacciatore, M.; Rutigliano, M. Molecular dynamics simulation of surface processes: oxygen recombination on silica surfaces at high temperatures. Experiments, Modelling and Simulation of Gas-surface Interactions for Reactive Flows in Hypersonic Flights; Educational Notes RTO-EN-AVT-142 Editor: Research and Technology Organisation (NATO) BP 25, F-92201 Neuilly-sur-Seine Cedex, France, ISBN 97892-837-0057-9, Paper 5 (2007). (4) Cacciatore, M.; Rutigliano, M.; Billing, G. D. J. Therm. Heat Transf. 1999, 13, 195−203. (5) Bedra, L.; Rutigliano, M.; Cacciatore, M.; Balat-Pichelin, M. Langmuir 2006, 22, 7208−7216. (6) Rutigliano, M.; Pieretti, A.; Cacciatore, M.; Sanna, N.; Barone, V. Surf. Sci. 2006, 600, 4239−4246. (7) Rutigliano, M.; Zazza, C.; Sanna, N.; Pieretti, A.; Mancini, G.; Barone, V.; Cacciatore, M. J. Phys. Chem. A 2009, 113, 15366−15375. (8) See for example: Cacciatore, M.; Billing, G. D. Surf. Sci. 1990, 232, 35−50. Billing, G. D. Dynamics of Molecule Surface Interactions; Wiley & Sons: New York, 2000. (9) See for example: Dovesi, R.; Saunders, V. R.; Roetti, C.; Causà, M., Harrison, N. M.; Apra, E. CRYSTAL 95 Usre’s Manual; University of Turin: Turin, 1996. (10) Nakamura, K. G. Chem. Phys. Lett. 1998, 285, 21−24. (11) Brivio, G. P.; Grimley, T. B. Surf. Sci. Rep. 1993, 17, 1−84. (12) Wyckoff, R. W. G. Crystal Structures; Wiley & Sons: New York, 1963, 13−25. (13) (a) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (b) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396−1396. (14) Adamo., C.; Barone, V. J. Chem. Phys. 1999, 110, 6158−6170. 1983

dx.doi.org/10.1021/jp205517j | J. Phys. Chem. A 2012, 116, 1975−1983