A thermodynamic view on the microsolvation of ions by rare gas: application to Li + with argon †

We present a thermodynamic perspective of the microsolvation of ions by rare gas atoms, which is based on parallel tempering Monte Carlo (PTMC) simulations. This allows to establish a clear relationship between the structure of the solvation shells and the heat capacity ( C V ) as function of the number of individual solvent species. The dependence of C V on the temperature allows to identify internal structure rearrangements and the onset of partial or total melting of the clusters. As an application, we have employed the PTMC technique to study the thermodynamic properties of clusters resulting from the microsolvation of Li + by argon atoms. Speciﬁcally, calculations have been carried out for the clusters Li + Ar n ( n = 4 − 18 , 33 , 34 , and 38) by applying two diﬀerent potential energy surfaces (PESs): one includes only two-body interactions, while the other also incorporates three-body contributions. Whenever possible, we compare the present thermodynamics results with global optimization studies carried out previously (Prudente, Marques and Pereira, Phys. Chem. Chem. Phys. , 2017, 19 , 25707; Jesus et al. , Int. J. Quantum Chem. , 2019, 119 , e25860). We conclude that the melting process arises for lower temperatures when the model PES accounts for three-body interactions. Additionally, we characterize the melting processes of the ﬁrst and second solvation shells. For some speciﬁc clusters, structural rearrangements of the most external argon atoms are observed at very low temperatures.


Introduction
Solvation plays a fundamental role in a great variety of phenomena within the broad area of physical chemistry.From a theoretical perspective, a popular strategy to look at solvation at the molecular level 1 consists in the step-wise addition of solvent molecules to the solute species 2 and, then, study the properties of the growing cluster.In comparison with continuum solvation approaches [3][4][5][6] , where the solvent is simply described by its dielectric constant, the microsolvation methodology is more computationally demanding, but, on the other hand, allows for a detailed understanding of the phenomena at the molecular level, because all interactions can be explicitly incorporated into the model.
Among the possible properties of the clusters resulting from a microsolvation process, energetics and structure are probably the most calculated ones.Specifically, the development of global optimization methods over the last decades have contributed for an enormous progress in the search of potentially stable chemical structures and, hence, allowed to explore the relevant energy landscapes of microsolvation clusters.By following the global optimization methodology, various research groups [7][8][9][10][11][12] have studied the ion microsolvation by using water and other molecules as solvents.In some works, the infrared spectrum has been also computed for the low-lying cluster structures 11,13 .Moreover, thermodynamic properties can be calculated by applying either Monte Carlo methods or molecular dynamics simulations.The former approach has a lower computational cost, since it only requires the calculation of the energy, while the equations of motion in molecular dynamics need the first derivatives of the PES.Nevertheless, it is worth noting that traditional Monte Carlo method suffers from problems related to slow relaxation and maintenance of ergodicity, which is essential to guarantee an adequate sampling of the phase space.These problems are minimized with the use of extended ensemble techniques, such as parallel tempering Monte Carlo (PTMC), which uses several replicas at different temperatures to build a new ensemble where attempts to exchange between replicas at different temperatures are performed.Although PTMC methodology has been applied to study several atomic and molecular clusters [14][15][16][17][18][19][20][21][22][23][24][25] , to the best of our knowledge, its use with microsolvation systems is very scarce and limited to model solvents 26,27 .
One of the most cumbersome stages of the microsolvation study is the development of the relevant potential energy surface (PES).In general, it implies the calculation of the interaction energy for several geometries at a high-level of theory, followed by a leastsquares fitting to an appropriate analytical function.Because of this difficulty in representing the accurate PES, theoretical insight on microsolvation has been acquired by using single-atom ions and atomic solvents.In particular, cases involving solvation of alkali-metal ions with rare-gas atoms have been subjected to many theoretical and experimental studies [28][29][30][31][32][33][34][35] .These have the advantage to be electronic closed-shell systems, which are easier to treat theoretically.Indeed, the relevant interactions depend essentially on the polarizability of the rare-gas, as well as on the charge and radius of the ion.In this context, we have performed global optimization studies for describing the microsolvation of Li + by argon [36][37][38] and krypton 38 , or a mixture of both 39 .
In this work, we employ our own implementation of the PTMC method to compute thermodynamic properties of the clusters resulting from the microsolvation of the Li + by argon atoms; calculations have been carried out for the Li + Ar n (n = 4 − 18, 33, 34, and 38) clusters.We aim at establishing how the microsolvation reveals itself in terms of the thermodynamic properties.In particular, the main features of the heat capacity curves can be understood as a consequence of the structure and energetics of the clusters.In addition, the influence of the PES, and specifically the importance of including three-body interactions, can be assessed by employing two different analytical potential functions for Li + Ar n clusters that were previously modeled with high-level ab initio energies 36,37 .We show by the present work that results for the heat capacity as a function of temperature present distinct features for the two PESs.Moreover, we are able to characterize how the melting processes of the first and second solvation shells occur for both small and large systems.

