Lutetium(iii) aqua ion: On the dynamical structure of the heaviest lanthanoid hydration complex

The structure and dynamics of the lutetium(iii) ion in aqueous solution have been investigated by means of a polarizable force field molecular dynamics (MD). An 8-fold square antiprism (SAP) geometry has been found to be the dominant configuration of the lutetium(iii) aqua ion. Nevertheless, a low percentage of 9-fold complexes arranged in a tricapped trigonal prism (TTP) geometry has been also detected. Dynamic properties have been explored by carrying out six independent MD simulations for each of four different temperatures: 277 K, 298 K, 423 K, 632 K. The mean residence time of water molecules in the first hydration shell at room temperature has been found to increase as compared to the central elements of the lanthanoid series in agreement with previous experimental findings. Water exchange kinetic rate constants at each temperature and activation parameters of the process have been determined from the MD simulations. The obtained structural and dynamical results suggest that the water exchange process for the lutetium(iii) aqua ion proceeds with an associative mechanism, in which the SAP hydration complex undergoes temporary structural changes passing through a 9-fold TTP intermediate. Such results are consistent with the water exchange mechanism proposed for heavy lanthanoid atoms.


I. INTRODUCTION
The heaviest element of the lanthanoid (Ln) series is lutetium, a metal whose upper crust abundance of 0.32 ppm is low among rare elements. 1 Due to its electronic configuration, Lu() has the smallest ionic radius, and thus the highest charge density, of all the Ln() ions. 2 Known application of lutetium, as its oxyorthosilicates doped with cerium, is in detectors for positron emission tomography, 3 and it has been proposed its usage in high refractive index immersion lithography as lens material 4 and in petroleum cracking as a catalyst. 5 Another application of lutetium, although only for its radioactive isotope 176 Lu, is in the determination of meteorites age, due to its half-life of 3.78 · 10 10 yr. 6,7 An important applicative reason for the study of Ln() hydration properties is the chemical analogy with the radioactively toxic actinoide series: [8][9][10][11] a thorough knowledge of such properties may be important to improve the efficiency of extracting procedures to separate lanthanoids from actinoids, and even to separate one particular lanthanoid from the other elements of the series.
To this aim, several studies regarding Ln() hydration properties by X-ray diffraction (XRD) or extended X-ray absorption fine structure (EXAFS) spectroscopy are reported in the literature. 12,13 While initially the Ln() series was believed to manifest a sudden change in hydration number from 9 to 8, an hypothesis known as "gadolinium break," [14][15][16] such hypothesis has been disproved by further studies. 13 The actual behavior is a smooth change in hydration number, with most elements of the series having non-integer numbers. A non-integer hydration number represents the concept that in a disordered system, such as aqueous solutions, the solute coordination properties can be better described as an ensemble of different hydration complexes weighted on their respective occurrence probabilities. 19 Many of the reported investigations agree on the mean coordination distance between the Ln() ions and the water oxygen atoms, while regarding the geometry of the hydration complexes results vary among studies. 13 Some works found the associated polyhedron to be a tri-or bicapped trigonal prism (TTP/BTP) for the whole Ln() series, [20][21][22] whereas other works found the heavier lanthanoids to have a square antiprism (SAP) geometry. [23][24][25][26][27][28] In particular, from the recent quantum mechanical charge field (QMCF) molecular dynamics (MD) simulation by Hitzenberger et al. 28 the Lu() aqua ion was found to be octahydrated with a solvation structure that continually changes between SAP and BTP geometries, with the former favored on average. Note that the QM approach imposes a limited simulation time hampering a proper statistical sampling of the hydration structures and therefore an assignment to a particular polyhedron is not possible.
As concerns the water exchange process for Ln() aqua ions, Helm and Merbach proposed a mechanism by which for light lanthanoid ions the solvent exchange is a dissociative process, passing through a SAP transition state, while for heavy lanthanoids it is an associative process with a 9-fold TTP transition state. 12 Instead, the work from Hitzenberger et al. 28 suggests for the Lu() ion the possibility of a dissociative water exchange mechanism passing through a seven-fold transition intermediate. However, this mechanism was proposed on the basis of a single dissociative event taking place in a 20 ps QMCF MD simulation. Clearly, in order to have a statistically significant determination of the water exchange process a much longer simulation time is needed (in the ns time scale). In this contest, classical MD simulations, even if based on semiempirical potentials, are the only suitable approaches to properly describe the water exchange process for highly charged cations. Therefore, more scientific effort is needed to prove, or disprove, the proposed models for both light and heavy lanthanoids also due to the lack of experimental methods that can accurately measure the water-exchange rate constant.
The aim of this work is to provide a conclusive description of the hydration properties of Ln() ions taking into account the dynamical character of the solvation complexes. In particular, we have focused on the study of the structural and dynamical properties of the lutetium aqua ion, being it the heaviest element of the series. In this work we used an approach based on polarizable force field MD simulations, where the theoretical description has been compared with the EXAFS data of a Lu() aqueous solution. The inclusion of explicit polarization in the theoretical framework has proven to give an improved description of the system when highly charged ions are involved, 13,29 while comparison with EXAFS experimental data allows one to validate the accuracy of the theoretical description of the system obtained from the MD simulation. In the present study several MD simulations have been carried out at four different temperatures (277 K, 298 K, 423 K, 632 K) with the aim of providing a reliable determination of the water exchange rate and activation parameters of the Lu() ion in water. Moreover, an accurate analysis of the geometry of the Lu() hydration complexes at room temperature has been carried out based on the MD results.

