Ab initio studies of O2-(H2O)n and O3-(H2O)n anionic molecular clusters, n12

. An ab initio study of gaseous clusters of O − 2 and O − 3 with water is presented. Based on thorough scans of conﬁgurational space, we determine the thermodynamics of cluster growth. The results are in good agreement with benchmark computational methods and existing experimental data. We ﬁnd that anionic O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters are thermally stabilized at typical atmospheric conditions for at least n = 5. The ﬁrst 4 water molecules are strongly bound to the anion due to delocalization of the excess charge while stabilization of more than 4 H 2 O is due to normal hydrogen bonding. Although clustering up to 12 H 2 O, we ﬁnd that the O 2 and O 3 anions retain at least ca. 80 % of the charge and are located at the surface of the cluster.


Introduction
Understanding the properties of clouds is one of the most challenging and important subjects within atmospheric chemistry. Recently, several studies have suggested that cloud formation is correlated to the influx of cosmic rays and that cosmic rays may contribute to the production of especially low clouds (Marsh and Svensmark, 2000;Carslaw et al., 2002;Svensmark et al., 2009). The primary effect of cosmic rays is ionization of the atmo-Correspondence to: N. Bork (nbor@space.dtu.dk) sphere, producing a variety of molecular cations and free electrons. It has therefore been speculated that the cloud forming effect of cosmic rays is arising from the ions produced. This mechanism seems plausible since both the physical and chemical properties of ions are so different from those of uncharged species. See for example reviews by Harrison and Carslaw (2003) and .
A free electron is very reactive and will most likely be accepted by an O 2 or O 3 due to the abundance and electron affinity of O 2 and O 3 . It is well known that small ions have a high affinity for water and that such ions quickly become hydrated under usual atmospheric conditions. The most important atmospheric hydration mechanism is the stepwise condensation of one H 2 O to the cluster, Both charge transfer and oxidation of a third species are then possible. Another likely event is reaction with a cation. For studying these reactions, knowledge of the structures of the most populated O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters is fundamental since any subsequent reactions may be dependent upon the degree of solvation (Kurtén et al., 2007;Bandy and Ianni, 1998;Gross et al., 2008;Nadykto et al., 2006Nadykto et al., , 2008. Currently, the distribution of the most stable and most important clusters remains uncertain and hence the geometries, thermodynamics, and electronic properties are all unknown. From both theory and experiment it is known that Reactions (R1) and (R2) are highly exothermic for small n.
Published by Copernicus Publications on behalf of the European Geosciences Union. It is however inevitable that the thermodynamics will converge towards that of ordinary hydrated water for large n ( G • 298K = −8.6 kJ/mol) (Lide, 1997). This convergence has not yet been demonstrated and is of interest for studying solvated ions in general.
A number of studies of both O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters have previously been published. Experimental studies in the 1970's determined the thermodynamics of cluster growth for both O − 2 and O − 3 , although only including up to 3 H 2 O (Arshadi and Kebarle, 1970;Fehsenfeld and Ferguson, 1974 (Yang and Castleman, 1990;Skalny et al., 2008;Hvelplund et al., 2010). Also recently, the thermodynamics of 1-5 H 2 O clusters has been studied using ab initio methods (Seta et al., 2003;Lee and Kim, 2002). Although in qualitative agreement with the experimental data, the difference is on the order of 5-10 kJ/mol.
To solve these issues we have undertaken an ab initio study of the initial clustering reactions of O − 2 and O − 3 with water, i.e. Reactions (R1) and (R2), to firmly establish the sizes and properties of these. We have included up to 12 H 2 O which constitute at least two hydration shells. Based on thorough scans of configurational space, we present thermodynamics of Reactions (R1) and (R2) and analyse the resulting structures. Further, we analyse the charge distributions of the clusters and determine the reactivity towards some relevant reactions.

Computational details
Due to the sizes of the clusters combined with the vast amount of possible configurations, optimization of the computational procedure was important. Initially, all structures were optimized using Hartree-Fock theory (HF). Subsequently, the most promising candidates were optimized using density functional theory (DFT).
The anionic charge of the clusters is a known challenge to DFT (Jensen, 2010). Addressing this problem, we utilized the Coulomb-attenuating method via the CAM-B3LYP functional (Yanai et al., 2004). This newly developed functional is based on the well known B3LYP functional, but is incorporating an increasing amount of exact HF exchange at increasing distances. By this, the CAM-B3LYP functional is optimized towards studying charged systems and results are generally much improved compared to standard B3LYP results (Peach et al., 2006;Yanai et al., 2004).
Also concerning the basis set, the negative charge of the clusters is known to require special attention by inclusion of diffuse, long-range functions on all atoms (Jensen, 2010;Kurtén et al., 2008). We utilised the Dunning style basis set aug-cc-pVDZ (Dunning, 1989) since this basis set is small enough to allow calculations on the clusters of interest but large enough to produce accurate results (Seta et al., 2003). However, for scanning configurations at the HF level, the smaller cc-pVDZ basis set was used.
All HF and DFT calculations were performed using the Gaussian 09 package (http://gaussian.com/). Standard convergence criteria were employed.
To confirm the energetics of the CAM-B3LYP calculations, single point calculations using explicitly correlated coupled cluster theory CCSD(T)-F12a  with the VDZ-F12 basis set  were performed. These calculations were performed using the MOLPRO 2010.1 package (http://molpro.net/). The explicitly correlated calculations used density fitting (Manby, 2003;Werner et al., 2007) and resolution of the identity approximations with the default auxiliary basis sets (Weigend, 2002;Weigend et al., 2002;Yousaf and Peterson, 2008). Test calculations on the O − 2 (H 2 O) 1 cluster indicate that the difference between VDZ-F12 and VTZ-F12 binding energies is less than 0.5 kJ/mol.
The CCSD(T)-F12 calculations are based on a restricted open-shell Hartree-Fock reference. Both restricted-open shell and unrestricted coupled-cluster calculations were performed; the difference between these was less than 4 kJ/mol in all cases. Values of the T1 and D1 diagnostics ranged between 0.011-0.039 and 0.03-0.24, respectively. This indicates that some of the studied systems (mainly the O − 2 (H 2 O) 1 cluster, and all the O − 3 (H 2 O) n clusters) have partial multireference character. However, the CCSD(T)-F12a results are still expected to be at least qualitatively reliable.
Since only including up to 5 H 2 O, previous ab initio studies of O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters have not applied systematic scans of configurational space. In the present study the high number of solvent water molecules induce countless possible configurations and thus require a fully systematic approach. Here, we have utilized the simulated annealing method (SA) (Corana et al., 1987) implemented via the Born-Oppenheimer molecular dynamics technique (Li et al., 2000). The main advantage of SA is that the algorithm is not just able to determine the nearest local minima, but also able to migrate from minima to minima. However, even considering mid-sized clusters it is practically impossible to scan all configurations. We chose to terminate the procedure when the same minimum had been confirmed by at least 3 independent calculations for the mid sized clusters (1-7 H 2 O) and 2 independent calculations for the 8-10 H 2 O clusters. For the 11 and 12 H 2 O clusters, the procedure was terminated after 25 SA runs despite no confirmed minimum. Therefore, for clusters containing more than 8 H 2 O but especially for 11 and 12 H 2 O clusters, the probability of undiscovered but more stable configurations is non-negligible.
The optimization were as follows: Initially, a number of randomly generated structures (between 5 and 15) were relaxed at the HF/cc-pVDZ level of theory and used as input for where ω i is the i'th vibrational mode and E F if the Fermi energy parameter. Using a Fermi-Dirac inverse temperature β>0, all frequencies above E F hereby are repressed while all frequencies below E F are enhanced (Santos et al., 2009). We chose E F = 700 cm −1 and β = 0.1 cm, since these values will effectively repress intermolecular vibrations while enhancing modes corresponding to molecular translation. After termination of the SA procedure, the 5-10 most stable structures were re-relaxed using CAM-B3LYP/aug-cc-pVDZ. In no cases were major discrepancies between then most stable HF structure and the most stable DFT structure found.
For evaluating the electrochemical properties of the O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters, knowledge of the distribution of the excess electron is important. The analysis was performed using the charge partitioning method by Bader, also known as the Atoms-In-Molecules approach (Bader, 1998(Bader, , 1990). This method is especially well suited in applications with electronic structure calculations since the required input is the electronic density and nuclear coordinates. The Bader method is based on partitioning the system into atomic regions separated by zero flux surfaces, i.e. surfaces consisting of stationary points with respect to the electronic density. The method, although computationally demanding, is rigid and has been demonstrated to work well, both for  (R1) and (R2). Regimes implying condensation and evaporation are indicated for 298 K and 50 % relative humidity. G ∞ denotes the value of hydrated H 2 O. Black spheres and squares are experimental data from Arshadi and Kebarle (1970) and Fehsenfeld and Ferguson (1974), respectively. charged systems (Bork et al., 2011) and water containing systems (Henkelman et al., 2006). The algorithm by Henkelman et al. was applied (Henkelman et al., 2006;Tang et al., 2009). 1 Since atomic charge is an ill defined property, charge assignation may be ambiguous. We are however interested in the total charge of the O 2 and O 3 species and since the electronic density of the boundaries between the molecules is small compared to the total electronic density, the method is in this case unambiguous (Cochran, 1961).

Thermodynamics
The thermodynamic data of the most stable configurations are summarized in Table 1 and illustrated in Fig. 1.
First, we note that the O − 2 anion is significantly more water affinitive than the O − 3 anion for n≤3. This is likely a consequence of the size of the ions. O − 3 has three atoms to distribute the charge amongst while O − 2 has just two. Consequently, the addition of one H 2 O to O − 2 yields a correspondingly larger decrease of electrostatic energy. However, uptake of the first few H 2 O is significantly exothermic and will occur at all relevant atmospheric conditions for both systems.
Adding more H 2 O, the reactions become smoothly less exothermic. At ca. 6-7 H 2 O, the Gibbs free energy become close to that of hydrated water, G ∞ = −8.6 kJ/mol, but addition of the 8'th H 2 O yield a strong stabilization. At even larger clusters the thermodynamics seems to fluctuate around the value of hydrated water within ca. 5 kJ/mol. Assuming thermal equilibrium, the law of mass action yields where G is the Gibbs free energy change, T is the absolute temperature and R is the gas constant. Activities are here approximated by partial pressures. Equation (2) is for Reaction (R1) but equivalent for Reaction (R2). Assuming a temperature of 298 K and a relative humidity of 50 %, i.e. p(H 2 O) = 16 mbar, yields that G<−10.3 kJ/mol favor condensation, while G>−10.3 kJ/mol favor evaporation. This value is denoted G critical and is illustrated in Fig. 1, separating the evaporation and condensation regimes. Based on Eq.
(2), the concentrations of the investigated clusters were examined. Here, thermal equilibrium is assumed and possible sinks are neglected. The results are shown in Fig. 2. Note the logarithmic scale. For simplicity, the concentrations of the un-hydrated O − 2 and O − 3 ions are set to 1 and all other concentrations evaluated relative to these. Hence, the relative concentrations of the O − 2 based and O − 3 based clusters are here not directly comparable. Although the predicted concentrations are sensitive to temperature and humidity, is is clear that under these conditions clusters of multiple sizes will co-exsist. In stead, many clusters containing more than 5 H 2 O may be found in significant concentrations and these clusters should all be thoroughly considered when assessing the further reactivity of these ions. For an analysis of the cluster size distributions as function of altitude, see Sect. 3.4. Also shown in Fig. 2, on the secondary ordinate, are the accumulated stabilities relative to G critical . These again demonstrate the absence of clusters of particular high stability. Note that the two ordinates are 1:1 linked through Eq. (2).  (R1) and (R2) from various methods. "CCSD(T)" refers to UCCSD(T)-F12/VDZ-F12 calculations on CAM-B3LYP structures. Units are kJ/mol. * and * * indicate data from Arshadi and Kebarle (1970) and Fehsenfeld and Ferguson (1974) respectively. Experimental data up to n = 3 have been produced, based on mass spectrometry (Arshadi and Kebarle, 1970) and the flowing afterglow technique (Fehsenfeld and Ferguson, 1974). These results are also shown in Fig. 1 and Table 2. The comparison with the experimental result for the O − 2 (H 2 O) system is quite poor, differing almost 10 kJ/mol, but the remaining results are in very good agreement with experiment, within 2.5 kJ/mol. Further, the results smoothly converge towards the correct value of hydrated water and hence, both ends of the data series are confirmed. Finally, the absence of "magic numbers", i.e. clusters of particular high stability, is in accordance with the findings of Yang and Castleman (1990) and Skalny et al. (2008). In general, the agreement between these results and previous experimental data is very satisfactory.
Although the CAM-B3LYP functional has been shown to perform well on charged systems both here and previously, DFT is arguably not the ideal method for studying charged systems due to problems with charge delocalization. We therefore confirmed the relative free energies of the smallest clusters using explicitly correlated UCCSD(T)-F12/VDZ-F12 calculations. Due to computational expense, no further structural optimizations were performed and neither was entropy determined. The resulting energies are thus a measure of the ability of CAM-B3LYP to describe the electronic energies of these systems. The resulting energy differences are shown in Table 2.
The agreement between these methods of calculations is in general surprisingly good, with an average deviation of just 1.  with experiment for the O − 2 (H 2 O) system as well, but otherwise, no apparent reason for the increased deviation for these particular systems was found. The remaining systems all deviate less than 1.5 kJ/mol, and especially for the O − 3 (H 2 O) n series, the agreement is excellent. We thus conclude that the obtained results consistently compare well with both experimental and benchmark ab initio data. This is a strong argument in favor of their general quality.

Structures
All structures are illustrated in Figs. 3 and 4.
The O − 2 (H 2 O) n , n≤5, structures have previously been described e.g. by Lee and Kim (2002) and Seta et al. (2003) relying on B3LYP and MP2 calculations with diffuse, double zeta basis sets. No fundamental discrepancies between these previously published structures and the structures presented here could be determined visually.
All O − 2 (H 2 O) n , n≥4 structures adopt a similar configuration around the O − 2 anion where 4 water molecules coordinate to the π* molecular orbital. This is the HOMO orbital facilitating the extra electron (Weber et al., 2000). In the n = 5 and n = 6 structures, this otherwise planar configuration is bent to coordinate the remaining water. For n≥7, the structures become increasingly difficult to describe as many 3, 4 and 5 membered rings are present but the basic n = 4 structure is conserved. Similar to Eq.
(2) we determine the relative populations of two configurations, A and B, by Since we found many structures separated by less than 5-10 kJ/mol we conclude that several structures, other than the configurational ground state, are found at normal atmospheric conditions.  Fig. 3, it is evident that the O − 2 ion does not becomes fully hydrated for n≤12. The O − 2 species remains located inside the cluster but always at or near the surface. The accessibility of the O − 2 species is important for the reactivity of the cluster, e.g. towards charge transfer to O 3 .
Although the structures are given in Fig. 3, more information regarding the hydrogen bonded network can be extracted using radial distribution functions. Here, the centreof-mass distances of O − 2 and H 2 O have been determined and are shown in Fig. 5 as function of n.
Immediately, a clear shell structure is seen. For n≥6, 5 H 2 O are forming a stable first hydration shell, distanced ca. 3Å from the central O − 2 ion. These consist of the 4 H 2 O described in the previous paragraph, coordinated to the π * orbital of the O − 2 species, and a 5'th H 2 O which is positioned above these H 2 O. The second hydration shell is found at distances greater than ca. 4Å and is not fully formed at n = 12. From Figs. 3 and 5, it is evident that growth of the O − 2 (H 2 O) n cluster is mainly by adding H 2 O to the existing structure with minor rearrangement of the existing network of hydrogen bonds.
The O − 3 (H 2 O) n structures have previously been described by Seta et al. (2003) for n≤4 and by Ryabinkin and Novakovskaya (2004) for n = 5. The structures found here are illustrated in Fig. 4. No major discrepancies were found, but generally, the structures are less ordered than the corresponding O − 2 (H 2 O) n structures and visual comparison is difficult. Similar to the O − 2 (H 2 O) n structures, the O − 3 (H 2 O) n structures are characterized by numerous 3, 4 and 5 membered rings and again several almost iso-energetic configurations were found.
Contrary to the O − 2 systems, we note that the O − 3 species is located almost outside the cluster and thus very little shielded by the clustering H 2 O towards external reactants. Also contrary to the O − 2 (H 2 O) n system, no first hydration shell could visually be identified. However, studying the radial distribution functions, illustrated in Fig. 5, clear shell structure was seen. The first hydration shell is found at distances between ca. 3-4Å containing 5 H 2 O. The second shell is found between 4.5 and 5Å containing up to 4 H 2 O. A third shell seems to begin around ca. 7Å.
In the O − 3 (H 2 O) n system, the shells are further from the central ion than in the O − 2 (H 2 O) n system. This is a consequence of both the size of the ion and the electronic structure. Since bond strength (Fig. 1) and bond length generally are proportional, it is not surprising that the O − 2 (H 2 O) n system is more tightly bound. It is also noteworthy that for both Reactions (R1) and (R2), G n G critical for all n≤5, signifying the importance of completing the first solvation shell. Further cluster growth, i.e. adding to the second hydration shells, is associated with significantly less energy gain.

Charge analysis
Studies on electrically neutral clusters of O 2 (H 2 O) n and O 3 (H 2 O) n have been published reporting binding energies much too low for clustering under atmospheric conditions (Tachikawa and Abe, 2005;Loboda and Goncharuk, 2009). It is therefore clear that the strong binding energies for the corresponding charged clusters is a direct consequence of the extra electron. Although formally located at the O 2 and O 3 species, some electronic delocalization occurs and the charge on the O 2 and O 3 species become less negative. We have used the Bader charge partitioning method, briefly described in the "Computational details" section, to analyse the most stable clusters. The Bader charges on the O − 2 and O − 3 molecules at varying degrees of hydration are illustrated in Fig. 6.
For n≤4, each added H 2 O is actively participating in dispersing the charge. For n≥5, significantly less charge is being dispersed upon addition of another H 2 O. This signifies that charge dispersion is the major cause of clustering stabilization up to just 4 H 2 O while normal hydrogen bonding is predominant thereafter. Previous studies of pure water clusters and hydrated electrons show that G is positive for small clusters but becomes negative for clusters containing ca. 4 H 2 O. Addition of further H 2 O leads to increased stabilization which, similar to the present findings quickly converge to values close to those of hydrated water (Lee et al., 2005;Maheshwary et al., 2001). We thereby conclude that the strong stablization of the clusters, up to n = 4, is due to dispersion of the extra electron, while the moderate additional stabilization is due to normal hydrogen bonding as  also found in pure water clusters and for the hydrated electron.
This is again apparent when investigating the HOMO orbital that facilitates the extra electron. As an example is the HOMO orbital of the O − 2 (H 2 O) 8 system illustrated in Fig. 7. As is seen from the spatial distribution of the orbital, only the 4 closest H 2 O molecules contribute significantly while the remaining H 2 O are stabilized by a network of hydrogen bonds.

Tropospheric size distribution
It is well known that for several gas phase reactions, the degree of hydration may be important both for the kinetic and thermodynamic properties (Kurtén et al., 2007;Bandy and Ianni, 1998). Hence, it is of interest to study the distribution of the clusters at varying temperatures and pressures corresponding to typical tropospheric conditions (Mhin et al., 1993). Therefore, we have determined the mole fractions at varying altitudes, assuming a constant lapse rate of 6.5 • C, a ground level temperature of 15 • C and a constant relative humidity (RH) of 25 %. Further, thermal equilibrium is assumed and both sources and sinks are neglected. Thereby, the mole fractions of the various clusters readily can be deduced via Eq. (2). The results are shown in Fig. 8. For clarity, only the 5 most probable clusters are shown.
First, we note that the mole fractions are very sensitive to the amount of water vapor. In more humid atmospheres, the clusters readily grow while smaller clusters are more probable at drier conditions. However, under all conditions will clusters of multiple sizes co-exsist. Also apparent from Fig. 8, is that smaller clusters dominate at lower altitudes due to the higher temperatures favoring water evaporation. However, only under extremely dry conditions (RH <0.1 %) are the distributions dominated by clusters smaller than n = 5.  Secondly, the distribution around the altitude of freezing is interesting, here found at ca. 2.3 km, whereafter gaseous H 2 O must equilibrate with ice rather than liquid water. This effect delays the tendency of cluster growth at increasing altitudes due to decreasing temperatures, and is seen as a kink in Fig. 8.

Simple chemistry
One possible evolution of the O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters is evaporation of neutral O 2 or O 3 resulting in a hydrated electron, i.e. the reaction The thermodynamics of Reactions (R4-R7) is given as the accumulated binding energies of Table 1, the electron affinity of O x , the cohesive energy of neutral water clusters (Lee et al., 2005), and the adiabatic electron affinity of the water cluster (Lee et al., 2003(Lee et al., , 2005, respectively. These values are tabulized in the Supplement. No in-depth discussion will be given, but it is thereby apparent that Reaction (R3) is at least 75 kJ/mol endothermic for O 2 evaporation and at least 210 kJ/mol endothermic for O 3 evaporation. Hence, it is clear that Reaction (R3) is without relevance for atmospheric chemistry.
We also note that previous studies of the vertical electronic detachment energies of these clusters show that this event is negligible under atmospheric conditions (Luong et al., 2001;Novakovskaya and Stepanov, 2004).
Of high relevance, however, is the charge transfer process between O − 2 (H 2 O) n clusters and O 3 . An electron produced by a cosmic ray impact will most likely attach to an O 2 molecule, due to the high atmospheric concentration of O 2 . Since H 2 O is much more abundant than O 3 , the initial O − 2 ion will hydrate to equilibrium before colliding with an O 3 molecule. From the electron affinities it is known that the charge transfer is ca. 160 kJ/mol exothermic for n = 0. From Figs. 1 and 2 it is however apparent that several of the O − 2 (H 2 O) n clusters are significantly more stable than the corresponding O − 3 (H 2 O) n clusters and that the charge transfer process will become correspondingly less exothermic. E.g. is the O − 2 (H 2 O) 1 cluster 13 kJ/mol more stable towards dissociation than the O − 3 (H 2 O) 1 cluster. Consequently is the charge transfer reaction only ca. 147 kJ/mol exothermic. However, evaluating the numbers yields that at no point is the reaction less than 100 kJ/mol exothermic and the Gibbs free energy converge to ca. −110 kJ/mol for large clusters. We thus conclude that transfer of an electron to an O 3 molecule from a O − 2 (H 2 O) n cluster is a likely process for an atmospheric electron generated by a cosmic ray.

Conclusions
An ab initio study of O − 2 (H 2 O) n and O − 3 (H 2 O) n clusters is presented. We have determined the thermodynamics of cluster growth via stepwise addition of H 2 O. The resulting structures show clear shell structure and filling up the first hydration shell, consisting of 5 H 2 O, is highly exothermic. Filling up the second shell is significantly less exothermic, and weather or not it will be filled depends on both temperature and humidity.
The results are in good agreement with experimental data, available for clusters containing up to 3 H 2 O (Arshadi and Kebarle, 1970;Fehsenfeld and Ferguson, 1974). Furthermore, the thermodynamics converges smoothly towards hydrated H 2 O as required, and is thus confirmed for both small and large clusters. The smooth convergence is also in accordance with experiment (Yang and Castleman, 1990;Skalny et al., 2008).
In all systems, the central ion is located at or near the surface of the hydrogen bonded network and the excess electron is in all clusters at least 80 % localized at the O 2 or O 3 species. This signifies that the central ion retains much of its reactivity despite hydration and is easily accessible to further chemical reactions.
Finally, we assessed the reactivity of the clusters towards a few simple reactions. Detachment of an electron or evaporation of a neutral O 2 or O 3 was found negligible, but charge transfer from a O − 2 (H 2 O) n cluster to an O 3 molecule is very exothermic, regardless of hydration level.
We conclude that these clusters are stable under atmospheric conditions, at least up to n = 5 and show no sign of chemical inactivation at increased hydration. These clusters thus remain of interest for ion induced atmospheric chemistry.
In order to investigate larger clusters and their impact on the scattering of electromagnetic radiation from the Sun and the Earth we will in future investigations utilize the methods presented by Poulsen et al. (2001Poulsen et al. ( , 2002 and Jensen et al. (2000).