Methodology
Within the canonical ensemble, where the number of particles (N), the temperature (T ) and the volume (V ) are constants, the average internal energy of the system may be calculated by the expression where the first (second) term represents the average kinetic (potential) energy; k is the Boltzmann constant.In addition, the calculation of the heat capacity is given by the expression: We should note that our discussion is always performed in terms of the reduced heat capacity quantity, c V , which is defined as c V = C V /Nk.The average terms in Eq. ( 2), i.e., V pot and V 2 pot , are determined by the multi-dimensional integral where R represents a 3N-dimensional vector with the position of the N particles of the system and Z is the configurational integral; note that q = 1 (q = 2) for V pot ( V 2 pot ).The interaction potential used to model the Li + Ar n clusters can be written as 36,37 V where the two first terms refer to the Li + Ar and Ar 2 two-body interactions, respectively, while V Li + Ar 2 and V Ar 3 are the three-body potentials; i-index refers to the Li + ion and indices j, l, and m label three distinct argon atoms.Similar to the original paper 36,37 , here we will use the designation of PES I for the potential with both two-body and three-body interaction terms, while the one that only includes the sum of pair potentials is denoted as PES II.The analytical expressions and the parameters of equation ( 4) are given in Ref. 36.
The integrals to obtain the average values, such as those in Eq. ( 3) can be numerically evaluated through a Monte Carlo simulation by using the Metropolis algorithm: where , is a random walk generated according to the probability distribution It should be mentioned that, although in the present implementation the Metropolis algorithm can be initiated by generating a random configuration, all the Monte Carlo simulations of this work begin at the global minimum structure reported in Refs.36,37.
Specifically, we have developed a computer code that implements a parallel tempering Monte Carlo (PTMC) strategy.The idea behind PTMC is the construction of an ensemble with M different independent simulations (using the Metropolis algorithm), each of them is in thermal contact with a reservoir at a distinct temperature, T m .These simulations are performed simultaneously (in parallel), and attempts to exchange configurations at different temperatures are made periodically according to the following acceptance probability, i ) is the i-th step of the Metropolis random walk associated with temperature T m (T p ).The success of the method mainly depends on the algorithms used to exchange the replicas.The most efficient sampling of phase space is related to the ease of overcoming energy barriers, which can be directly associated with those replica-exchange protocols.
One of the problems arising in Monte Carlo simulations of clusters is the evaporation of the less bounded atoms when the temperature increases.Usually 14 , an attempt to reduce this phenom-ena consists of adding the following rigid-wall function to the interaction potential of Eq. ( 4): where R cm is the vector of the center of mass coordinates and R cut is a cutoff distance; the effect of the magnitude of R cut on the calculated properties is investigated in Section 3.1.Two crucial aspects in PTMC are the choice of the {T m } temperature-set to be used by each independent simulation and, as already mentioned, the exchange process of structures between different temperature simulations.For the first, we have chosen one of the most commonly employed methods, which is based on a geometric distribution algorithm 40,41 .In this method, the temperature schedule is a geometric progression given by with m = 2, • • • , M − 1; T 1 and T M are, respectively, the lowest and highest temperatures considered in the simulation.
To perform the exchange of configurations between Metropolis simulations at different temperatures, we have employed the efficient deterministic even-odd (DEO) algorithm 42,43 .In this scheme, the set of all neighbor-temperature pairs is divided into two subsets: one with all even (T 2 j , T 2 j+1 ) pairs and other with all odd (T 2 j−1 , T 2 j ) pairs; note that j varies between 1 and mod(M, 2), and the exchanges are alternately attempted for the even and odd subsets.
The initial step size is selected by using a geometric progression similar to the choice of temperature [Eq.(9)].However, in order to maintain a 50% acceptance rate for the Monte Carlo (MC) moves carried out at each temperature, we have used the stepsize adjustment proposed in Ref. 44.According to this procedure, the maximum displacement size (∆X(T )) is adjusted after a predefined number of MC steps by using the following expression: In addition to thermodynamic properties, a more detailed analysis of the cluster systems may be achieved by calculating some relevant structural features as a function of temperature.One quantity of interest is the radial probability distribution function (PDF), which is the probability of finding an atom at a distance r from another one.In present work, we define the PDF for the distribution of argon atoms of the cluster in relation to Li + (P Li + Ar (r; T )).The P Li + Ar (r; T ) is normalized, so that the distancedependent coordination number N Li + −Ar (r; T ), obtained by integration of the PDF, converges to the total number of Ar atoms in the cluster.
The Lindemann indices constitute another set of relevant structural quantities that have been considered.Specifically, atom re-solved (δ i ) and the mean ( δ ) Lindemann indices are defined as the root-mean-square bond length fluctuations as follows 45 : These indices are a measure of the fluidity of a given cluster 45 : in general, δ < 0.1 indicates a quite solid cluster, whereas a clear fluid behavior corresponds to δ > 0.2.
In order to identify the main structures arising in each simulation, we have performed an additional PTMC run to store a thousand of configurations for selected temperatures; this is carried out at each 1000 MC steps (or 2000 MC steps for larger clusters) in order to guarantee the storage of non-correlated configurations.Then, we proceed with a local optimization of each configuration, which is thus relaxed to a stable structure within the same energy-landscape basin.Relaxed structures with energies differing by less than 5 × 10 −4 mE h were considered as being equal.