A. Molecular dynamics simulations
We performed, at each of the four temperatures (277 K, 298 K, 423 K, 632 K), six independent MD simulations for the Lu() ion in aqueous media using the polarizable potential developed by us. 29 In this approach the non-electrostatic ion-water parameters are dependent on the ionic radius. The total interaction potential is given by where V elec is the electrostatic term comprised of a Coulomb and a polarization term following the Thole's induced dipole model, 30 V LJ O−O is the Lennard-Jones potential describing the oxygen-oxygen interaction, and V Lu−O accounts for the non-electrostatic lutetium-oxygen interaction. For waterwater interaction potential we used the TIP3P model with polarization included. For lutetium-oxygen interactions we employed the following pair potential: where the parameters are obtained using D'Angelo ionic radius for lutetium. 2 MD simulations have been carried out in the NVE ensemble with an in-house developed software, using the extended Langrangian method to obtain atomic induced dipoles during the dynamics. The system is composed by 1 Lu() ion and 216 rigid water molecules in a 18.64 Å edge cubic box at room temperature (298 K). For other temperatures, we have kept fixed the number of water molecules and modified the box size in order to have the corresponding densities, 31 as reported in previous simulations, 32,33 i.e., 1.000, 0.9169, and 0.589 kg/dm 3 for 277, 423, and 632 K, respectively. Periodic boundary conditions have been applied to the simulation box. Long-range interactions have been calculated by using smooth particle mesh Ewald method. Simulations have been performed using a velocity-Verlet-based multiple time scale with a time step of 1 fs. Each of the six simulations was equilibrated for 2 ps, followed by production runs of 3 ns.

B. Computational analyses
All computational analyses have been performed by means of in-house written codes. The structural properties have been calculated starting from a total trajectory of 18 ns that has been obtained by concatenating the six independent trajectories at room temperature (298 K). The radial distribution functions for Lu-O and Lu-H pairs are defined by the following equation: where N t is the number of time steps of the MD simulation, ρ is the average number density of oxygen or hydrogen atoms in the system, and V (r) is the volume unit at distance r. The g(r) for first shell O-O pairs has also been calculated following Eq. (3), however the summations run only on the first hydration shell oxygen atoms where C N is the coordination number of the Lu() ion.
where r i j is the Lu-O* distance, θ j ik is the O*-Lu-O angle, and V (r, θ) is the volume unit at distance r and angle θ. The water mean residence time (MRT) was determined employing the method proposed by Impey, 34 which is based on the calculation of the following correlation function: where N t is the number of time steps of the MD simulation and P j runs on the water molecules j being either 0 or 1. It takes the value 1 if molecule j lies within the first coordination shell of the ion at both time steps t, and t + t n , and in the interim does not leave the coordination shell for any continuous period longer than t * . Under all other circumstances, it takes the value 0. This function is then fitted with a single exponential function from which we obtain the MRT as the characteristic time, τ. Water exchange kinetic rate constants (k T ex ) have been determined as the reciprocal of the average MRTs.
By fitting the rate constant at different temperatures on the Eyring equation activation enthalpy, ∆H ‡ , and activation entropy, ∆S ‡ , of the solvent exchange process have been obtained.

