Scattering of N2 Molecules from Silica Surfaces: Effect of Polymorph and Surface Temperature

The inelastic scattering of N2 molecules from silica surfaces, taken at 100 K, has been investigated by adopting a semiclassical collision model in conjunction with the appropriate treatment of the long-range interaction forces. Such forces promote the formation of the precursor state that controls all basic elementary processes occurring at the gas–surface interphase. The probabilities for the different elementary surface processes triggered by quartz are determined and compared with those recently obtained for another silica polymorph (cristobalite). In addition, the final roto-vibrational distributions of N2 molecules undergoing inelastic scattering have been characterized. N2 molecules, impinging on both considered surfaces in low-medium vibrational states, preserve the initial vibrational state, while those inelastically scattered are rotationally excited and translationally colder. The surface temperature effect, investigated by raising the temperature itself from 100 K up to 1000 K, emerges more sharply for the cristobalite polymorph, mainly for the molecules impinging in the ground roto-vibrational state and with low collision energies.


Introduction
The inelastic scattering of molecules from surfaces triggers several energy exchanges that can have repercussions on the energetics and collision dynamics in gaseous and plasma environments. Moreover, this process exhibits some selectivity and propensity strictly related to the alignment of impinging molecules and to the surface roughness seen by the gaseous molecules.
Silica is a highly versatile and very common material in natural, laboratory and industrial environments. Traces of SiO 2 are found in grains of spatial origin [15], but silicabased materials are also of interest for the thermal protection panels of space vehicles [16], micro-nano electronics [17], and nano-medicine [18]. Accordingly, the investigation of the interaction of small molecules, such as H 2 , O 2 , N 2 , CO with this kind of surface, appears mandatory for the advancement of knowledge in the different fields of scientific and technological interest.
In ref. [19], a strong correlation has been found between the incidence speed, angle, surface topology and the probability distribution of accommodation coefficients for the The complete V N2-silica potential is obtained by adding a Coulomb term to the interaction component, represented by the ILJ function, accounting for the electrostatic contribution, due to the quadrupole moment of the nitrogen molecule interacting with the anisotropic charge distribution of the surface, according to what has been performed in ref. [20].
In Figure 1a,b, the pure long-range ILJ potential contribution and the complete (including the electrostatic component) V N2-silica potential are reported as a function of the Z CM , the distance along the normal to the surface plane of the molecule centre of mass (CM). The different curves refer to N 2 , with the molecular axis aligned perpendicular and parallel to the surface plane, impinging on a Si atom (top layer). In Figure 1c,d, the same curves are reported for N 2 impinging on an O atom (first sublayer). For an easy and immediate comparison in each panel, the curves for the interaction on the quartz (blue lines) and the cristobalite (red lines) are reported [29]. Looking at Figure 1, we can immediately appreciate that the N 2 molecule interacts more strongly with quartz than with cristobalite. Furthermore, for a given polymorph, a stronger interaction manifests when the molecule interacts with the molecular axis parallel to the surface plane and the interaction directly involves an oxygen atom in the first sublayer (see Figure 1c,d).
Molecules 2022, 27, x FOR PEER REVIEW 3 details are given on the procedure adopted and on the comparison of determined va with data available in the literature. The complete VN2-silica potential is obtained by adding a Coulomb term to the inte tion component, represented by the ILJ function, accounting for the electrostatic contr tion, due to the quadrupole moment of the nitrogen molecule interacting with the an tropic charge distribution of the surface, according to what has been performed in ref.
In Figure 1a,b, the pure long-range ILJ potential contribution and the complete cluding the electrostatic component) VN2-silica potential are reported as a function of the the distance along the normal to the surface plane of the molecule centre of mass (C The different curves refer to N2, with the molecular axis aligned perpendicular and pa lel to the surface plane, impinging on a Si atom (top layer). In Figure 1c,d, the same cu are reported for N2 impinging on an O atom (first sublayer). For an easy and immed comparison in each panel, the curves for the interaction on the quartz (blue lines) and cristobalite (red lines) are reported [29]. Looking at Figure 1, we can immediately ap ciate that the N2 molecule interacts more strongly with quartz than with cristobalite. thermore, for a given polymorph, a stronger interaction manifests when the molecul teracts with the molecular axis parallel to the surface plane and the interaction dire involves an oxygen atom in the first sublayer (see Figure 1c,d). Therefore, for the quartz, N2 in this molecular orientation interacts more stron with its CM impinging on an O atom in the first sublayer than in the case of incidenc a Si atom.
It is also worth noting that the contribution of the Coulombian potential is very s in the case of interaction with the silicon atom, whereas for oxygen it is not even ap ciable. Therefore, for the quartz, N 2 in this molecular orientation interacts more strongly with its CM impinging on an O atom in the first sublayer than in the case of incidence on a Si atom.
It is also worth noting that the contribution of the Coulombian potential is very small in the case of interaction with the silicon atom, whereas for oxygen it is not even appreciable.

Potential Energy Surface Determination
The complete PES has the same functional form as that used in ref. [29], obtained by combining the term V N-silica , describing the interaction of a single N atom with the surface, with that for the interaction of N 2 molecule with the same surface, V N2-silica , through a weight function depending on the interatomic distance r N-N . V N-silica has been obtained by adding the interaction potential for short-range distances, resulting from the fitting procedure of the ab initio DFT calculations in ref. [30], to that for the long-range region modelled as an ILJ interaction potential (see previous Section and Section 3). Since the interaction of the N atom with silica clusters, cleaved by the crystallographic unit cell, has been proven to be a local component involving only the N impinging atom and Si (in the cluster) partners, we assumed the same short-range potential for both polymorphs. Furthermore, here we are interested in studying the interaction of molecules with the two polymorphs and, in the investigated collision energy range, we are expecting a minor influence of atomic potential on the collision dynamics. For all these reasons, the strong short-range potential for one N atom interacting with a Si atom on the surface is that used in ref. [29]. In Figure 2, the weak long-range interaction of N with an O atom in the first sublayer and with a Si atom in the second sublayer is reported as a function of the distance along the normal to the surface of the nitrogen atom, Z N . In the assumed reference frame, Z = 0 is on the top layer.

Potential Energy Surface Determination
The complete PES has the same functional form as that used in ref. [29], obtained by combining the term VN-silica, describing the interaction of a single N atom with the surface, with that for the interaction of N2 molecule with the same surface, VN2-silica, through a weight function depending on the interatomic distance rN-N. VN-silica has been obtained by adding the interaction potential for short-range distances, resulting from the fitting procedure of the ab initio DFT calculations in ref. [30], to that for the long-range region modelled as an ILJ interaction potential (see previous Section and Section 3). Since the interaction of the N atom with silica clusters, cleaved by the crystallographic unit cell, has been proven to be a local component involving only the N impinging atom and Si (in the cluster) partners, we assumed the same short-range potential for both polymorphs. Furthermore, here we are interested in studying the interaction of molecules with the two polymorphs and, in the investigated collision energy range, we are expecting a minor influence of atomic potential on the collision dynamics. For all these reasons, the strong short-range potential for one N atom interacting with a Si atom on the surface is that used in ref. [29]. In Figure 2, the weak long-range interaction of N with an O atom in the first sublayer and with a Si atom in the second sublayer is reported as a function of the distance along the normal to the surface of the nitrogen atom, ZN. In the assumed reference frame, Z = 0 is on the top layer. Looking at Figure 2, it appears that the N atom long-range interactions can be active and not negligible up to the second layer placed at 1.78 Å below the surface top layer. Figure 3 shows two projections of a cut (parallel molecular approach to the surface) of the complete PES on the XY surface plane, corresponding to the molecular center of mass placed at two selected distances from the surface: ZCM = 2 Å and ZCM = 3.5 Å. We can appreciate that with increasing ZCM, the specificity of the sites corresponding to the Si and O atoms building the surface is lost and each SiO2 unit manifests almost as a single center of attraction. Looking at Figure 2, it appears that the N atom long-range interactions can be active and not negligible up to the second layer placed at 1.78 Å below the surface top layer. Figure 3 shows two projections of a cut (parallel molecular approach to the surface) of the complete PES on the XY surface plane, corresponding to the molecular center of mass placed at two selected distances from the surface: Z CM = 2 Å and Z CM = 3.5 Å. We can appreciate that with increasing Z CM , the specificity of the sites corresponding to the Si and O atoms building the surface is lost and each SiO 2 unit manifests almost as a single center of attraction.

The Dynamics of N2 Interacting with Quartz
The interaction of N2(vi,ji) molecules in the ground vibrational (vi = 0) and rotational (ji = 0 for ortho and ji = 1 for para symmetries) states with the quartz surface leads to inelastic scattering, whose probability of occurring rises sharply as Ecoll increases and then it assumes almost constant values at the highest collision energies, independently of vi and ji value. This trend can be seen in Figure 4, where the scattering probability for ortho and para nitrogen molecules, approaching along the normal to the surface plane, at the two considered surface temperatures, is reported. For the three lowest collision energy values (Ecoll ≤ 0.1 eV), the integration of most of the trajectories is stopped as the molecule reaches the edge of the crystal lattice by diffusing on the surface through single or multi-rebounds. These trajectories are those for which it is not possible to assign with certainty the final reaction channel, as they are neither completely adsorbed nor completely scattered off in the gaseous-phase. Instead, for Ecoll ≥

The Dynamics of N 2 Interacting with Quartz
The interaction of N 2 (v i ,j i ) molecules in the ground vibrational (v i = 0) and rotational (j i = 0 for ortho and j i = 1 for para symmetries) states with the quartz surface leads to inelastic scattering, whose probability of occurring rises sharply as E coll increases and then it assumes almost constant values at the highest collision energies, independently of v i and j i value. This trend can be seen in Figure 4, where the scattering probability for ortho and para nitrogen molecules, approaching along the normal to the surface plane, at the two considered surface temperatures, is reported.

The Dynamics of N2 Interacting with Quartz
The interaction of N2(vi,ji) molecules in the ground vibrational (vi = 0) and rotational (ji = 0 for ortho and ji = 1 for para symmetries) states with the quartz surface leads to inelastic scattering, whose probability of occurring rises sharply as Ecoll increases and then it assumes almost constant values at the highest collision energies, independently of vi and ji value. This trend can be seen in Figure 4, where the scattering probability for ortho and para nitrogen molecules, approaching along the normal to the surface plane, at the two considered surface temperatures, is reported. For the three lowest collision energy values (Ecoll ≤ 0.1 eV), the integration of most of the trajectories is stopped as the molecule reaches the edge of the crystal lattice by diffusing on the surface through single or multi-rebounds. These trajectories are those for which it is not possible to assign with certainty the final reaction channel, as they are neither completely adsorbed nor completely scattered off in the gaseous-phase. Instead, for Ecoll ≥  For the three lowest collision energy values (E coll ≤ 0.1 eV), the integration of most of the trajectories is stopped as the molecule reaches the edge of the crystal lattice by diffusing on the surface through single or multi-rebounds. These trajectories are those for which it is not possible to assign with certainty the final reaction channel, as they are neither completely adsorbed nor completely scattered off in the gaseous-phase. Instead, for E coll ≥ 1.0 eV the possibility for molecular adsorption also opens up, so that the N 2 molecule remains on the surface at the end of the trajectory. For N 2 /quartz interaction, unlike cristobalite, even at extremely low collisional energies, the second most likely process is the inelastic scattering, albeit with a low probability value. In the complete range of considered E coll values, the initial vibrational state is preserved after the scattering from the surface.
In Figure 5 the final rotational distributions, obtained for the scattering of N 2 (0,0) molecules, are reported for T S = 100 K. It appears that for E coll = 0.06 eV the distribution has a sharp peak in the range between j f = 1 and j f = 4. The distributions widen and the peak height decreases as the collision energy increases. As for the lowest collisional energy value considered, E coll = 0.01 eV, although the number of trajectories is too small to be statistically significant, it is observed that almost all trajectories end in j f = 1.
Molecules 2022, 27, x FOR PEER REVIEW 1.0 eV the possibility for molecular adsorption also opens up, so that the N2 molecu mains on the surface at the end of the trajectory. For N2/quartz interaction, unlike c balite, even at extremely low collisional energies, the second most likely process inelastic scattering, albeit with a low probability value. In the complete range of co ered Ecoll values, the initial vibrational state is preserved after the scattering from th face.
In Figure 5 the final rotational distributions, obtained for the scattering of N molecules, are reported for TS = 100 K. It appears that for Ecoll = 0.06 eV the distributio a sharp peak in the range between jf = 1 and jf = 4. The distributions widen and the height decreases as the collision energy increases. As for the lowest collisional e value considered, Ecoll = 0.01 eV, although the number of trajectories is too small statistically significant, it is observed that almost all trajectories end in jf = 1. Increasing the value of ji to 1, the obtained final rotational distributions, show Figure 6, are very similar to those obtained for the lowest ji value, only showing a propensity to increase by 1 the value of jf corresponding to the maximum of the dis tion in comparison with ji = 0 case. This propensity is also confirmed at the lowest col energy, despite the uncertainty about the probability values corresponding to this en due to the low trajectories number ending with the scattering.
These findings differ from those obtained for N2 molecules' interaction with a g ite surface [13]. Firstly, the scattering of N2 molecules from a graphite surface takes s icant values even at the lowest collision energy (Ecoll = 0.01 eV) allowing us to defin final rotational distributions, while for the quartz, at that energy, the scattering even very few. Secondly, the scattering of N2 from a quartz surface produces final rota distributions with the peak in jf = ji + 2 for Ecoll ≤ 0.3 eV, while for Ecoll > 0.3 eV the jf  Increasing the value of j i to 1, the obtained final rotational distributions, shown in Figure 6, are very similar to those obtained for the lowest j i value, only showing a slight propensity to increase by 1 the value of j f corresponding to the maximum of the distribution in comparison with j i = 0 case. This propensity is also confirmed at the lowest collision energy, despite the uncertainty about the probability values corresponding to this energy, due to the low trajectories number ending with the scattering.
These findings differ from those obtained for N 2 molecules' interaction with a graphite surface [13]. Firstly, the scattering of N 2 molecules from a graphite surface takes significant values even at the lowest collision energy (E coll = 0.01 eV) allowing us to define the final rotational distributions, while for the quartz, at that energy, the scattering events are very few. Secondly, the scattering of N 2 from a quartz surface produces final rotational distributions with the peak in j f = j i + 2 for E coll ≤ 0.3 eV, while for E coll > 0.3 eV the j f value corresponding to the distribution peak is independent of j i , appearing more closely related to the value of collision energy. Finally, no secondary maxima are observed in the region of high j f , bringing out the absence, for the interaction under study, of rotational rainbow scattering, already highlighted for the scattering of N 2 from metal [31][32][33] and graphite [13] surfaces.
corresponding to the distribution peak is independent of ji, appearing more closely re to the value of collision energy. Finally, no secondary maxima are observed in the re of high jf, bringing out the absence, for the interaction under study, of rotational rain scattering, already highlighted for the scattering of N2 from metal [31][32][33] and grap [13] surfaces. Increasing the values of vi to low-medium-lying values, vi = 5 and vi = 10, the int tion dynamics do not change as well as the rotational distribution features of scatt molecules, suggesting a direct correlation between the initial rotational state and the rotational distributions of the scattered molecules as already observed for N2 mole interacting on graphite [13].
To better understand the microscopic mechanism underlying the interaction w determines the results presented above, we have examined in detail some prototy trajectories. From this analysis, it emerges that the interaction occurs mainly throu direct mechanism, in the sense that the molecule is immediately backscattered in the phase, although an indirect mechanism is also observed for some trajectories, most the lowest collision energies or when the molecule impinges on the void space betw pairs of Si atoms on the surface. Under such conditions, the molecule remains on the face from a few tens to a few hundreds of femtoseconds before being backscattere relationship exists between the motion of the molecule after the interaction with the face and the rotational state in which it is reflected in the gas-phase, as found for th scattering from a graphite surface [13]. In fact, molecules leaving the surface moving a cartwheel-like motion are scattered in high jf, while molecules moving with a helico like motion have a low jf. Trajectories are also observed for which the internal motio the molecule changes as an effect of the interaction with the surface so that N2 inci Increasing the values of v i to low-medium-lying values, v i = 5 and v i = 10, the interaction dynamics do not change as well as the rotational distribution features of scattered molecules, suggesting a direct correlation between the initial rotational state and the final rotational distributions of the scattered molecules as already observed for N 2 molecules interacting on graphite [13].
To better understand the microscopic mechanism underlying the interaction which determines the results presented above, we have examined in detail some prototypical trajectories. From this analysis, it emerges that the interaction occurs mainly through a direct mechanism, in the sense that the molecule is immediately backscattered in the gasphase, although an indirect mechanism is also observed for some trajectories, mostly at the lowest collision energies or when the molecule impinges on the void space between pairs of Si atoms on the surface. Under such conditions, the molecule remains on the surface from a few tens to a few hundreds of femtoseconds before being backscattered. A relationship exists between the motion of the molecule after the interaction with the surface and the rotational state in which it is reflected in the gas-phase, as found for the N 2 scattering from a graphite surface [13]. In fact, molecules leaving the surface moving with a cartwheel-like motion are scattered in high j f , while molecules moving with a helicopter-like motion have a low j f . Trajectories are also observed for which the internal motion of the molecule changes as an effect of the interaction with the surface so that N 2 incident with a helicopter-like motion leaves the surface moving with a cartwheel-like motion. In Figure 7, an example of this trajectory kind is shown. Here, we report the time evolution along a typical trajectory of Z (the distance along the normal direction to the surface plane of both N atoms in the molecule), of the rotational number j, of the translational kinetic energy of the molecule center of mass (E CM ), of the energy exchanged with the surface phonons (∆E ph ) and of the effective potential (V eff ), describing the perturbation of impinging molecule to the motion of surface atoms. The N 2 molecule impinging in the lowest ground-state is inelastically reflected in the gas-phase in j f = 17. Figure 7, an example of this trajectory kind is shown. Here, we report the time evo along a typical trajectory of Z (the distance along the normal direction to the surface of both N atoms in the molecule), of the rotational number j, of the translational k energy of the molecule center of mass (ECM), of the energy exchanged with the su phonons (ΔEph) and of the effective potential (Veff), describing the perturbation of im ing molecule to the motion of surface atoms. The N2 molecule impinging in the l ground-state is inelastically reflected in the gas-phase in jf = 17. The translation-rotation coupling is effective when the molecule comes close surface. For N2 molecules interacting with the quartz surface, such coupling has also highlighted in ref. [19] by using a different model potential. Therefore, we can claim the translation-rotation energy exchange is strictly related to the long-range interac governing the scattering of molecules weakly interacting with surfaces, through sisorption states.

Effect of Surface Polymorph: N2/Cristobalite and N2/Quartz
The effect of the surface polymorph on N2 molecule interaction is highlight comparing the results presented in the previous paragraph with those obtained fo tobalite in ref. [29]. Firstly, it is found that when N2 molecules impinge on a cristo silica surface, the dominant process, regardless of the initial internal state and of th lision energy, is the molecular adsorption process, whereas the inelastic scattering h extremely low probability of occurring. Instead, as shown in the previous Section, The translation-rotation coupling is effective when the molecule comes close to the surface. For N 2 molecules interacting with the quartz surface, such coupling has also been highlighted in ref. [19] by using a different model potential. Therefore, we can claim that the translation-rotation energy exchange is strictly related to the long-range interactions, governing the scattering of molecules weakly interacting with surfaces, through physisorption states.

Effect of Surface Polymorph: N 2 /Cristobalite and N 2 /Quartz
The effect of the surface polymorph on N 2 molecule interaction is highlighted by comparing the results presented in the previous paragraph with those obtained for cristobalite in ref. [29]. Firstly, it is found that when N 2 molecules impinge on a cristobalite silica surface, the dominant process, regardless of the initial internal state and of the collision energy, is the molecular adsorption process, whereas the inelastic scattering has an extremely low probability of occurring. Instead, as shown in the previous Section, when N 2 molecules collide with a quartz silica surface the inelastic scattering plays a relevant role, becoming the dominant process at the intermediate and high values of probed E coll . Different findings resulted from the investigation of N 2 interaction on a graphite surface [13] in which molecular adsorption occurs only at the lowest collision energies. Therefore, it appears that both the chemical composition and the crystallographic structure of the surface material have a great influence on the dynamics of a given interacting species. The only difference between quartz and cristobalite is the way of binding the tetrahedral SiO 4 unit in the 3D structure of the surface, while both have a different composition with respect to graphite. Therefore, while the cristobalite exhibits a free region on the surface in which molecules can be trapped, the quartz, due to the displacement of the SiO 4 tetrahedra within different layers, has a more compact structure, thus limiting the chances for the molecule to precipitate below the first layer (see Figure 8). [13] in which molecular adsorption occurs only at the lowest collision energies. Therefore, it appears that both the chemical composition and the crystallographic structure of the surface material have a great influence on the dynamics of a given interacting species. The only difference between quartz and cristobalite is the way of binding the tetrahedral SiO4 unit in the 3D structure of the surface, while both have a different composition with respect to graphite. Therefore, while the cristobalite exhibits a free region on the surface in which molecules can be trapped, the quartz, due to the displacement of the SiO4 tetrahedra within different layers, has a more compact structure, thus limiting the chances for the molecule to precipitate below the first layer (see Figure 8).  [34], β -cristobalite surfaces [35].
The distribution for the graphite surface [36] is also reported for comparison, being useful in the discussion. A break is inserted on the frequency axis between 20 cm −1 and 100 cm −1 . (right) Top view of crystal structures of the corresponding material.
These structural and chemical differences lead, as a consequence, to a different distribution of the phonon states, explicitly considered in our calculations, and reported in Figure 8. It emerges that cristobalite exhibits a band at extremely low frequencies (<20 cm −1 ) that is missing in the spectra of quartz and graphite. Indeed, the first significant band in the graphite phonon state density is centered around the 20 THz (~667 cm −1 ): this spectrum region corresponds to the maximum and intermediate region of the quartz and cristobalite spectrum, respectively. Therefore, the silica surfaces exhibit a greater density of states at low frequencies than graphite. This circumstance is of primary importance for gas-surface heterogeneous reactions, meaning that this frequency region corresponds to the phonons associated with the surface atoms. Thus, the different behavior of quartz is These structural and chemical differences lead, as a consequence, to a different distribution of the phonon states, explicitly considered in our calculations, and reported in Figure 8. It emerges that cristobalite exhibits a band at extremely low frequencies (<20 cm −1 ) that is missing in the spectra of quartz and graphite. Indeed, the first significant band in the graphite phonon state density is centered around the 20 THz (~667 cm −1 ): this spectrum region corresponds to the maximum and intermediate region of the quartz and cristobalite spectrum, respectively. Therefore, the silica surfaces exhibit a greater density of states at low frequencies than graphite. This circumstance is of primary importance for gas-surface heterogeneous reactions, meaning that this frequency region corresponds to the phonons associated with the surface atoms. Thus, the different behavior of quartz is attributable to the different structures of the phonon spectrum. In fact, N 2 molecules interacting with cristobalite undergo a strong coupling with the surface contrary to what occurs on quartz where the coupling is much lower and the molecules are more efficiently reflected.
The influence of the phonon dynamics manifests not only in the occurring surface processes but also in the rotational distributions of the N 2 molecules inelastically scattered by the two silica polymorphs considered here. In Figure 9, the determined distributions for quartz are reported in comparison with those recently obtained for cristobalite [29]. As mentioned in the previous paragraph, looking at Figure 9, it appears that at low E coll the rotational distributions of molecules scattered by quartz exhibit relative peak heights with a maximum for j f < 5 and then gradually decline by extending towards higher j f values as E coll increases. On the contrary, those for cristobalite have a more symmetrical trend around a maximum, that moves to higher values of j f as E coll increases. The distributions here determined for N 2 (v i , j i ) scattered from quartz are very similar to those obtained in ref. [13] for N 2 (v i , j i ) molecules scattered from graphite, although in the case of quartz the secondary peaks at high j f values, ascribed to rotational alignment and rotational rainbow already discussed in refs. [12,13], are missing. by the two silica polymorphs considered here. In Figure 9, the determined distributions for quartz are reported in comparison with those recently obtained for cristobalite [29]. As mentioned in the previous paragraph, looking at Figure 9, it appears that at low Ecoll the rotational distributions of molecules scattered by quartz exhibit relative peak heights with a maximum for jf < 5 and then gradually decline by extending towards higher jf values as Ecoll increases. On the contrary, those for cristobalite have a more symmetrical trend around a maximum, that moves to higher values of jf as Ecoll increases. The distributions here determined for N2(vi,ji) scattered from quartz are very similar to those obtained in ref. [13] for N2(vi,ji) molecules scattered from graphite, although in the case of quartz the secondary peaks at high jf values, ascribed to rotational alignment and rotational rainbow already discussed in refs. [12,13], are missing.  Despite the emergence of differences, both in the probabilities for the various surface processes and the final rotational distributions, the mechanism underlying collision events is the same for the three surfaces, essentially based on the exchange between rotational and translational molecular internal degrees of freedom, as N2(vi, ji) approaches the surface. How efficient this exchange is at promoting molecular adsorption or inelastic scattering depends on the topology and energy put into play by the phonons participating in the interaction. The effect of the polymorph on the rotational inelasticity becomes evident in the final rotation distributions that are generally more specific to cristobalite, whereas for quartz they extend to higher values of jf.  Despite the emergence of differences, both in the probabilities for the various surface processes and the final rotational distributions, the mechanism underlying collision events is the same for the three surfaces, essentially based on the exchange between rotational and translational molecular internal degrees of freedom, as N 2 (v i , j i ) approaches the surface. How efficient this exchange is at promoting molecular adsorption or inelastic scattering depends on the topology and energy put into play by the phonons participating in the interaction. The effect of the polymorph on the rotational inelasticity becomes evident in the final rotation distributions that are generally more specific to cristobalite, whereas for quartz they extend to higher values of j f .

N 2 /Quartz
Looking at Figure 4 in Section 2.3, the interaction dynamics of N 2 molecules in the two roto-vibrational ground states appears independent of surface temperature. A slight surface temperature effect for ortho and para N 2 molecules in the ground vibrational state is evident only for E coll ≤ 0.1 eV and manifests with an increase in the population probabilities for the final rotational states placed to the right of the main peak, which roughly corresponds to the same value of j f for the two surface temperature values. The trend of rotational distributions is preserved and no temperature effect is observed by increasing the initial vibrational state.
A careful examination of the trajectories enables us to ascribe the increment in the population of higher rotational states to the energy contribution of surface phonons, depending on the T S value. Indeed, the same trajectory integrated at the two considered T S values exhibits the main difference in the potential term, V eff , as can be inferred by looking at Figure 10. Here, for the same trajectory, propagated at the two considered surface temperatures, the evolution time of the rotational number and V eff are reported. N 2 (0,1) is scattered in N 2 (0,4) for T S = 100 K and in N 2 (0,5) for T S = 1000 K. Again, from Figure 10, it emerges that the j increment is attributable to the positive and larger contribution of V eff with respect to what happens at the lowest temperature for which this term is mainly negative. In particular, the net increment/decrement of V eff occurring around 1800 fs and 2800 fs correspond to the instant in which the molecule comes close and leaves the surface. Instead, in the range 1800-2800 fs, V eff remains almost constant and the difference in the values assumed by this latter at the two temperatures is almost imperceptible. This trend is due to the fact that in this time-lapse the molecule comes close to the surface and changes its motion, beginning to rotate. Then, the molecule moves away from the surface by moving with a cartwheel-like motion, different from the one with which it had impacted. Instead, in the range 1800-2800 fs, Veff remains almost constant and the difference in the values assumed by this latter at the two temperatures is almost imperceptible. This trend is due to the fact that in this time-lapse the molecule comes close to the surface and changes its motion, beginning to rotate. Then, the molecule moves away from the surface by moving with a cartwheel-like motion, different from the one with which it had impacted.

N2/Cristobalite
In Figure 11, the results obtained for inelastic scattering of ortho and para N2 molecules in the different initial vibrational states from the cristobalite surface, taken at TS = 1000 K, as a function of Ecoll are displayed in comparison with those obtained in ref. [29] for TS = 100 K.

N 2 /Cristobalite
In Figure 11, the results obtained for inelastic scattering of ortho and para N 2 molecules in the different initial vibrational states from the cristobalite surface, taken at T S = 1000 K, as a function of E coll are displayed in comparison with those obtained in ref. [29] for T S = 100 K. The temperature effect is evident in Figure 11 for ortho and para nitrogen ground vibrational state. However, even at the highest surface temperature, the dif in the scattering probability for the two species is perceptible only for the lowest ered collision energy values. Instead, increasing vi values, this effect is lost. The temperature effect is observed only for vi = 0 because molecules, although mostly ing the same trajectories, at the lowest TS values undergo a more efficient interactio the surface. Therefore, some trajectories ending with scattering at the highest TS v TS = 100 K, evolve by molecular diffusion on the surface, terminating when the m reaches the lattice edges or adsorbs. Furthermore, since the interaction is stronge lowest temperature, the molecules show a greater energy exchange which causes tribution of energy at the play, leading to a decrease in the CM translational ener an increase in the molecule rotational energy content. These findings highlight onc the role of translation-rotation coupling.
On the other hand, a nitrogen molecule initially in vi = 5 or vi = 10 holds an content sufficient to counteract the interaction with the surface, even at the lowe perature. Accordingly, the molecule is scattered inelastically with the same probab that at the highest temperature. The higher scattering probability for vi = 0, for values, can be directly related to the corresponding lower energy content respons molecule trapping in the physisorption well from which, through rebounds on the sive potential wall, it is re-emitted. Instead, increasing vi to 5 or 10 the highest av energy determines that the molecule does not remain trapped in the well and ca much closer to the surface. It follows that a large number of trajectories end with m lar adsorption, perhaps also due to the increased effect of the electrostatic potentia P scatt E coll (eV) E coll (eV) Figure 11. Probability (P scatt ) for inelastic scattering of N 2 (v i ,0) and N 2 (v i ,1) from cristobalite for T S = 100 K (blue solid dots) [29] and T S = 1000 K (red solid dots) as a function of the collision energy.
The temperature effect is evident in Figure 11 for ortho and para nitrogen in the ground vibrational state. However, even at the highest surface temperature, the difference in the scattering probability for the two species is perceptible only for the lowest considered collision energy values. Instead, increasing v i values, this effect is lost. The surface temperature effect is observed only for v i = 0 because molecules, although mostly following the same trajectories, at the lowest T S values undergo a more efficient interaction with the surface. Therefore, some trajectories ending with scattering at the highest T S value, at T S = 100 K, evolve by molecular diffusion on the surface, terminating when the molecule reaches the lattice edges or adsorbs. Furthermore, since the interaction is stronger at the lowest temperature, the molecules show a greater energy exchange which causes a redistribution of energy at the play, leading to a decrease in the CM translational energy and an increase in the molecule rotational energy content. These findings highlight once again the role of translation-rotation coupling.
On the other hand, a nitrogen molecule initially in v i = 5 or v i = 10 holds an energy content sufficient to counteract the interaction with the surface, even at the lowest temperature. Accordingly, the molecule is scattered inelastically with the same probability as that at the highest temperature. The higher scattering probability for v i = 0, for both T S values, can be directly related to the corresponding lower energy content responsible for molecule trapping in the physisorption well from which, through rebounds on the repulsive potential wall, it is re-emitted. Instead, increasing v i to 5 or 10 the highest available energy determines that the molecule does not remain trapped in the well and can come much closer to the surface. It follows that a large number of trajectories end with molecular adsorption, perhaps also due to the increased effect of the electrostatic potential.

Computational Setup
The computational setup, used to describe the collisions of N atoms/N 2 molecules with silica surfaces, is based on a semiclassical collisional approach [23]. The method has already been used in the past [12,13,30,34,[37][38][39] to study molecule formation and inelastic scattering on surfaces. More recently, it has been adopted to describe the collision dynamics of N and N 2 with a silica (cristobalite) surface [29]. However, a brief outline of the most important and specific features of the method is here provided.
The dynamics of gas-phase particles is followed by solving the relevant 3D Hamilton's equations of motion self-consistently with the dynamics of the lattice phonons.
The complete Hamiltonian for the two gas-phase atoms interacting with the surface is given by: where Pi is the momentum of atom i having mass mi, V(r) is the interatomic potential of the gas-phase species, ∆E ph is the energy exchanged with the surface phonons. V eff (t, T S ) is the effective potential of the mean field type which depends upon the interaction time and the surface temperature. According to the Ehrenfest theorem, V eff is obtained as the expectation value of the interaction potential between the atoms in the gas-phase and the surface (V int ) on the phonon wave function, that is: R being the two-body distance between each gas-phase atom and the lattice atoms. Ψ (t,T S ) is the total wave function of the phonon modes given as a distribution of quantum states at a fixed surface temperature T S : The time evolution of the phonon wave function is given by analytically solving the time-dependent Schrödinger equations of motion under the harmonic oscillator approximation, that is assuming that a set of M = 3N at − 6(N at is the number of atoms in the 3D slab) independent harmonic oscillators are perturbed by a linear (and quadratic) force exerted between the atom/molecule approaching the surface from the gas-phase and the solid substrate. Thus, by expanding the interaction potential in the normal mode coordinates up to the linear term, the following expression is obtained for the effective potential of Equation (2): where V 0 is the "static" interaction potential between the atoms in the gas-phase and the lattice atoms in their equilibrium positions. V (1) is the linear driving force exerted on each k-th phonon mode Q k . η k (t) are the "phonon excitation strengths" given in terms of the Fourier components I i,k (t) of the external force: With Θ k (t, ω k ) ≈ ω k t, ω k being the frequency of the k-th phonon mode. ∆E ± k ω k , ρ ± k ; T S is the energy exchanged between the chemical particles in the gas-phase and the solid substrate due to the phonon creation (∆E + k ) and phonon annihilation (∆E − k ) processes in the k-th phonon mode.
The energy exchanged with the phonons can be obtained from the transition probabilities P n 0 k →n k for the excitation/de-excitation phonon processes and is given by: where E n k − E n 0 k is the energy exchanged in the transition n k ← n 0 k between the quantum state n 0 k and n k of the k-th phonon mode, and p n 0 k is the Boltzmann distribution of the phonon energies.
Finally, the dynamics of the nuclear motions of each gas-phase atom interacting with the surface is obtained by integrating the classical Hamiltonian equations of motion: Then, the complete dynamical simulation is worked out through three main steps: (1) A 3D slab of considered material is built up from crystallographic data or from ab initio calculations and the corresponding phonon dynamics is determined. (2) An accurate interaction potential for the [gas-phase species/surface] system is obtained as sum of pair-wise interaction starting from ab initio calculations. (3) The dynamics is carried out by self-consistently solving Hamilton's equations of motion of the two interacting atoms [Equation (8)] and the dynamics of the phonons, i.e., by computing the phonon excitation strengths given in Equation (7) at each time step of the classical trajectory.
The surface models, adopted here for the two silica polymorphs, are those of refs. [34,35] for quartz and cristobalite, respectively. The phonons dynamics for these 3D surface models has been determined there.
Recently, the need to include a proper formulation of the long-range attraction in the gas-surface interaction has been demonstrated, since it is promoting the formation of the precursor state that controls all basic elementary processes occurring at the gas-surface interphase [10,[12][13][14]24,25]. Therefore, in the present study, we have evaluated the longrange interaction potentials, represented by the improved Lennard-Jones ILJ model [26], for N/N 2 impinging on quartz, inserted in the new PES formulation (see next Section 3.2), whilst the PES for the interaction of N/N 2 with cristobalite was lately determined in ref. [29]. The adopted method treats the molecular translational motion classically, whereas the rotational and vibrational degrees of freedom of scattered molecules (the latter, represented as Morse oscillators [40]) are analyzed in terms of the action-angle variables by using the semiclassical quantization rules [41]. Accordingly, at each time, the internal energy of the interacting molecule is partitioned into vibrational and rotational contributions by using equations that include rotational angular momentum and vibrational spectroscopic constants. We are unable to predict some features due to quantum effects (tunnelling) and selection rules (forbidden transitions between the different roto-vibrational levels); however, at least in the choice of the initial molecular states, we distinguished between ortho and para nitrogen.

PES Formulation Details
As for N,N 2 /cristobalite [29], the multidimensional PES describing the interaction N,N 2 /quartz has been obtained by combining the interaction potential of the molecule, V N 2 -quartz , with that of the atom, V N-quartz , through a weight function, f sw (r N−N ), that permits a smooth transition from N 2 -quartz to N-quartz interaction potential as r N−N , the intramolecular distance, increases. In particular, f sw (r N−N ) provides quick switching between 1 and 0, preserving the values of V N 2 -quartz and V N-quartz for r N−N lower and higher, respectively, of the intramolecular distance (r diss = 3.5 Å) considered critical for N 2 molecule dissociation. In the range of distances where the switch occurs, the potential is obtained as a combination of both components, weighted according to the f sw (r N−N ) value. To avoid discontinuity points in the transition from N 2 -quartz to N-quartz interaction representation, a continuity check of the obtained potential function has been carried out.
The adopted PES is represented as a sum of pair-wise interactions between the free atom or each effective atom (i.e., in a state different from the isolated one) in the molecule and each effective atom of the crystal lattice. In particular, it is formulated as: where V(r N−N , R) is the total interaction potential. In the first term, on the right side of Equation (1), the variable R is the distance between the i-th atom in the lattice and the k-th N atom in the molecule, while, in the second term, it represents the distance between the i-th atom in the lattice and the gaseous N atom. The analytical expression of f sw (r N−N ) is the same used in ref. [13]. According to the ILJ formulation, each (atom, effective atom)−(effective atom) pair interaction contribution, between N or N(in N 2 ) and each atom in the lattice, is obtained as: The first term of Equation (10) describes the size repulsion, whereas the second one the dispersion attraction associated with each interacting pair. The parameters ε and R m , which represent the potential well depth and its location, respectively, associated with each considered pair, define the strength of both terms in Equation (10). An extended discussion on their choice is presented in the next Section 3.3. Moreover, β is an additional parameter depending on the "hardness" of the two partners. For the neutral-neutral interactions, as the present ones, β has been fixed to 7.5 and m = 6 has been used. For each interacting pair, the ILJ function provides an asymptotic dispersion attraction contribution, related to a partial C 6 coefficient defined as C 6 = ε·R m 6 . The combination of all partial C 6 coefficients determines the value of the global atom/molecule-surface attraction coefficient C 3 [42] which controls capture efficiency and formation of the precursor state of gas-surface elementary processes.

Polarizability and Potential Parameters
Electronic polarizability is considered the key property controlling strength and range of the fundamental components of non-covalent long-range interactions. Basic correlation formulas between polarizability and potential parameters, corresponding to the potential well ε the minimum position R m and the long-range attraction coefficient C 6 , have been proposed some years ago (see ref. [43] and references therein) and exploited here. Since the polarizability distribution in the SiO 2 unit was not known, estimated effective-atomic components of polarizability have been initially assumed. In more detail, the starting values of the ILJ potential parameters, namely ε, R m and related C 6 attraction coefficient (see the previous Section 3.2), were obtained by partitioning the polarizability of SiO 2 into the two contributions due to O and Si atoms in the monomer unit. Then, the polarizability contributions were modified, leaving that of the O atom unchanged and decreasing that of Si atom more uncertain due to the polar bond. We obtained a sequence of ε values in agreement with those reported in the literature, where it is found that the most important centers of interaction are located on the O atoms [21,22,44]. Therefore, the value of ε for the interaction with each O atom must be greater than that with the Si atom. This assumption is also confirmed in ref. [45], where the effective polarizability of Si in SiO 2 is smaller than that of O due to the polarity of the bonds. Accordingly, the used final values for the N-SiO 2 interaction were selected by comparing the value of the so-derived dispersion coefficient C 6 with that describing the attraction of the single monomer unit. The adopted ε and R m parameters, equal to 7.30 meV and 3.67 Å for N-O (in SiO 2 ) and 4.86 meV and 3.67 Å for N-Si (in SiO 2 ), give partial C 6 , respectively equal to 1.78 × 10 4 meVÅ 6 and 1.18 × 10 4 meVÅ 6 . The global (N-SiO 2 ) C 6 value of 4.74 × 10 4 meVÅ 6 for the single monomer unit is in good agreement with the value of 4.49 × 10 4 meVÅ 6 , that results by combining the C 6 for the single coefficient contributions reported in ref. [20], equal to 0.59 × 10 4 meV Å 6 for the N-Si interaction (Si in SiO 2 ) and 1.95 × 10 4 meV Å 6 for the N-O interaction (O in SiO 2 ), respectively.
The same method has been adopted to obtain the parameters for the interaction of N 2 molecule with the SiO 2 unit. A partial value of C 6 = 1.48 × 10 4 meV Å 6 , with R m = 3.62 Å and ε = 6.59 meV, is estimated by us for the interaction of N (in N 2 )-O (in SiO 2 ). Previously [20], for N 2 -O (in SiO 2 ) it was found a global value of C 6 = 2.97 × 10 4 meV Å 6 with R m = 3.79 Å and ε = 10.02 meV (if, as indicated in ref. [20], R m = 3.66 Å is assumed, ε = 12.35 meV would be obtained). Similarly, for N (in N 2 )-Si (in SiO 2 ) we obtain partial C 6 = 0.99 × 10 4 meV Å 6 with R m = 3.62 Å and ε = 4.38 meV. Therefore, the global C 6 for each N 2 -SiO 2 (unit) pair proposed by us amounts to 7.90 × 10 4 meV Å 6 However, results in the literature do not show uniqueness in the reciprocal relationship of the interaction components between N and the Si/O atoms in the lattice crystal. We determined the parameters by assuming a ε value for the interaction with the oxygen atoms greater than that due to the interaction with silicon atoms, but not as drastically as reported in ref. [20]. In addition, for the global interaction of the SiO 2 monomer unit with each N atom, a value of C 6 ≈ 6.00 × 10 4 meV Å 6 is reported in ref. [19], obtained by considering a value of ε for the interaction of N with O atoms (in SiO 2 ) equal to about half of that for the interaction with Si atom (in SiO 2 ). However, in this latter work, the assumed potential is a modified Morse function that works well to describe van der Waals interactions in the region of minima but not at long-range distances. In fact, the MSV (Morse-Spline-van der Waals) model [46,47] was used for a long time to partially overcome these limits, obviated with the formulation of the ILJ potential [26].
Furthermore, the Lorentz and Berthelot combination rules, often exploited to obtain the parameters for pair interactions should be used with great caution as they are effective only for very similar pairs of partners interacting through pure van der Waals forces. The hierarchy of values of ε for the interactions with the two elements in the crystal structure is respected in ref. [48], although perhaps too unbalanced in the interaction with the oxygen atoms, and this leads to a partial C 6 value of 7.75 × 10 4 meV Å 6 for N-O (in SiO 2 ). Instead, in ref. [49], the hierarchy for the values of ε is not respected and the parameters of the LJ potential have been defined using the Lorentz and Berthelot combination rules. In this case, the C 6 coefficient for the N atom interacting with the SiO 2 unit amounts to 3.41 × 10 4 meV Å 6 .
It should also be noted that when the traditional LJ potential is used for an interacting pair, if the depths of the physisorption wells are described correctly, the obtained C 6 is double the correct value. Instead, if the long-range dispersion coefficient is the correct one, the value of the depth of the potential well must necessarily be either underestimated or moved to shorter distances from the surface. This inconsistency results from the fact that the relation for C 6 = 2· ε ·R m 6 , adopted by the classical LJ (12,6) model, is not correct, while in the ILJ model the proper relation is given above. A confirmation of these claims comes from the fact that in most of the considered cases the C 6 obtained using the parameters of LJ potential are higher than the correct ones determined by us and those reported in ref. [20].

Molecular Dynamics Simulations
MD simulations, as described in Section 3.1, involve the integration of the Hamiltonian equations (Equation (8)) describing the interaction for fixed initial conditions and for each interacting species. For each considered set of initial roto-vibrational states and collision energy, we propagate and examine 30,000 trajectories. The molecule CM impinge with polar angle θ = 0 • (defining the selected normal approach) and azimuthal angle (ϕ) of the molecular axis randomly determined for each trajectory. The surface temperature (T S ) is taken at two different values T S = 100 K and T S = 1000 K. The initial coordinates of both atoms in the molecule are chosen randomly in an aiming area on the surface. This area is large enough and at the same time such as to prevent the edge effects during the trajectory propagation. The integration step is 0.20 fs and the accuracy required in the integration procedure is 10 −8 . In addition to the roto-vibrational ground state, we considered the medium excited vibrational states v i = 5 and v i = 10 for the two ground rotational states j i = 0, 1. The CM collision energy (E coll ) has been spanned on a wide range of values, in the range [0.01-2.0] eV.
We can describe and follow the different elementary surface processes occurring when the N 2 molecule impinges on the surface (i.e., molecular scattering, molecular adsorption, dissociative adsorption), assumed as the PES reactive. The impact of a molecule on a surface can give rise to the adsorption of both atoms, to the desorption of just one atom with the other adsorbed on the surface, or to the desorption of both atoms.

Conclusions
The dynamics of collision between nitrogen molecules and a quartz surface has been studied by MD simulations adopting a new determined PES accounting for a suitable treatment for long-range interactions. Impinging N 2 molecules are mainly inelastically scattered in the same vibrational state and in rotational states producing distributions exhibiting sharp peaks for j f < 5 at low collision energies. Instead, the peaks attenuate themselves and the distributions widen to higher values of j f , increasing E coll . The comparison of these results with those reported in ref. [29] for another silica (cristobalite) surface highlights the role of surface polymorph in molecular scattering. Indeed, for silica (cristobalite), the leading surface process is molecular adsorption, whereas for silica (quartz) the inelastic scattering plays a relevant role, becoming the dominant process at intermediate and high collision energy. This finding represents an effect of the different surface's atom packing.
The effect of surface temperature, evident only for silica (cristobalite) and low collision energies, is closely connected with the dynamics of surface phonons. In fact, the spectrum of the density of phonon frequencies only in the case of cristobalite has a band at very low frequencies, those precisely attributable to surface atoms.
For all considered surfaces, the coupling between the translational and rotational molecular internal degrees of freedom, which manifests itself by producing distributions pumped on higher levels than the initial one, plays a primary role in the scattering process.
Author Contributions: Conceptualization, methodology, formal analysis, investigation, writingoriginal draft preparation, writing-review and editing, M.R.; Conceptualization, methodology, investigation, writing-original draft preparation, writing-review and editing, F.P. All authors have read and agreed to the published version of the manuscript.