Volume of the available configurational space
Prior to the systematic calculation of the heat capacity for different cluster sizes, we investigate how the cutoff distance of the rigid-wall potential [Eq. ( 8)] influences the values of the relevant quantities.Thus, we have calculated the reduced heat capacity and the mean Lindemann index for Li + Ar n (n = 6 and 12), with various values of R cut ; in this investigation, only the PES II was employed.In Figure 1, c V and δ for both Li + Ar 6 and Li + Ar 12 are represented as a function of temperature for different values of the parameter R cut , which is obtained by the expression: Here D max is the distance of the farthest atom from the centerof-mass of the global minimum structure of the cluster and α varies between 1.1 and 2.0.The results displayed in Figure 1 were obtained from 10 PTMC runs with N total = 10 6 , whereas the temperature-exchange procedure is performed every 25 MC steps; M = 801 temperatures between T 1 = 1 K and T M = 500 K were considered in Eq. ( 9) for these calculations.We can observe in Figure 1 that c V is strongly dependent on the α parameter, which is related to the volume considered by the integration in Eq. (3).Indeed, small values of α lead to c V curves that diverge from those for larger simulation boxes, even at low temperatures (around T ≈ 50 K for Li + Ar 6 and T ≈ 20 K for Li + Ar 12 ).It is also clear that, for both aggregates, simulations with small α values tend to underestimate the mean Lindemann index.This appears to indicate that small simulation boxes do not allow for significant changes in the structure of the cluster and, hence, the calculated thermodynamic properties reflect the rigid-type potential at the walls.
In turn, Figure 1   tural changes are apparent (i.e., for large values of δ ).This happens because larger volumes favor the evaporation of atoms from the cluster when the temperature increases.Thus, we conclude that the choice of α is a delicate process: by one hand, a sufficiently small volume should be used in the integration to prevent evaporation that tends to occur, mainly, at high temperatures; by the other hand, the size of the simulation box should be sufficiently large to reduce, as much as possible, the spurious wall effect on the phase transitions.Actually, we should emphasize that the exact influence on the simulation of finite systems using a spherical rigid-wall potential is not totally understood, even for the well-studied cases of atomic and molecular clusters (see Refs 17,46-48 and references therein).Based on the present results, we have chosen α = 1.5 for the calculations involving the Li + Ar n clusters up to n = 18.However, as D max tends to increase with the number of argon atoms, it is necessary to reduce the α-factor in Eq. ( 13) to avoid an excessive integration volume.Then, we have considered α = 1.4 for the largest clusters, that is, Li + Ar n (n = 33, 34, and 35).This procedure is consistent with the choice of Cezar et al. 24 for large Lennard-Jones clusters.
To further verify the reliability of our choice for D max , we represent in Figure 2 the calculated Li + -Ar PDF (left panels) and the corresponding coordination number (right panels) for both Li + Ar 6 (top panels) and Li + Ar 12 (bottom panels) at some selected temperatures.All curves in Figure 2 were calculated with PES II, but analogous results are expected for PES I.By selecting temperatures in the range that covers the relevant features of c V (including the peaks in Figure 1), we are able to capture the evolution of the PDF and the expected Li + -Ar distances as T increases.For Li + Ar 6 , the PDF shows a very pronounced and straight peak at the lowest temperature, since all the six argon atoms are kept essentially at the same distance from the ion (cf. the corresponding integrated PDF curve).As the temperature increases for T = 295 K (c V peak), however, the PDF curve shows, now, a tail that extends at larger Li + -Ar distances.Thus, the argon atoms are no more "fixed" near the ion (i.e., forming a solvation shell), but rather have the ability to move away from Li + well beyond the equilibrium distance.This indicates that the solvation shell is beginning to "melt".At T = 500 K, the PDF tail significantly increases and the PDF-integrated curve shows that only a reduced number of argon atoms are at Li + -Ar distances compatible with the solvation shell.Indeed, the values of N Li + Ar (r; T ) indicate that approximately one Ar-atom is farther away than R cut = 6.97 a 0 (i.e., a value 150% greater than the equilibrium Li + -Ar distance for the first solvation shell).
In the case of Li + Ar 12 , the PDF has two distinct peaks at T = 39 K, which correspond to the six argon atoms in the first solvation shell (the one at smaller distances) and the other six in the second solvation shell; the corresponding integrated-PDF curve shows that the two sets of six argon atoms are essentially localized at the corresponding equilibrium distances, i.e., the position of the two PDF-peaks.We should note that the peak associated to the second solvation shell is broader than the one for the first shell, which is consistent with a more fluid-like behavior, as we had anticipated in a previous work 39 .At T = 70 K, the smallest peak becomes broader due to a larger mobility of the secondshell atoms, while the straight peak associated to the first shell keeps essentially unchanged.As the temperature reaches 198 K, the second solvation shell is already melted.In turn, the atoms of the first solvation shell increase their mobility and "melt" only at higher temperatures (T ≥ 324 K), which is associated to the second peak of the c V curve in Figure 1.Specifically, we have performed PTMC calculations on the ionic clusters Li + Ar n (n = 4 − 18) as modeled by PES I and PES II.Such calculations were identical to those performed in subsection 3.1, with α = 1.5 for all cases.The standard error of c V was also computed and the corresponding values fall in the range between 0.11% and 0.57% for small clusters (such as Li + Ar 4 with PES I), or between 0.22% to 1.4% for larger clusters (such as Li + Ar 18 with PES II).For the latter, nonetheless, the standard error was less than 0.5% for most temperatures considered in this work.

