Size-dependent solvation of p-H(2) in (4)He clusters: a quantum Monte Carlo analysis.

Variable-size (4)He(N) clusters doped with a single p-H(2) molecule are studied here using variational and diffusion Monte Carlo calculations that show the highly quantum nature of the dopant and the solvent. Energetic and structural features extracted from our analysis reveal that the p-H(2) molecule behaves as a gentle perturber: The He droplets remain essentially liquidlike, with no evident structural change with respect to the pure ones. The p-H(2) dopant represents a kind of "competitor" for helium in the smaller droplets because it can replace the solvent adatoms; it also remains immersed in the cluster as N increases although located off-center within the droplet, while, finally, getting fully solvated in the larger droplets. The calculations are carried out up to N=100 as the largest number of solvent (4)He atoms and clearly show no evidence of either shell structuring or of "magic" numbers in the size of the smaller droplets.


I. INTRODUCTION
Nanodroplets of 4 He of various size constitute an ideal spectroscopic matrix because of the very low temperature which can be reached inside them ͑below 0.5 K͒. [1][2][3][4] Helium clusters have several unusual physical properties 5 which distinguish them from all other clusters: A notable exception is provided by molecular hydrogen clusters, which are also surprisingly similar to those with 4 He. Many of the differences are related to the macroscopic quantum properties of bulk helium which include superfluidity. Furthermore, He droplets are able to pick-up any species with which they collide 2,6 and which can then reside in the bulk of the droplet or at its surface depending on the strength of the helium-dopant interaction with respect to the He-He one. Helium nanodroplet isolation 7 is thus especially well suited for the study of a wide variety of molecules ranging in size from a few atoms to large organic molecules of biological interest ͑e.g., amino acids and nucleobases͒. 2 Further advances in the preparation of solvated molecular species 8,9 were able to reveal a host of interesting features connected to the quantum nature of the solvent environment: One typical example is provided by the OCS molecule, for which the infrared spectrum when solvated in 3 He nanodroplets closely resembles the spectra of heavy molecules immersed in liquids, with their unstructured broad shapes, while the corresponding spectrum in a bosonic 4 He nanodroplet is similar to that of the molecule in the gas phase, with its well-defined P-and R-branches. 3,10 This intriguing observation further led to the hypothesis that the dopant molecule undergoes essentially free rotations in the bosonic solvent and could be interpreted as indication of superfluid behavior within a finite system. 11 Several experimental and theoretical works have been also dedicated to the understanding of the properties of hydrogen clusters and of the H 2 bulk. The goal was to verify the possible presence of superfluidity in parahydrogen systems ͑p-H 2 , with total nuclear spin I =0͒ since the possibility that liquid parahydrogen might also exhibit superfluid behavior was raised some time ago. 12 Although the p-H 2 molecule obeys the Bose statistics, and has a lighter mass than 4 He, the low-temperature bulk phase is not a superfluid but a hcp solid. 13 Attempts to produce a superfluid phase by supercooling the normal liquid below the triple point of 13.8 K have been thus far unsuccessful. 14 A different approach was proposed by Grebenev et al. [15][16][17][18] who, in their experiments at Göttingen, were able to introduce a tiny quantity of parahydrogen and orthodeuterium ͑o-D 2 , I Ͼ 0͒ within a helium droplet also containing the OCS molecule. Each of the droplets contained between 14 and 16 parahydrogen molecules and, in the stable state, the OCS molecule is at the center of the droplet surrounded by a thin layer of hydrogen. The latter is, in turn, surrounded by a relatively thick shell of liquid helium-4 to which it was also possible to add an outer shell of liquid helium-3. 15 In the pure 4 He droplets ͑0.38 K͒ both the systems with the p-H 2 and the o-D 2 exhibit spectral features that indicate the excitation of angular momentum around the OCS axis. In the colder 4 He / 3 He droplets ͑0.15 K͒ these features remain in the o-D 2 cluster spectra but disappear in the p-H 2 spectra, where no excitation of the molecular axis angular momentum was detected. These results are consistent with the onset of superfluidity in finite size hydrogen systems. Our own analysis of the structuring of p-H 2 around OCS ͑Ref. 19͒ using diffusion Monte Carlo ͑DMC͒ revealed the clear presence of shell formation around the dopant molecule.
Many theoretical papers have also investigated the nature of hydrogen. Superfluidity properties in pure p-H 2 clusters were studied by path integral Monte Carlo ͑PIMC͒ calculations 13 and by a similar method based on the so-called worm algorithm. 20 In both cases the superfluid fraction is not a trivial function of the size dimension of the system: This finding was referred to as "quantum melting." 20 No superfluid response was obtained for the N = 33 cluster, 13 while the large fluctuations in the superfluid fraction were related to important structural changes, as those which occur from N =25 to N = 26. 20 Other PIMC calculations 21 agree with the experimental findings reported previously, showing anisotropic local superfluidity for the hydrogen layer around the OCS molecule. Barnett 24 a result confirmed by our later computational findings of a bound state for that system. 25 We therefore present below VMC and DMC calculations for the p-H 2 ͑ 4 He͒ N system with N from 1 to 100: For brevity of presentation, in the following discussion 4 He will simply become He and p-H 2 H 2 . In Sec. II a short description of the computational methods and their implementation is reported. Our findings are explained in Secs. III and IV, where our calculations are seen to exclude the possible existence of magic numbers related to those suggested for the pure He clusters. 26,4 Finally, our conclusions will be reported in Sec. IV.