C. X-ray absorption measurements
A 0.2 M Lu() aqueous solution was obtained by dissolving in distilled water the stoichiometric amount of hydrated lutetium trifluoromethanesulfonate [Lu(H 2 O) n ](CF 3 SO 3 ) 3 . In order to avoid hydrolysis, trifluoromethanesulfonic acid was used to acidify the solution to about pH = 1. The spectrum was collected at the Lu K-edge in transmission geometry, on the BM23 beamline at the European Synchrotron Radiation Facility of Grenoble. The storage ring was operating in 16-bunch mode, with a typical current of 80 mA after refill. A Si(511) double-crystal monochromator was employed for the data collection, with the second crystal detuned by 80% for harmonic rejection. The aqueous solution was kept in a cell with Kapton film windows and Teflon spacer of 3 cm.

D. EXAFS data analysis
The main contribution to the EXAFS signal of ionic aqueous solutions is associated with the ion-oxygen single scattering signal. In the case of disordered systems rather than using the usual discrete form of the EXAFS equation, the signal is modeled as a function of the g(r) as (9) where A(k,r) and φ(k,r) are the amplitude and phase functions, respectively, and ρ is the density of the scattering atoms. χ(k) theoretical signals are calculated by introducing in Eq. (9) the model g(r) functions derived from MD simulations. The GNXAS program has been employed for the EXAFS data analysis. 35 Phase shifts have been calculated using muffin-tin (MT) potentials and advanced models for the exchange-correlation self-energy (Hedin-Lundqvist). The values of the MT radii are 0.2, 0.9, and 1.7 Å for hydrogen, oxygen, and lutetium, respectively. Both Lu-O and Lu-H g(r)'s obtained from the simulations have been used to calculate the χ(k) first shell signal, as the hydrogen atoms have been found to provide a detectable contribution to the EXAFS spectra of cations in aqueous solutions. [36][37][38][39][40][41][42][43][44] Two nonstructural parameters have been optimized, namely, E 0 (core ionization threshold) and S 2 0 , and the quality of the fit has been determined by the goodness-of-the-fit parameter. 35