Li
Results for the reduced heat capacities, c V , are displayed in Figure 3, while the mean Lindemann indices, δ , are shown in Figure 4. To better highlight some details, the c V and δ curves are shown by individual plots for each cluster in the Supplementary Information (see Figures S1-S15).The reduced heat capacity for each cluster size coincides, within the statistical error, for both PES I and PES II at very low temperatures; at T = 1 K, the value of c V increases with increasing cluster size.As a general trend in Figure 3, we observe the appearance of two main peaks (or, just one for n ≤ 6) for c V as a function of T , which is independent of the PES.As discussed for Li + Ar 12 in subsection 3.1, the c V -peak at lower-temperature (higher-temperature) is associated to the melting of the second (first) solvation shell in clusters with n ≥ 7.
It is apparent from Figure 3 and Table 1 that three-body interactions have a relevant effect on the phase transitions of the Li + Ar n clusters: peaks of the c V curves occur at lower temperatures for PES I than for PES II, which is consistent with the less deep poten- tial energy landscape of the former 36,37 .It is also worth noting that, around the phase-transition temperature, the mean Lindemann index tends to show a faster variation with T in the case of PES I than for PES II, which is in agreement with a more straight c V -peak associated to the melting of the second solvation shell of the clusters modeled with PES I (cf. Figure 3).Moreover, the variation of the mean Lindemann index depends on the cluster size, as well as on the potential employed in the simulation: the Li + Ar 4 cluster is essentially rigid ( δ < 0.1), while Li + Ar 5 and Li + Ar 6 are intermediate between rigid and fluid ( δ is slightly above 0.1).However, until completing the first solvation shell (i.e., n ≤ 6) the variation of the Lindemann index upon increasing temperature is relatively small.For n ≥ 7, the mean Lindemann index tends to show only one transition at low temperature, which corresponds to the melting of the second solvation shell.Indeed, the strong increase of the Lindemann index (i.e., reaching values greater than 0.3) can be attributed to the large-amplitude motion of the argon atoms that are farther apart from Li + .This corroborates the idea of a "fluid"-like second solvation-shell, which has been suggested in our previous global optimization study 39 .
Although not always clear from Figure 3 (see Figures S1-S15 for further details), we must point out that a small c V -peak appears for some cluster sizes at low temperatures (up to 5 K for PES I and up to 15 K for PES II); also similar features have been reported for ion-Stockmayer clusters 26 .In the present work, the low-temperature c V -peak can be observed for n = 5, 8 and 9 (n = 5) for PES I (PES II).In addition, there is a "shoulder" in the  same temperature range for n = 10 with both PESs.Such peculiar features of c V deserve additional analysis in order to characterize the physical process associated to this low-temperature regime.
In the case of Li + Ar 5 , the small c V -peak arises at T ≈ 2.8 K (T ≈ 12.9 K) for PES I (PES II).By looking at the corresponding PDFs around those temperatures, which are represented in Figure 5, we conclude that the c V -peaks may be attributed to a large-amplitude motion of the argon atoms.This is similar to an "umbrela" motion that drives the system from the global minimum 36,37 (where four Li + -Ar distances are slightly larger than the other one) to a degenerated and indistinguishable one, i.e., corresponding to the same point in the PES.Indeed, the structures obtained by local geometry-minimization at such low temperatures are all coincident with the global minimum of Li + Ar 5 .In turn, Figure 5 shows that the two maxima of the PDF occurring at T ≈ 2.8 K (T ≈ 12.9 K) for PES I (PES II) coalesce to a single broad-peak at higher temperatures, which is a clear fingerprint of a no more hindered "umbrela" motion.
As for Li + Ar 8 and Li + Ar 9 with PES I, the low-temperature c Vpeaks arise at T ≈ 2.8 K and T ≈ 3.4 K, respectively.Conversely, no peak is observed in the case of PES II at such low temperatures.This difference in the behavior of c V (T ) for PES I and PES II may be attributed to specific features of the low-energy landscapes of the two PESs.Whereas we have detected only the global minimum structures of Li + Ar 8 and Li + Ar 9 during the simulations with PES II at such temperature regime, several other  low-energy local-minima were found when the PES I was employed.In the case of Li + Ar 8 , four low-energy structures are assigned as I (global minimum), and II, III and IV (local minima) in Figure 6; for each structure, bars indicate the corresponding frequency, after geometry relaxation, at different temperatures.We have detected a stable structure for Li + Ar 8 (structure II) whose energy is only 0.02 mE h above the global minimum of PES I.It is important to refer that about 18% of the configurations explored in the simulation at T = 2.8 K can be associated to such low-lying local minimum; at T = 11.0K, structure II is already the most important motif (∼ 44%) present in the simulation, though other two minima (structures III and IV) also appear within more than 1% (cf. Figure 6).In turn, Li + -Ar PDF displayed in Figure 7 (left panel) shows only one peak associated to the argon atoms of the second solvation shell at T = 1.0 K, but it is split into two peaks (i.e., compatible with structure II in Figure 6) at T = 2.8 K, which corresponds to a structural transition.
It is worth noting that a similar situation occurs in Figure 8 for Li + Ar 9 , but in this case the first local minimum (assigned as structure II), which is just above the lowest-energy structure (i.e., ∆E = 0.002 mE h ), becomes the most prevalent structural motif even at T = 1 K (ca.69% against 30% for the global minimum); another local minimum (structure III) appears 1% of the times in the simulation at this temperature.The presence of these lowest-energy minima (structures I, II and III in Figure 8) is apparent in Figure 7 (right panel) by the two PDF-peaks shown at larger Li + -Ar distances, i.e., at 8.7 a 0 and 9.5 a 0 .Accordingly, we should emphasize that the three argon atoms of the second solvation shell of structure I (II) are all located at about 9.5 a 0 (8.65 a 0 ) from Li + , whereas structure III shows two of those Li + -Ar distances with ∼ 8.65 a 0 and the third one with ∼ 9.4 a 0 .Because the PDF represented in Figure 7 (right panel) is an weighted average over the Li + -Ar distances of all structures arising during the simulation, we found a fractional value for the number of argon atoms associated to each one of those peaks: 1.62 (1.38) atoms for the peak at 8.7 a 0 (9.5 a 0 ) when T = 1 K.At the temperature of transition (T = 3.4 K), structure III is already the second most frequent minimum (∼ 23%), behind structure II which has a frequency of 50%.In addition, structures I, IV and V also appear at such temperature, with frequencies of ∼ 15%, ∼ 3% and ∼ 6%, respectively, and the presence of several low-energy minima leads the PDF to spread out across the two large-distance peaks (as observed in the right panel of Figure 7).
In a different way, the Li + Ar 10 cluster in both PESs shows essentially a constant value of c V , rather than a prominent peak, at very low temperatures.In the case of PES I (PES II), c V ≈ 2.87 (c V ≈ 2.95) in the temperature range from 4.0 K (9.0 K) up to 6.4 K (16.0 K).Although the structural "phase transition" pointed out for Li + Ar 8 and Li + Ar 9 is not now evident for Li + Ar 10 , it becomes apparent when monitoring the structures obtained by relaxation of the configurations along the PTMC simulations carried out at those low temperatures.For instance, the calculations with PES I show four low-energy minima for the Li + Ar 10 cluster.At T = 4.4 K, frequencies of 82% and 18% were obtained for the global minimum and the other three local minima, respectively; nonetheless, the contribution of the latter increases to about 30% at T = 6.2 K.
It is important to underline, at this stage, the common feature associated to the "structural transition" occurring in all of the three clusters, which we have just discussed: Li + Ar 8 , Li + Ar 9 and Li + Ar 10 show low-energy minima (other than the global one) that were already reported in our previous global optimization study 39 .Nonetheless, the existence of low-energy minima cannot be seen as a sine qua non condition to observe a c V -peak at low temperatures.Indeed, we have found in Ref. 39 that the cluster Li + Ar 14 modeled with PES I has a local minimum whose energy is slightly above the global minimum (i.e., ∆E = 0.007 mE h ), but no c V -peak is observed in Figure 3 at very low temperatures.In fact, the global (local) minimum shows a frequency of 95% (5%) at T = 5 K and this situation does not change too much as the temperature increases (the frequency of the local minimum becomes 10% at T = 15 K).This result appears to indicate that, in contrast to Li + Ar 8 and Li + Ar 9 , the potential well of the Li + Ar 14 local minimum in PES I may be narrower than the one corresponding to the global minimum.Now, we move the discussion to the higher-temperature c Vpeak in Figure 3, which may be attributed to the melting of the first solvation shell (as already pointed out above).For the smaller clusters, this c V -peak tends to be higher than the one associated to the melting of the second solvation shell, but the situation is reversed for n ≥ 10 with both PESs.Moreover, the c V -peak associated to the melting of the first solvation shell tends to disappear for n ≥ 14 (n ≥ 17) when employing PES I (PES II).This is easily understood if one has in mind that the number of atoms of the first solvation shell becomes very small in relation to the total number of atoms as the size of the cluster increases.However, the "melting" of the first solvation shell may be identified for the larger clusters by representing the Li + -Ar PDF for various temperatures.In the case of PES I, this is shown in Figure 9 for the clusters with a completed first solvation shell at T = 100, 136 and 170 K.It can be observed in Figure 9 that the prominent PDFpeak associated to the first solvation shell extends to large distances as the temperature increases (see the inset of the figure).Clearly, Figure 9 shows that the first PDF-peak tends to collapse at T ≥ 136 K.This indicates that, as the temperature increases, the six argon atoms of the first solvation shell are able to escape from the surroundings of Li + .