A. VMC
Quantum Monte Carlo ͑QMC͒ methods 27 represent a powerful tool to solve many-body problems thanks to their favorable scaling with the number of particles of the system. VMC and DMC can be considered the simplest approaches based on a random walk procedure which allows one to reliably approximate the exact, but unknown, solutions of the many-body ͑electronic or nuclear͒ Schrödinger equation. In our case, we apply our QMC implementation to extract the energetics and several geometrical information from the theoretical study of the highly delocalized H 2 in He clusters.
An initial variational approach that optimizes the parameters of the chosen trial wave function with respect to the variance of the mean of the local energy E L is combined with a DMC procedure. We have employed the "rigid body" 28 scheme for the treatment of the H 2 molecule whose bond length has been kept frozen at 0.74144 Å. The details of that implementation have been given by earlier publications of our group 29,30 and therefore will only be briefly outlined in this Section. The masses we have used are 2.015 650 070 and 4.002 603 24 for H 2 and 4 He, respectively.
For all the systems we have used 1500 walkers ͑500 for N = 100͒ propagated in blocks of 3000 time steps ͑6000 for N = 2, 1000 for N = 20, 21, 30, 50, and 100, 2000 for N =19 and 22͒, according to the pure Metropolis algorithm; the number of blocks was kept to 10 or 20. The VMC simulation was repeated during 3 ͑2 for N =1͒ optimization cycles: After a random walk, in which the ensemble of walkers spans the configurational space of the system, a minimization technique ͑using the Powell method 31 or a genetic algorithm 32 ͒ is applied to obtain the best set of parameters which minimizes the chosen cost function for that specific VMC distribution. In addition to the variance, the cost function can be the mean of the local energy, its second moment, the mean of the absolute deviation of E L with respect to a reference energy E R and other quantities related to a possible Lorentzian shape of the distribution of E L . However, one standard approach uses the variance of E L , which implies a Gaussian distribution of E L around the "true" value. 33 The trial wave function of a cluster, composed of a homonuclear diatomic impurity and N solvent atoms, is given by a product of nodeless exponential forms, where R i is the distance between the center of mass of the impurity, i is the usual Jacobi angle and one He atom, and R jk represents the He-He pair distance. The pair function used to recover the He-He distribution is described by a typical two-body Jastrow factor, with ͑R͒ = exp͑f͑R͒͒, characterized by a set of variational parameters p i , which are then optimized. For the H 2 -He we have instead used a more complicated exponential which takes into account the anisotropy of the potential energy surface ͑PES͒, 34 where Z = cos and P n are the Legendre polynomials of order n, with n even. The f Ϯn are Jastrow factors with the same functional form as that given by Eq. ͑2͒. In the optimization step we need to span a 20-dimension configurational space in order to obtain the best set of parameters for the more accurate DMC calculations.