A. Lu() aqua ion structural properties
The hydration shell around the Lu() ion features a sharp and well-defined first peak as evident from the Lu-O and Lu-H g(r)'s shown in Figure 1. The Lu-O g(r) first peak maximum falls at 2.32 Å, in agreement with experimental values from EXAFS spectroscopy (2.29-2.34 Å), 17,23,45,46 and XRD data (2.34-2.35 Å). 47,48 Moreover, a second welldefined peak, with a maximum at 4.5 Å, is visible in the Lu-O g(r) and it corresponds to the second hydration shell. Integration of the first peak results in a first hydration shell CN of 8.1 ± 0.2, in agreement with previous results reported in the literature. 12,29 The Lu-H g(r) shows a maximum at 3.00 Å, a region where the Lu-O distribution is zero, thus indicating that the first shell water molecules are strongly structured by the Lu() ion. It is interesting to compare our results with those obtained from a previous QMCF MD simulation. 28 The Lu() hydration complex structure obtained from this theoretical study is in good agreement with our determination as an 8-fold coordination with Lu-O and Lu-H distances of 2.35 and 3.0 Å, respectively, was found.
In order to gain a more direct evidence of the reliability of the MD structural findings it is useful to compare the computational results with the EXAFS experimental data. To this aim the Lu-O and Lu-H g(r)'s have been used to calculate theoretical two-body χ(k) signals starting from Eq. (9). The total first shell contribution to the EXAFS signal has been compared with the experimental spectrum and during the analysis the structural contribution has been kept fixed while the E 0 and S 2 0 parameters have been optimized (the best fit values are 4.5 ± 0.5 eV above the first inflection point for E 0 , and 1.0 ± 0.1 for S 2 0 ). The minimization procedure has been carried out in the k range 2.4-15.1 Å −1 and the best-fit results are shown in Figure 2. The agreement between the experimental and theoretical curves is very good thus proving the validity of the potential functions used in the simulations. This finding is confirmed by the Fourier transform (FT) moduli of the EXAFS signals shown in the lower panel of Figure 2. The FT's have been computed in the k range 3.7-15.0 Å −1 with no phase shift correction applied. Also in this case there is an excellent agreement between the theoretical and experimental determinations. As previously reported the MD simulations provide a CN = 8.1 ± 0.2 for the Lu() first hydration shell, and this indicates the coexistence of 8-and 9-fold coordination clusters. A deeper insight into the dynamical properties of the first hydration shell can be obtained by defining an instantaneous coordination number, as the number of oxygen atoms at a distance shorter than the Lu-O g(r) first minimum, and analyzing its variation along the simulation. Analysis of the instantaneous CN distribution finds an 8-fold coordination for the 97 ± 0.1% of the configurations explored by the MD simulation, and a 9-fold coordination for the 3 ± 0.1% of the configurations. This finding indicates the existence of a stable 8-fold structure for the Lu() coordination shell, while the absence of a significant population of 7-fold configurations is consistent with an associative pathway for the first shell ligand exchange reaction. The existence of an 8-fold hydration complex for the Lu() ion is well accepted in the literature but it has been reported to have either a SAP [23][24][25][26][27] Table I and these values are also highlighted in Figure 3 by vertical green lines. Inspection of this figure shows that the O-O g(r) of the 8-fold configurations is consistent with a distorted SAP geometry as it contains only three main peaks with maxima and intensities in reasonable agreement with the values expected for an ideal SAP geometry.
As concerns the BTP geometry, the quadrilateral faces of the prism should have a square shape if the polyhedron were regular. However, in many crystal structures the prism is distorted into having rectangular faces, with an height longer than the side of the triangular faces. This kind of BTP geometry has 6 groups of vertex pairs. The prism-prism vertex pairs form the 6 sides of the triangles (l BTP 2 ), the 3 long sides of the rectangles (l BTP 3 ), and the 6 diagonals of the rectangular faces (l BTP 5 ). A capped vertex has 4 shorter distances (l BTP 1 ) with the vertices of the capped rectangular face (8 pairs for the 2 capping vertices), 2 longer distances (l BTP 6 ) with the farthest vertices of the prism (4 pairs), and another long distance (l BTP 4 ) with the other capping vertex (see Figure 4). Through geometric relations all O-O distances for a BTP geometry can be calculated in terms of 3 parameters: the Lu-O distance of the prismatic molecules (r P Lu−O ), the Lu-O distance of the capping molecules (r C Lu−O ), and the ratio a between the height of the prism and the side of the triangular faces. For our calculations we consider a BTP geometry with the same a ratio as the crystal structure of [Lu(H 2 O) 8 ](CF 3 SO 3 ) 3 . 17 The r P Lu−O and r C Lu−O values can be extracted from the simulation as the maxima of two Lu-O g(r)'s for the 6 closest oxygen atoms and for the 7th-8th closest oxygen atoms, respectively. These two g(r)'s are shown in Figure Table I, and their positions are evidenced in Figure 3 by vertical red lines. It can be seen that both the position and the intensity of the O-O g(r) maxima for the 8-fold configurations are not consistent with these values, and the peak at 3.96 Å, highlighted in the inset of Figure 3, is not accounted for.
In order to further characterize the 8-fold configurations we evaluated a CDF distribution between the Lu-O* distances and O*-Lu-O angles, where O and O* are oxygen atoms of the 8 first shell water molecules (see Figure 6). CDFs are bidimensional probability distribution functions that combine a radial pair distribution with an angular distribution. It is important to stress that the CDFs are more informative than angular distribution functions (ADF's) as the latter can be obtained from the CDFs by projection on a vertical plane. The CDFs provide additional information on the distanceangle correlation that is essential to obtain a conclusive determination of the hydration complex geometry. Note that the CDF for a SAP geometry should feature three distinguishable peaks centered at a Lu-O distance of 2.32 Å, as all water molecules should be at about the same distance. A BTP geometry should instead show two wide and less symmetric peaks, shifted towards larger Lu-O distances, due to the contribution of the capping water molecules. The CDF evaluated for the 8-fold configurations is reported in Figure 6   BTP polyhedron the contribution in that angular region comes only from the capping-capping water pairs and hence should fall at longer Lu-O distances. The white spots in Figure 6 mark the positions of the theoretical contributions for the ideal SAP structure, while the black spots represent the BTP geometry. Overall, the CDF seems more consistent with a slightly distorted SAP geometry than a BTP. Regarding the 9-fold complex, the most accepted geometry in the literature is TTP. However, another possible coordination could be a capped square antiprism (CSA), where the oxygen atoms form a gyroelongated square pyramid. In order to determine which of the two structures is dominant, White and black spots mark the theoretical contributions for the SAP and BTP geometries, respectively. In particular, white spots represent the SAP structure that is more consistent with the CDF. with a maximum at 4.45 Å. If compared to the 8-fold O-O g(r) in Figure 3, the asymmetric first peak is shifted towards slightly shorter distances. The O-O distances for the CSA and TTP model clusters are shown as green and red vertical lines in the upper and lower panel of Figure 7, respectively. Position and shape of the two peaks in the O-O g(r)'s are more consistent with the TTP geometry. Note that the asymmetry of the first peak of the O-O g(r) is well accounted for by the two sets of distances at 2.70 and 2.83 Å present in the TTP polyhedron. Conversely, in the SAP geometry there is only one short distance at 2.80 Å that is not very compatible with the shape of the O-O g(r) first peak. Furthermore, in the CSA structure there should be a broad low intensity peak at about 3.96 Å, while the O-O g(r) shows a minimum in this region as highlighted in the inset of Figure 7. The peak is therefore significantly shifted towards longer distances for lower angles, but not for angles greater that 73 • , resulting in an L-shaped peak (the L shape is even clearer for the white spots). Such explanation is not possible for a CSA geometry where the contribution of a single capping molecule at longer distance is outweighed by the eight molecules at shorter distance. Finally, the typical peak of SAP and CSA geometries at about r Lu−O = 2.32 Å and θ O−Lu−O = 118.5 • is absent. The visible shift towards longer distances of all contributions in that angular region (115 • -125 • ) is instead typical of BTP and TTP structures, due to the fact that those contributions come only from capping molecules. The accurate CDF analysis therefore confirms that the 9-fold hydration complex has a TTP geometry. Taken together these results lead us to conclude that the Lu() ion hydration complex has a slightly distorted SAP structure and during the dynamics it evolves in a 9-fold TTP intermediate complex.

B. Lu() aqua ion dynamical properties
In order to study in detail the water exchange dynamics of the Lu() aqua ion, we characterized the system at four different temperatures (277 K, 298 K, 423 K, and 632 K) performing, at every temperature, six independent MD simulations of 3 ns each. From each simulation we have evaluated the water MRT, by using the Impey method as   explained in the Sec. II, 34 and the solvent exchange kinetic rate constant k T ex . As it is well known, the MRT values depend on the t * values used. We report in Table III the MRT and kinetic rate constant obtained for two values of t * , 1.0 ps and 0.1 ps from the six simulations at room temperature (298 K), while in Table IV we show the average values (with the associated standard deviation) obtained at different temperatures. The relatively strong dependence of the MRT values on t * is observed when recrossing occurs, as discussed in details by Laage and Hynes. 49 By fitting the values of the kinetic rate constant obtained at different temperatures using Eq. (8), we obtained the activation enthalpy ∆H ‡ and activation entropy ∆S ‡ of the water exchange process of Lu() aqua ion. The results of fitting procedure carried out for t * = 1.0 ps and t * = 0.1 ps are shown in Figure 9, while the resulting ∆H ‡ and ∆S ‡ are reported in Table IV Our results suggest that the water exchange dynamics is a complex phenomenon that occurs on time scales of hundreds of picoseconds at room temperature. In particular, if compared to previous results reported in the literature, the water exchange mechanism for the Lu() aqua ion shows a rate constant k 298 ex lower than those of central elements of the lanthanoid series, such as gadolinium, but similar to the early ones, like lanthanum or cerium. 12,13 Moreover, we obtained a different result as compared to what found by Hitzenberger et al. 28 from a QMCF simulation, where a dissociative mechanism was suggested for the water exchange process. However, the authors state that a reliable conclusion could not be drawn due to the short time interval of the simulation (20 ps), during which only one exchange event was observed. In this contest it is important to stress that a major advantage of classical over QM simulations is that in the former case it is possible to extend the simulation time to the ns time scale, thus allowing a proper sampling of the dynamic and structural properties of the system. Therefore, if one is interested in characterizing the water exchange properties of highly charged cations, classical MD simulation is the only affordable approach. Note that use of classical MD allowed us FIG. 8. Combined distribution function between Lu-O* distances and O*-Lu-O angles evaluated for the 9-fold configurations only. All oxygen atoms included in the calculation belong to the Lu() first hydration shell. White and black spots mark the theoretical contributions for the TTP and CSA geometries, respectively. In particular, white spots represent the TTP structure that is more consistent with the CDF. to carry out six long simulations at four different temperatures from which it was possible to calculate the activation enthalpy and entropy. This would not have been possible with a QM approach. Furthermore, we investigated the instantaneous first shell CN as a function of time obtained by the six independent MD simulations carried out at 298 K, which are shown in Figure 10. This analysis shows that the presence of a 7-fold complex is a very rare event and that both the 7-and 9-fold complexes are short-lived configurations for the Lu() aqua ion. As concerns the mechanism of the water exchange process, as previously mentioned it has been proposed to be associative for the Lu() aqua ion, meaning that first a ninth water molecule enters the hydration shell of an 8-fold complex, forming a transient 9-fold complex, and afterwards another water molecule exits the shell restoring an 8-fold complex. 12,29 Since we found that the 8-fold complex has a SAP geometry, while the 9-fold complex has a TTP geometry, an associative exchange process requires a temporary structural reorganization of the first shell water molecules. Our results also show that 8-fold SAP and 9-fold TTP configurations are in thermodynamical equilibrium, although this equilibrium is almost completely in favor of the SAP clusters (97% 8-fold, 3% 9-fold). Such a gap in populations means that, for the Lu() aqua ion, the free energy of a 9-fold TTP complex is significantly higher than that of a 8-fold SAP complex. As the 9-fold TTP configuration can be considered a high energy intermediate in the structural rearrangement, and thus in the water exchange process, the transition state should be similar to it in structure and energy, meaning that the exchange process would have a high ∆G ‡ . The relatively high ∆G ‡ value obtained from our analysis, along with the low k 298 ex and the short lifetime found for 9-fold configurations are very consistent with the associative model proposed.

IV. CONCLUSION
In this work we focused on the study of structural and dynamical properties of the Lu() aqua ion by the use of polarizable forcefield MD simulations. A theoretical EXAFS spectrum has been calculated from the simulations and compared to the experimental one acquired for a 0.2 M aqueous solution of Lu(), obtaining a good agreement. A thorough geometrical analysis of the ion first shell has been carried out on the simulations, showing solid proof of a SAP geometry for the dominant 8-fold configurations and of a TTP geometry for the 9-fold configurations. In particular, the investigation of the ion first hydration shell has been conducted through the use of O-O g(r) and O*-Lu-O CDF, whose results have been compared with the ideal polyhedron associated with the most probable hydration geometries. This approach has proven to give extensive information on the geometry of a particular complex, improving the commonly applied method of analysing separately the ionsolvent g(r) and the solvent-ion-solvent angular distribution function. Furthermore, the time evolution of the Lu() aqua ion first shell coordination number in the simulations has been studied and the water exchange rate constants at several temperatures have been calculated. Finally, by fitting the kinetic rate constants on the Eyring equation, we have extrapolated the activation parameters, ∆H ‡ , ∆S ‡ , and ∆G ‡ , of the water exchange process between the Lu() first and second hydration shells. The results presented in this paper are all consistent with the associative water exchange mechanism proposed by Helm and Merbach for heavy lanthanoids, 12 furtherly proving the reliability of such model in describing the dynamical behaviour of Lu() aqua ions.