Large clusters
To investigate how the thermodynamic features described in the previous section evolve for larger clusters, we have performed PTMC calculations on the Li + Ar n (n = 33, 34 and 38) systems as modeled with both PES I and PES II.We should mention that the aggregates Li + Ar 34 and Li + Ar 38 exhibit high-symmetry global minimum configurations (T and O h point groups, respectively), while the lowest-energy structure of Li + Ar 33 in PES I (PES II) belongs to the C 1 (C s ) point group 37 .For these large systems, Fig. 9 FDPs for Li + Ar n (n = 6−18, 33, 34 and 38) at temperatures around the "melting" temperature of the first solvation shell (as indicated in each panel).All curves are for clusters modeled with PES I.Note that the prominent peak at ∼ 4.9 a 0 refers to the first solvation shell, whereas the one corresponding to the second solvation shell is already collapsed at these temperatures.
12 PTMC runs with N total = 2 × 10 6 have been performed in each case and the temperature-exchange procedure is performed every 50 MC steps.Moreover, we have considered M = 230 (M = 278) temperatures between T 1 = 0.5 K and T M = 250 K for Li + Ar 34 and Li + Ar 38 (Li + Ar 33 ).In these cases, the set of temperatures, {T m }, was chosen from constant steps in predefined intervals, instead of the geometric progression given by Eq. ( 9).This procedure allows to reduce the number of sampled temperatures (hence, reducing the computational effort) and it has been also used in other PTMC studies 27 .All calculations have employed α = 1.4 (see Subsection 3.1).
The reduced heat capacity (c V ) and the mean Lindemann index ( δ ) are displayed in Figure 10; as above mentioned for the smaller clusters, the c V and δ curves are also shown by individual plots for each cluster in Figures S16-S18 of the Supplementary Information.The standard error of c V was also computed and the corresponding values fall in the range between 0.22% and 2.4% for the Li + Ar n (n = 33, 34 and 38) clusters.For the latter, nonetheless, the standard error was less than 1.0% for most of temperatures considered.In all cases, we observe in Figure 10 one prominent c V -peak, which corresponds to the melting of the external solvation shell.As shown in Table 1, the melting temperature is around 40 K (50 K) for PES I (PES II).The difference in the values of T melt for the two PESs may be attributed to distinct welldepths, as above mentioned for lower-size clusters.It is worth noting in Figure 10 that, in general, the mean Lindemann index is similar for both PESs.Just before T melt , δ is around 0.2 for Li + Ar 33 and Li + Ar 34 , while the corresponding value for Li + Ar 38 is less than 0.1.This indicates that the latter is more solid-like than the former ones at such temperatures.
Although the c V curves do not display a maximum (or a shoulder) corresponding to the melting of the first solvation shell, the analysis of the Li + -Ar PDF at T = 136 K (middle panel of Figure 9) shows that the argon atoms next to Li + (i.e., associated to the first PDF peak) may be delocalized up to large distances.As we have observed for lower-size clusters, the temperature for destroying the first solvation shell is high and essentially invariant with the cluster size.Then, it appears that the Li + ion solvated by a firstshell of Ar atoms constitutes an independent sub-cluster from the thermodynamic point of view.
We should also emphasize that the c V values at the lowest temperatures (e.g., T = 1.0 K) are equal, within the statistical error, for PES I and PES II.In contrast to the smaller clusters discussed in previous section, the low-temperature c V values now slightly vary with the size of the aggregate.This may be rationalized from a classical perspective by having in mind that, in the lowtemperature limit, the heat capacity, C V , depends on 3N − 6 (with N = n + 1) vibrational degrees of freedom.Hence, c V is proportional to 3 − 6/(n + 1) and, thus, it becomes less dependent on n as the size of the cluster increases.
In the low-temperature regime for Li + Ar 33 , we observe in Figure 10 a distinct behavior in the c V curves obtained with PES I and PES II.For the latter, there is a small c V -peak at T = 2.6 K, which corresponds to a structural reorganization.This is apparent from Figure 11 which displays the frequency of the main relaxed structural-motifs as a function of temperature for the large clusters.It is clear from this figure that, for Li + Ar 33 with PES I, the global minimum (C 1 structure) is the most frequent up to ∼ 36 K, while there is an exchange in the frequencies of the two lowestminima around T = 3 K for PES II.It is interesting to note that The green curve refers to the global-minimum structure, while the blue one is for the second lowestenergy minimum.In turn, the yellow curve for n = 38 refers to another minimum structure with a frequency-maximum around 7%.Finally, the magenta curve represents the total frequency of the remaining structures (all having small frequencies), i.e., those not included in the other curves.
the second lowest minimum of PES II corresponds to the global minimum structure of PES I. We may conclude that, from a thermodynamics point of view, the C 1 structure is the most stable one up to T ≈ 35 K (between ∼ 3 K and ∼ 44 K) for PES I (PES II).In turn, we may observe in Figure 10 that the mean Lindemann index for Li + Ar 33 with PES II is less than 0.1 in the temperature range where the global minimum is the dominant structure.Thus, it appears that the global minimum is a more solid-like than the C 1 structure.For higher temperatures, there are several structural motifs appearing whose collective frequency becomes larger than the one for the above mentioned C 1 structure.
Conversely, there is not so much difference between PES I and PES II for Li + Ar 34 .We observe in Figure 11 that the global minimum is dominant at low temperatures, but a great diversity of structures emerges at temperatures slightly below the melting temperature, T melt , associated to the c V -peak appearing in Figure 10.In particular, there is a C 1 structure (i.e., the secondlowest minimum) with a maximum frequency reaching about 5% for both PESs.In turn, all other minimum structures show fre-quencies below 3% and, hence, we collect them in a single curve.
Moreover, the intersection between the dominant-structure frequency-curve and the collective-curve arises at approximately the melting temperature for Li + Ar 33 and Li + Ar 34 , whereas such crossing occurs at lower temperatures for Li + Ar 38 , due to the appearance of different structures with a high frequency at lower temperatures (cf. Figure 11).Indeed, the Li + Ar 38 global minimum is dominant at low temperatures, but a C 3v structure has a frequency maximum of ∼ 19% at T ≈ 21 K in PES I (∼ 5% at T ≈ 33 K in PES II).Besides these two low-energy minima, we have observed three higher-energy minima with a maximum frequency greater than 4%; all these minima are low-symmetry structures in both PES and we represent in Figure 11 the corresponding ones with the highest frequency.