B. DMC
The DMC procedure relies on the well-known short-time approximation, 27 in which the Schrödinger equation is solved through a relaxation process in imaginary time where the starting point is generally given by a trial wave function. The diffusion-drift process is simulated by iterating the Langevin equation to the first order: A standard numerical approach, which is not computationally expensive. For the branching part we follow the approach suggested by Assaraf et al. 35,36 in which the number of walkers is maintained constant through a minimal reconfiguration step where we minimize the fluctuation of the weights associated to each walker during the stochastic propagation. At each lth iteration we calculate the average branching factor for the M walkers where the branching factors are calculated as in Ref. 28, i.e., and where E L ͑x i ͒ and E L ͑y i ͒ are the local energies before and after the move, ͗E͘ l−1 is the average energy in the previous ͑l −1͒th iteration and eff is the effective time step which depends on the acceptance/rejection ratio. 37 In the original work 36 a constant value, E ref , close to the eigenvalue, was used instead of ͗E͘ l−1 . We then create a relative weight for each walker, given by The population of M walkers is reconfigured by dividing it in two subsets: There are M + walkers with i,l Ն 1 and M − walkers with i,l Ͻ 1. In order to maintain a constant population, a number of the M + walkers substitute the same number of M − walkers. This number is given by 35 where is a uniform random number in ͓0,1͔.
The expectation value of an observable is then obtained as an average over the ensemble of the configurations. 30 In our calculations the parameters of DMC simulation depend on the size of the system and can be summarized as follows.
• Apart from the largest clusters ͑N = 19, 20, 21, 22, 30, 31, 50, and 100͒ we have employed a "double" Metropolis scheme, in which translational ͑relative to the solvent͒ and rotational ͑relative to the dopant͒ moves are treated separately. Several numerical tests for N = 50 and N = 100 with the "double" random walk gave us essentially the same results: In these cases, the small anisotropy of the dopant is not important.
• Time step of 200 a.u.: This choice guarantees an acceptance of 99% for all the clusters in both moves.
• 10, 20, or 30 blocks have been used according to the size of the system.
• 1000, 3000, 5000, or 9000 time steps for each block; these values represent the best compromise between viable computational effort and acceptable accuracy in the calculations.
We always obtain a small error in the energy, oscillating between 10 −3 and 10 −1 cm −1 ͓the error is of the order of cm −1 only for H 2 ͑He͒ 100 ͔, as seen in Table I. This result is further an important piece of evidence on the high quality of the optimized trial wave functions generated by our VMCoptimization method. The bias on the geometric distributions due to the DMC mixed estimate turned out to be completely negligible when corrected with a forward walking algorithm. The size dependence on allows us to have a systematic 99% acceptance ratio in the Metropolis step of the DMC calculations. The dependence on of the expectation values turns out to be largely negligible for all our systems.