Conclusions
We have employed PTMC calculations to describe ion microsolvation in isotropic solvent media.This allowed to characterize the first stages of microsolvation phenomena from a thermodynamic perspective and, clearly, it sheds light over previous optimization studies involving such kind of systems.The most relevant quantities are the reduced heat capacity (c V ), the mean Lindemann index ( δ ), the ion-solvent PDF, and the frequency of the main relaxed structures.Specifically, we have studied the microsolvation of Li + with argon atoms by employing two PESs.One of the PESs describes only two-body interactions (PES II), while the other includes, in addition, the three-body terms of the manybody expansion (PES I).In general, the major difference between the c V results from PES I and PES II are related to the melting temperature, which is expected to fall at lower temperatures for the less-attractive PES I.
In the present thermodynamics work, we are able to characterize the melting temperatures of the first and second solvation shells.The latter arises at low temperatures and, in some cases, it may be preceded by structural transitions corresponding to rearrangements of the external argon atoms at even lower temperatures.These structural rearrangements show that a thermodynamics study is fundamental to assess the most stable structures of the clusters as a function of temperature.In turn, the global optimization analysis is able to identify the energy differences among the stable configurations, but it cannot account to the structural rearrangements.
In contrast, the melting of the first solvation shell occurs at much higher temperature, which corresponds approximately to the T melt value found for the Li + Ar 6 cluster.However, such transition becomes difficult to detect for larger clusters (n ≥ 14 for PES I and n ≥ 17 for PES II) through the observation of the c V curve.In these cases, we have looked at the Li + -Ar PDF for temperatures around the above mentioned T melt value.At such temperature, it is expected that the tail of the PDF-peak associated to the six argon atoms of the Li + neighborhood would extend up to large Li + -Ar distances, thus showing a certain delocalization of the solvent atoms.Accordingly, this results clearly show that the structure formed by Li + and its six nearest-neighbor Ar-atoms can be viewed as an independent sub-aggregate of the growing microsolvation cluster.This picture is compatible with a "rigid"-like

)
Parameters a and b are chosen so that ∆X(T ) is increased (decreased) by a factor in the interval [1, c] ([1/c, 1]).In particular, we have employed the values indicated in Ref. 44: a = 0.672924 and b = 0.0644284 with P ideal = 0.5 for c = 3.
shows that c V increases with the value of the α parameter in the range of temperatures where significant struc-J o u r n a l N a me , [ y e a r ] , [ v o l .] ,