III. INTERACTIONS, STRUCTURES, AND ENERGETICS
The global interaction forces are treated as a sum of two-body potential contributions, which we call the sum-ofpotentials approximation, where the first sum runs over the number of He atoms and takes into account the H 2 -He interaction ͑described in Jacobi coordinates͒, while the second sum runs over the He-He pairs. Neglecting higher-order contributions represents in this case a good approximation: The presence of a gentle perturber as H 2 does not induce any appreciable electric moment for the He atoms, whose electronic density remains essentially unchanged. 38 Structural features and energetics within this potential scheme represent therefore a realistic description of the H 2 ͑He͒ N clusters. The H 2 -He interaction is described by the very accurate ab initio PES from Boothroyd et al. 34 The H 2 molecule is Ϫ4.640͑4͒¯Ϫ7.29͑4͒ 1.55͑5͒ 10 Ϫ5 fixed at its equilibrium distance of 0.74144 Å and the full angular dependence is included. As one can see from Fig. 1, the spatial anisotropy of this interaction is not very strong: Potential cuts at = 0 and = / 2 show substantially the same characteristics, both in well depth and minimum distance. For He-He pairs we have employed a recent version of the Tang-Toennies model. 39 The comparison in Fig. 1 points out the differences between the pairs: H 2 -He is characterized by a slightly deeper well ͑ϳ10 versus ϳ7 cm −1 ͒ and by a larger minimum distance ͑ϳ3.4 Å at = 0 against ϳ3.0 Å for the He-He case͒. This means that the H 2 location in the small He droplets studied here will depend on the balance between the stronger interaction with the solvent and its lighter mass, which is related, in principle, with a larger zero point motion and a greater delocalization of the impurity.
In Fig. 1 we also show the potential used by Barnett and Whaley 22 in their calculations which is the spherical average ͑V 0 ͒ of an anisotropic potential 40 fitted through a Lennard-Jones-type formula. The curve is very similar to the / 2-cut of the full Boothroyd-Martin-Peterson surface, but it underestimates the overall interaction, hence affecting their VMC results. 22 In Fig. 2 we report three quantities which are useful for the understanding of the role of H 2 in He clusters. They are defined in the following way: The DMC energies for the pure He N are taken from Refs. 41 and 42. E ex is the exchange energy, i.e., the amount of energy released from the system when a He atom is replaced by a H 2 molecule. Its values are always negative: H 2 gives more stability to the cluster, even when oscillations occur, indicating that the He replacement with H 2 could be energetically less favorable. We see clearly that the change in E ex with N becomes negligible from N = 10 on, where it assumes the nearly constant value of ϳ−1.25 cm −1 ͑left panel of Fig. 2͒, perhaps hinting at a saturation effect. The binding energy E bind represents the gain in energy for the nanodroplet when a H 2 molecule is attached to it. The general behavior is similar to that of E ex ͑left panel of Fig.  2͒, although shifted by one unit on the abscissas, although the energy values are now slightly larger, as expected.   Finally, we report on the right panel of Fig. 2 the evaporation energy ⌬E, the energy needed to extract one He atom from the cluster as defined by Eq. ͑11͒. No definite shell structure appears in the series, but we see that the bulk regime is not yet reached and ⌬E is always found to be less than 2.5 cm −1 ͑apart the N = 100 case͒. Indeed, applying the simple formula ⌬E mn = ͑E͑H 2 ͑He͒ m ͒ − E͑H 2 ͑He͒ n ͒͒ / ͑n − m͒ it is possible to estimate the evaporation energy for H 2 ͑He͒ 50 and H 2 ͑He͒ 100 which is 2.3͑3͒ and 3͑3͒ cm −1 , respectively, less than the bulk value of about 5 cm −1 . 43 In the N = 100 cluster, the large statistical uncertainty ͑Ϯ3 wavenumbers͒ is due to the intrinsic difficulty to obtain an accurate wave function from the time-consuming minimization procedure. The oscillations of the evaporation energy ⌬E up to N =20 are probably due to a small structural effect which stabilizes preferentially the species with N = 13 and N = 20. For larger clusters the uncertainty of the data does not allow to draw any further conclusion.
In Table I 22 and to our present DMC findings: Keeping in mind the differences in the potential model employed, we find that our values are always lower than those of Ref. 22 and that the expected agreement could be only qualitative because these authors studied the H 2 ͑He͒ N system at the VMC level, which is not expected to be accurate enough for this very delocalized, and weakly interacting, system. As an example, our "exact" discrete variable representation result for H 2 He is −0.038 cm −1 ͑to be compared with −0.036 cm −1 in Ref. 25͒, while our DMC eigenvalue is −0.034͑3͒ cm −1 , i.e., 90% of the exact result and in good agreement with the experimental estimate 24 of the binding energy of about Յ−0.028 cm −1 , which demonstrates the accuracy of the computational methods employed. In the last column of Table I the single atom evaporation energies ⌬E are also reported.
The above analysis on the energetics has shown no dramatic change in the He environment for the presence of H 2 , which therefore can be classified as a "gentle" impurity. Thanks to its light mass and its electronic configuration, a bosonic ͑I =0͒ p-H 2 molecule could replace a He atom within the droplet without strong perturbation of the latter. The mass effect would tend to push the molecule out toward the droplet surface, while the slightly stronger interaction potential would collocate H 2 inside the cluster as a weak coordinator center. Replacement effects are already known from experimental and theoretical works ͑Refs. 16, 21, and 23͒: Our geometrical analysis represents an attempt to clarify the role of the very quantum H 2 particle in the bosonic He environment of a confined liquidlike structure. Figure 3 collects the radial distributions for the H 2 -He pair distance, averaged over the angular variable: The curves are normalized to the number of He atoms in each cluster. Across our series of cluster sizes, their growth leads to curves of very similar shapes, all characterized by a long tail in their distributions which is an indicator of a large delocalization of the light molecular partner. We also note that their maxima are shifted toward shorter distances as N increases: This last finding can be interpreted as a sign of increasing compactness of the system with the addition of solvent at-oms. The findings for the larger clusters are also characterized by the appearance of a shoulder at larger distances: The droplet increases its size to allow for the addition of more He atoms which are not all "geometrically" equivalent. The shoulder becomes a second well-distinguishable maximum in the H 2 ͑He͒ 30 system. The second peak, dominant for N = 50 and N = 100 and located at ϳ8 and ϳ10 Å, respectively, cannot yet be really considered as a probe of a shell closure, as also shown by the earlier calculations reported in Ref. 22, which emphasized the extensive delocalization of H 2 in small helium clusters. In that work, 22 the rms distance of the H 2 dopant with respect to the center of the cluster was always larger than that for the He atoms, for all N values. 22 That work also found a H 2 resistance to reside inside the cluster since it was mostly attached to the surface, although a progressive embedding was seen as the cluster size increased: Our DMC calculations seem to show instead a more pronounced H 2 solvation, as further analyzed below. We report in Fig. 4 the normalized densities in Å −3 for the solvent atoms and for H 2 for the species with N = 19 and N = 50. The former can be compared directly with Fig. 7c of Ref. 22 and we can see that our H 2 is located well inside the cluster. This  effect is further enhanced in the larger cluster as shown in the same figure by the N = 50 profiles. In Fig. 5 two-dimensional density maps for H 2 ͑He͒ 5 , H 2 ͑He͒ 20 , H 2 ͑He͒ 50 , and H 2 ͑He͒ 100 are presented in order to understand the spatial arrangements of the following clusters.
• The He density around the impurity is initially not symmetric with respect to a plane perpendicular to the molecular axis, even if the dopant molecule is homonuclear. Although the H 2 -He pair trial wave function is symmetric by construction ͓see Eq. ͑3͔͒, the resulting total distribution around the molecule can be asymmetric. This effect can be pointed out by carefully choosing an appropriately unambiguous reference frame during the two-dimensional binning procedure. We choose a reference frame in which the z-axis always coincides with the molecular one and we reorient all the walkers in such a way as to have a positive value for the z-coordinate of the cluster center-of-mass.
• Even for N =20 H 2 is located off center and the increased size of the cluster does not show any shell structuring within it.
• Clusters with 50 and 100 He atoms are characterized by a spherical distribution of the solvent around the dopant, now at the center of the droplet. A simple explanation for this effect is provided by considering the similarities in relative strength between the He-H 2 interaction and the He-He network of interactions, and the corresponding differences in reduced masses which cause the zero point energy ͑ZPE͒ to be larger for the solvent atoms with respect to the solute. Hence the H 2 dopant is drawn to the center of the droplet by its larger binding energies while there are also repulsive forces with the solvent network that tend to push it out of the cluster. As a result, the smaller clusters are characterized by an asymmetric distribution which is overcome in the larger ones by the stronger binding forces acting on the dopant: Hence the final prevalence of spherical distributions beyond N = 50.
• All the cases shown here indicate the presence of a small "bubble" of few angstroms in diameter induced by the repulsive forces between the He and H 2 closed shells due to the Pauli exclusion principle, as evidenced in the larger value of the potential turning point for H 2 -He with respect to He-He ͑see Fig. 1͒. The bubble progressively migrates to the center of the droplets as N increases.
The fact that H 2 eventually gets embedded in the He droplet is confirmed by the results summarized by Fig. 6. H 2 -He and He-He average distances are reported on the left panel while H 2 and He average distances with respect to the geometric center ͑GC͒ of the cluster are shown on the right panel. Some specific features can be observed from the following distributions.
• One should note, first of all, that the monomer complex, H 2 -He, exhibits the largest average distance between the partners because of the weakness of the interaction and of the highly quantum character of He-H 2 which is not balanced by further attraction centers.
• The R values for all He-He pairs are the larger ones, confirming the highly quantum nature of this special solvent despite the presence of an impurity.
• In both panels one can see that the results for N = 2 are consistent with a description of it in terms of a very floppy equilateral triangle: In this case the dopant replaces one of the solvent atoms. The same is observed for N = 3 but now it is not possible to extract any "classical" analogy.
• In the right panel of Fig. 6, a clear size expansion of the droplet, with a contemporary increased solvation of H 2 , is shown up to N = 15; beyond this value the oscillatory behavior may be explained in terms of the cluster at- tempt at structural stabilization. The He-GC mean distance for N = 50 is about 7.5 Å while that for H 2 -He assumes the same value it had in the smaller clusters ͑about 2.5 Å͒.
• H 2 is therefore located inside the cluster, as the results reported in the right panel of Fig. 6 clearly emphasize, and resides in a vacuum "bubble" with a diameter of about 6 Å that varies little with cluster size; • on the left panel of Fig. 6 the computed He-He distances in pure He droplets, obtained from the DMC calculations of Ref. 41, are also reported. The comparison with the He-He distances obtained here for our doped clusters clearly indicates that the solvent is scarcely perturbed by the presence of the H 2 dopant and closely resembles its pure liquid situation.
• Calculations in Ref. 22 pointed out an enhanced resistance to H 2 penetration in the 13-He cluster; our results instead indicate no special behavior at this size, the mean distances being very similar in both frames to the ones of the previous and following sizes of clusters.
It is interesting to note here that we loosely introduce the concept of a bubble to describe the structureless void surrounding the dopant molecule bound within the cluster. The physical origins of such empty areas could be related to at least two different features: ͑i͒ The onset of the repulsive walls in the H 2 -He potentials ͑see Fig. 1͒ which occurs around 3 Å, and ͑ii͒ the large ZPE of the light H 2 dopant which places it, within the overall cluster potential, very close to the zero-energy threshold. Thus, it is not surprising to see that our DMC calculations provide average bubble diameter of about 6 Å, in keeping with the above considerations.
Finally, Fig. 7 represents a pictorial summary of our analysis. The figure takes into account the cluster evolution with the number of He atoms and adds to the plot of the evaporation energy ⌬E as a function of N with additional information coming from our three-dimensional ͑3D͒ representation of solvent and impurity densities, discussed in previous work. 44 Up to N = 10 the dopant density ͑red color online͒ is progressively embedded in the solvent droplet: For the smaller sizes H 2 still resides on the surface ͑at our cutoff value for the He density of about 70%͒ but for N = 10 the molecule is fully immersed ͑in this case the isosurface is rendered transparent for reasons of clarity͒. For N = 20 and 22 the molecule is again surrounded by the solvent layer and gets clearly solvated into it. It is useful to point out here that because of the lightness and quasispherical shape of the H 2 dopant, our 3D method cannot establish the preferential orientation of H 2 on the cluster surface or inside it, i.e., whether the dopant is placed along the surface or perpendicularly to it. However, given the strong delocalization of H 2 and the fairly large size of its cavity within the droplet, that information would be of little significance and would not modify the general findings of the present study.

IV. PRESENT CONCLUSIONS
A Monte Carlo study on the H 2 ͑He͒ N system has been reported in this work, where we carried VMC and DMC calculations in order to extract nanoscopic information about energetics and geometrical features of these very weakly bound clusters. Several features were unearthed by our present calculations and by the analysis of our results.
• Exchange and binding energies indicate that the impurity is bound to the solvent, giving rise to a more stable aggregate than the pure one, but only by a very small amount of energy ͑less than 3 -4 cm −1 ͒.
• The evaporation energy shows a slight oscillatory behavior amplified by the smallness of its values. ⌬E is however below the bulk value of about 5 cm −1 , also for H 2 ͑He͒ 100 and no shell structuring is found to be present across the size range which we have examined in our study.
• Radial distributions indicate a large delocalization of both dopant and adatoms, in keeping with the highly quantum nature of both partners.
• The dopant is progressively embedded within the cluster as N increases, although it does not initially reside at the exact center of the droplet. The final, full solvation is here at odd with previous VMC results; 22 it may be due to our different choice of interaction potentials and to our increased accuracy of calculation when using the DMC procedure after the initial VMC results.
In conclusion, the study of mixing two highly quantum bosonic particles such as 4 He and p-H 2 , indicates that the DMC treatment is an essential tool for reliably extracting their nanoscopic behavior and that the method provides useful indications on the competitive adaptability of the two species into stable aggregates.