Fig. 2
Fig. 2 PDF (left) and integrated PDF (right) of the Li + Ar 6 (top) and Li + Ar 12 (bottom) clusters for different temperatures related to the main features of the corresponding c V curves.

+
Ar n (n = 4 − 18) clusters We have conducted a detailed thermodynamic study on the initial stages of the microsolvation of the Li + ion by argon atoms.

Fig. 5
Fig. 5 Li + -Ar PDF of the Li + Ar 5 cluster for three low temperatures comprising the small c V -peak in Fig. 3: PES I, left panel; PES II, right panel.

Fig. 7
Fig. 7 Li + -Ar PDF of the Li + Ar 8 (left panel) and Li + Ar 9 (right panel) for temperatures around the c V -peak.

Fig. 10
Fig. 10 Reduced heat capacity and Lindemann index for the Li + Ar n (n = 33, 34 and 38) clusters modeled with PES I (magenta) and PES II (green).

Fig. 11
Fig. 11 Configurations for the Li + Ar n (n = 33, 34 and 38) clusters modeled with PES I (left panel) and PES II (right panel).The green curve refers to the global-minimum structure, while the blue one is for the second lowestenergy minimum.In turn, the yellow curve for n = 38 refers to another minimum structure with a frequency-maximum around 7%.Finally, the magenta curve represents the total frequency of the remaining structures (all having small frequencies), i.e., those not included in the other curves.

Table 1
Melting temperatures (a) (T melt ) for the Li + Ar n (n = 4 − 18, 33, 34 and 38) clusters obtained from the corresponding c V curves.For n ≥ 7, T melt refers to the melting of the second solvation shell.
(a)Temperatures are given in Kelvin.