Coronene molecules in helium clusters : Quantum and classical studies of energies and configurations

Coronene molecules in helium clusters: Quantum and classical studies of energies and configurations Rocío Rodríguez-Cantano,1 Ricardo Pérez de Tudela,1,a) Massimiliano Bartolomei,1 Marta I. Hernández,1 José Campos-Martínez,1 Tomás González-Lezana,1,b) Pablo Villarreal,1 Javier Hernández-Rojas,2 and José Bretón2 1Instituto de Física Fundamental, IFF-CSIC, Serrano 123, 28006 Madrid, Spain 2Departamento de Física and IUdEA, Universidad de La Laguna, 38205 Tenerife, Spain


I. INTRODUCTION
Polycyclic aromatic hydrocarbons (PAH) are components of the interstellar medium, formed by means of photoevaporation processes in carbon-richened clouds.Their presence in the stellar dust, 1 the interaction with cosmic rays, 2 formation and destruction mechanisms 3 and the role played as initial precursors of organic globules in meteorites 4 have been subject of numerous investigations.
Besides the importance as finite size precursors or prototypes for graphene and related carbon layered materials, the adsorption of single atoms and molecules on the surface of these PAH molecules has been studied in detail.In particular, most of the interest lies on the understanding of these processes making special emphasis on the energy barriers probed by adsorbed atoms and the possible role played by quantum mechanical (QM) effects. 5Experimental investigations of helium clusters on aromatic molecules [6][7][8][9][10][11][12] have been accompanied by theoretical work, 5,[13][14][15][16][17] usually consisting in the application of Monte Carlo (MC) techniques, to characterise their structure and energetics.
The specific localization of the He atoms on the PAH molecules and excitations on the adsorbed layers can explain anomalous features on the observed spectra. 9,10,16,17ine splittings in hole-burning measurements of tetracene within droplets of approximately a hundred He atoms were interpreted in terms of occupation of two nearly equivalent sites by localized atoms or due to tunneling of one or two He atoms through a double well potential on the tetracene a) Present address: Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, 44780 Bochum, Germany.b) Electronic mail: t.gonzalez.lezana@csic.essurface. 10Spectral shifts observed for He 2 -anthracene were found to differ noticeably depending on whether the He atoms were localized in the same side of the PAH molecule or were placed almost equally in both sides. 9For small benzenedoped helium clusters, path integral Monte Carlo (PIMC) calculations determined that atoms occupying the potential wells are excluded from the exchange paths associated with superfluid fractions on the surface of the PAH molecule. 14ndications of superfluid behaviour were no expected in films formed with less than three thick layers of He atoms on a graphite surface. 18In this sense, properties of monolayer helium on these molecular substrates have been studied with some detail. 19Thus, for example, studies on films of 4 He on graphite revealed that in the first layer one-third of the available substrate adsorption sites are occupied and MC calculations with no explicit inclusion of the particle permutation are capable to reproduce the most significant features of the clusters. 20Despite melting and phase changes depend in general on the density and the temperature, a low-density monolayer of helium on graphite has been found to consist of solid clusters with no traces of superfluidity. 21,22The different environments offered to the solvating He atoms by aromatic molecules leads to the coexistence of a variety of situations involving different degrees of localization.
Recent investigations on PAH cations 23 revealed that the graphitic planes establish relatively strong bond with the helium atoms which thus form a frozen-like first solvation layer which can be then coated in successive layers exhibiting fluid-like behaviour.Around the peripheral regions, usually explored once the preferred minima location on the molecular substrate have been already occupied, the solvating He atoms are less firmly bound and therefore their positions are more delocalized presenting a slushy behaviour.The coexistence between confined and delocalized phases in He-PAH clusters, which may have its origin in the highly anisotropic character of the interaction between the He atoms and the aromatic molecule, is responsible for unusual features on the corresponding spectra such as zero-phonon lines (ZPL) and zero-phonon wings.
The finite extent of PAH substrates introduces the additional interest of size effects whose magnitude is connected with the actual dimensions of the corresponding cluster.In this sense, coronene, C 24 H 12 , is a medium-sized example, bigger than some other species in the same family such as pyrene (C 16 H 10 ) or benzene (C 6 H 6 ), of the same order as tetracene (C 18 H 12 ) or pentacene (C 24 H 14 ) but smaller, for instance, than circumcoronene (C 54 H 18 ).Given its structure with a central carbon ring (ZPL) bound only to other C atoms, coronene is considered as the simplest reduced representation of a graphene layer.Crucial to explain the infrared emission spectrum of interstellar dust, 24 coronene has been investigated as a molecular substrate for the adsorption of H, 5 a water molecule 25,26 or Ar. 26 Helium aggregates on a cationic version of coronene, C 24 H + 12 , have also been recently studied by means integral molecular dynamics (PIMD) calculations. 23Classical binding energies, obtained by means of the basin-hopping (BH) approach and quantum-corrected values via estimates of the zero-point energy (ZPE) from harmonic contributions for He N -coronene + clusters with N up to 120, were compared.Various delocalization factors were examined through the PIMD simulations in order to characterize the different layers coating the coronene molecule as successive He atoms were added.
A crucial aspect to take into account for the theoretical investigations of doped PAH systems is the description of the interaction between the adsorbed atoms and the molecular substrate.Thus, anisotropic and spherical Lennard-Jones potentials for the He-C interactions in the PIMC study of 4 He on a single graphene sheet yield different predictions regarding the existence of mobile vacancies on the surface inducing finite superfluid fractions. 27In addition, the importance of the surface corrugation has been analyzed in some detail for the case of helium monolayer graphene [28][29][30] and graphite. 31A global potential for the physisorption of He on the surface of coronene has been reported by some of us 32 using improved Lennard-Jones (ILJ) functions.The comparison with ab initio calculations by means of density functional theory-symmetry adapted perturbation theory (DFT-SAPT) and "coupled" second order Møller-Plesset perturbation theory (MP2C) reveals a very good agreement (within 0.1 meV) especially at the top and bridge positions in which the He atom locates over a C atom or over a C-C mid-bond, respectively.Slight deviations for the potential energy well depth corresponding to the hollow position with He over the center of the inner carbon ring were observed, with the ab initio calculations predicting minima about 1 meV deeper than the analytical potential.In this work, we have employed the He-coronene analytical potential of Ref. 32.
The He N -PAH clusters are in general very special systems in which the internal energy is ruled by a very delicate balance between the He-coronene and He-He interactions.As a result, regions of phase coexistence in clusters formed with layers of the He atoms on graphene 22 have been observed.In particular, for coronene, the equilibrium distance for V He−He , 2.97 Å, is close to the existing separation between hollows of consecutive hexagons (∼2.57Å).Inspection of local and global minima which differ in some few meV thus becomes a complicate task which requires the choice of a robust minimization method.In this work, two classical approaches, an evolutive algorithm (EA) and the BH method are employed to explore the equilibrium configurations on the analytical potential of Ref. 32.In addition, both the energetics and structures of clusters formed by the PAH and up to 69 He atoms at T = 2 K, are investigated by means of a PIMC method and ZPE corrected BH estimates.The role played by thermal effects is further investigated by comparing the PIMC results with diffusion MC (DMC) calculations for the smallest clusters.
The paper is sectioned as follows: First, we present the theoretical approaches employed in the calculations (Section II); details of the potential energy surface (PES) and the ab initio calculations performed in order to test specific minima are explained in Section III; results are shown in Section IV and discussed afterwards in Section V and finally concluding remarks are listed in Section VI.

A. Path integral Monte Carlo method
4][35][36][37] Briefly, the density matrix at a temperature T is replaced by the product of M density matrices at higher temperatures MT, where β = 1/k B T, with k B being the Boltzmann constant and η = β/M.R α is the vector which collects the 3N positions of the N He atoms: R α ≡ {r α 1 , . . ., r α N }, being r α i the position vector of the ith He atom at the time slice or imaginary time α.The center of mass (CM) of the nonrotating coronene molecule, with a much larger mass than He, is supposed to be fixed to the origin of coordinates.Under such assumptions the total Hamiltonian Ĥ is written as follows: where is the total potential, m He is the helium atom mass, is the Planck constant, and r i j = | r i − r j | are helium-helium distances.
The internal energy can be obtained by means of the estimators proposed in Refs.38 and 39 as where r C i = M −1  M α=1 r α i defines the centroid of the M beads of the ith-particle.The first term in Eq. ( 4) accounts for the classical kinetic energy, the second one is a quantum correction where F α i is the force experienced by the i-particle on the α slice and the third term describes the interaction between each pair of particles on that α slice.Finally, the integration is carried out via a Metropolis MC algorithm, as an average over a number of paths {R 1 , R 2 , . . ., R M , R M+1 } sampled according to a probability density proportional to the factorised product of M density matrices of Eq. ( 1).
The PIMC calculation was performed at T = 2 K, a sufficiently large temperature to neglect exchange permutation symmetry effects in Eq. ( 1) due to the bosonic character of the 4 He. 33,40,41Indeed previous studies on 4 He clusters on PAH also suggest that particle permutation is explicitly included only at much lower temperatures 14,15,42 and it is not required for low coverage densities of helium. 14,18,20,22,23,43,44he number of M beads vary between 50 and 500, including the M = 1 case, referred as the classical MC (CMC).The MC simulation involves between 10 6 -10 8 steps and, in order to guarantee the ergodicity of the sampling algorithm, we employ a combined procedure consisting on: (i) staging involving 8 beads in each movement and (ii) a classical displacement of all beads describing a particular atom as a whole every 10 5 steps.][47] Information about the geometry and the structure of the different clusters formed combining He atoms with the coronene molecule has also been obtained by means of probability densities as a function of the Cartesian coordinates describing the coronene plane D(x, y); as functions of y and the distance between the He atoms and such a molecular plane z, D( y, z), or simply with respect to the distance to such a plane, D(z).These probability densities are calculated by means of a histogramming procedure on the coordinates of interest.
The analysis of the spreading of the wave functions describing the systems as N changes can be performed via the fluctuations of the position of each atom i with respect to the corresponding centroid: and the rigidity of the system is measured with the Lindemann index, defined as that is, a root mean square of the bond lengths between pairs of He atoms considering their respective centroids. 23

B. DMC method
The DMC [48][49][50] method is a quantum algorithm that leads to the ground state wave function and energy.4][55][56][57][58] In this method, the time-dependent Schrödinger equation is transformed by means of the τ = it substitution, yielding a diffusion equation that for a system of N particles m k (k = 1, . . ., N ) including both the N He atoms and the coronene, reduces to where V is the PES of Eq. ( 3) and both the kinetic and potential terms are expressed in Cartesian coordinates.With the above transformation, the propagation in imaginary time changes the usual oscillatory behavior of the Schrödinger equation (with terms e −it E n / ) to a real exponential solution In the long-time limit, for large τ, if the energy values are taken as positive all terms will vanish and the ground state will be the last decaying term.Thus, the longest transient state will be proportional to the ground state wave function.Since we consider the intramolecular distances within coronene fixed, we apply the rigid body treatment developed by Buch 59 and assisted with the descendant weighting method 49 for the determination of the ground state energy.Differences with respect to PIMC results due to the rotation of the coronene are expected to be negligible since the corresponding rotational contribution is about ∼10 −3 meV.
The calculations have been carried out using QCLUS-TER, 60 the DMC code for rigid monomers developed by Sandler and Buch.Depending of the number of He atoms included, we have performed about 2-5 simulations in order to decrease the statistical error.Each simulation involved the propagation of 4000-8000 replicas over 1900-4200 time steps to generate a reliable ground state distribution, and with typical time steps between ∆τ = 40-80 a.u.

C. Evolutive algorithm
One of the methods employed to search for the equilibrium configurations of the He N -coronene clusters is the EA. 61The algorithm starts with the generation of initial populations of M individuals or clusters consisting on N He atoms aggregated to the C 24 H 12 structure.For each individual i the pair of vectors ( xi , ηi ) represents the 3N Cartesian coordinates of the atoms and the "strategy parameters" or standard deviations for Gaussian mutations, respectively.Initial values of η i = 1 and random choices for the positions within a specific range (0, ∆) are considered and each parent evolves creating a single offspring according to the following law: where j = 1, . . ., 3N; τ and τ ′ are adjustable parameters depending on the value of N; N (0, 1) is a random number from a Gaussian distribution of mean µ = 0 and standard deviation σ = 1, and N j (0, 1) means that the randomly generated number is for each component j.
Pairwise comparisons of the energy of each individual with random choices as opponents over the union of 2M elements composed of parents (x i , η i ) and offsprings (x ′ i , η ′ i ) are performed.Points are awarded in case of the energy turns to be lower than the energy from the corresponding opponent.Finally, those M individuals out of the union of parents and offsprings with a larger amount of winning points are selected as survivals to the next generation and thus becoming new parents.The procedure is repeated until the difference between the potential energies of consecutive generations is lower than certain tolerance value.

D. Basin-hopping method
The BH 62 or MC minimization 63 is another stochastic algorithm used to locate the global minimum on a PES.This method transforms the PES into a collection of basins and explore them by hopping between local minima.The BH method can be considered as an iterative algorithm composed by a finite number of steps.Each BH step consists of three parts: (1) random perturbation of the coordinates, (2) local minimization, and (3) accept or reject the new minimum according to a Metropolis criterion.This method has been used successfully for both neutral 62,64,65 and charged atomic and molecular clusters, [66][67][68][69][70][71] along with many other applications. 72uitable parameters for the present BH simulations were determined for all cluster sizes based on preliminary tests on He 14 -coronene.These benchmarks consisted of 10 4 minimization steps and were initiated from independent random geometries varying the optimization temperature.
The results below were obtained at a constant optimization temperature of T = 6 K.A total of 8 runs of 3 × 10 4 BH steps each were performed for all clusters sizes.The same global minima were obtained in all trajectories.This ensures a reasonably high degree of confidence.
In order to include nuclear quantum effects, we have calculated the ZPE in the harmonic approximation to each of the classical potential energy minima.We build a database of local minima close to the global minimum for each cluster size.It has been obtained by means of MC calculations.Hence, the "quantum" energy minimum is defined by where V BH (R e ) is the local energy minimum calculated by the BH technique and V ZPE (R e ) is the ZPE in the harmonic approximation calculated as where ω i is one of the angular frequencies of the cluster obtained by means of the diagonalization of the Hessian matrix.

III. INTERACTION POTENTIAL A. Potential energy surface
As shown in Eq. ( 3), the PES is expressed as the sum of the corresponding pairwise He-He and He-coronene interactions.For the V He−He (r i j ) interaction, we have used the potential by Aziz and Slaman 73 with a well depth of ∼0.94 meV and a equilibrium distance of ∼2.97 Å.The V He−cor (r i ) interaction has been described by means of the PES recently obtained 32 by some of the authors and based on the atom-bond pairwise additive scheme. 74,75In this approach, the van der Waals interaction between He and coronene is obtained as a sum of atom-molecular bond pairwise contributions, where the chemical bonds C-C and C-H of coronene are considered as the interaction centers on the molecular frame.The He-coronene interaction is thus given as a sum of atom-bond pair contributions (l), each one represented by a potential term V l , which depends on the distance R l between the atom and the bond reference point, as well as on θ l , the Jacobi angle describing the orientation of the atom relative to the bond axis, Each bond is assumed to be an independent diatomic sub-unit with electronic charge distribution of cylindrical symmetry.
Given that the dispersion center and the bond center for C-C are the same and almost coincident for C-H, 74,75 for the present application, the reference point for each bond is set to coincide with the geometric bond center.The pair potentials V l of Eq. ( 13) are parameterised by means of the atom-bond potentials 75 where x l is a reduced distance with ε l (θ l ) and R ml (θ l ) being the well depth and the location of the atom-bond interaction, respectively.The dependence of these parameters with the θ l angle is given by: 75 so that ε ⊥ l , ε ∥ l , R ⊥ ml and R ∥ ml are, respectively, well depths and equilibrium distances for the perpendicular and parallel approaches of the atom to the bond.With this definition, f (x l ) is taken to be the same for all relative orientations.Finally, in Eq. ( 14), the n exponent is given by the following where κ is a parameter defining the shape of the potential which depends on the nature and hardness of interacting particles.
Values for all the parameters involved in the above discussed expressions employed are shown in Table I.These parameters, derived from the polarizability of the interacting partners, were fine-tuned in Ref. 32 through comparison with high level electronic structure methods 77,78 using large basis sets.
The interparticle distances between the C and H atoms forming the coronene molecule are described as follows: the experimental C-C bond lengths and C-C-C angles appropriate for graphite were employed (1.420 Å and 120 • , respectively, 79  Despite a test of the accuracy of the present potential by comparison with experimental results, as in the case of graphene or graphite, for instance, was not possible in Ref. 32, several indications allow us to rely on its quality.Thus, the description of the equilibrium distance (z e ) and the binding energy (D e ) provided by the corresponding atom-bond curves for the top and bridge geometries is remarkably good and only slight overestimations of about 0.1 Å for z e and about 1 meV for D e , respectively, are found at the hollow configuration.In addition to this, the extrapolated interaction correctly describes the low lying bound states for He-graphite as revealed by the comparison with benchmark experimental findings (see Table 4 of Ref. 32).A similar good agreement was also found for the low lying bound states for the He-graphene system as reported in a recent investigation. 80

B. He N -coronene ab initio calculations
For some specific He N -coronene clusters ab initio calculations were also performed in order to obtain complexion energies of the corresponding minimum configurations predicted by the optimization methods.To this scope the supermolecular MP2C theory 77 was employed as implemented in the Molpro2012.1 package. 81The accuracy of MP2C has been recently assessed by extensive calculations of Rg-fullerene and Rg-coronene interaction energies 32,82 and it was concluded that the MP2C results are in good agreement with reliable DFT-SAPT 78 calculations, with the advantage of a lower computational cost.
For all calculations the aug-cc-pV5Z(He)/aug-cc-pVTZ(coronene) basis sets 83 were used, following the same procedure reported in Ref. 32.The obtained complexion energies were corrected for the basis set superposition error by the counterpoise method of Boys and Bernardi. 84

A. Classical energies and configurations
The potential minima and equilibrium configurations of the He N -coronene have been investigated by means of the EA and BH approaches of Sections II C and II D, respectively.The EA approach, slightly less demanding in computational terms, has been implemented in our PIMC code as a standard procedure to localize minimum configurations to be used as starting points for the MC simulation.However, the BH method provides more precise energy values E cl for the equilibrium geometries for the different clusters here investigated.
The configurations associated with these minima are shown in Figure 2. The equilibrium geometries suggest that He atoms locate on the coronene following certain general tendencies.For instance, if available, the hollow position of the central C-hexagon is one of the preferred sites.Thus, the He 2 -coronene cluster is formed with He atoms at both sides of the molecule plane according to a (1, 1)  II.The comparison between the global and the local minima has certainly revealed small differences (for N = 3 there is only 0.07 meV between the (2, 1) and (3, 0) configurations), an indication of the complex energy landscape of the system.
The repetition of specific locations as more He atoms are added has also an effect on the energies of the corresponding clusters.Thus the classical energy found for the (3, 1) configuration of N = 4, −57.5 meV is fairly well reproduced as the combination of the energies obtained for the (3, 0), −42.9 meV, and (1, 0), −14.6 meV, structures for the N = 3 and N = 1 cases, respectively.Analogously, the energy of the (3, 3) structure for He 6 -C 24 H 12 , −85.9 meV (not shown here) is roughly twice the value found for (3, 0), thus indicating the equal contribution coming from the three He atoms located in each side of the coronene plane.
In an attempt to test the role of the parametrization of the PES employed, ab initio calculations were performed as described in Section III B for the minima structures displayed in Table II.E ab energies for these small clusters are more negative than the E cl , a result which was somehow expected for those configurations involving He atoms in the central hollow due to the existing differences in these positions 32 between the ab initio He-coronene potential and its ILJ analytical representation.This effect is observed for the two structures of the N = 3 case, with less pronounced differences between the analytical and MP2C potential energies seen for (3, 0), in which three He atoms occupy almost bridge positions.However, not much can be globally concluded from the comparison with the ab initio results given that an accurate description of the He-He interaction still constitutes a challenge for the first principles approach employed here.The effect becomes more noticeable for configurations (N, 0) with an increasing number of He atoms, N, at the same side TABLE II.Energies in meV obtained by means of the classical minimization procedures (in column three), the ab initio calculations (in column four) and with the CMC (fifth column), DMC (sixth column), PIMC (seventh column) and BH+ZPE (eight column) approaches for He N -C 24 H 12 clusters with N = 1-4, for global and local minima.In second column, the number of He atoms above, n a , and below, n b , the plane of the coronene.of the coronene where the number of He-He pair interactions is large.Another source of discrepancy may lie in the possible role of three body effects (required to describe the interaction among the coronene and a given pair of He atoms) which are not taken into account in the strictly pairwise potential employed for the minimization algorithms.
Calculations by means of the classical approaches have been carried out up to N = 69 and results are shown in Figure 3.Total energies become more negative as N increases, with values of about approximately −700 meV for the largest He N -coronene molecules.No significant differences were observed between EA and BH results.
A much finer detail about the effect of the addition of He atoms on the energy is gained from the analysis of the energy per He atom, E He = E/N, and the evaporation energy, E evap = E N − E N −1 , both shown in Figure 4.It is worth noticing the scale for the energies in the ordinate axis (meV), which means that variations between results for the smallest cluster, N = 2, and the largest, N = 69, are in fact within the range of ∼5 meV and ∼10 meV for E He and E evap , respectively.Even within such reduced intervals for the energies, Figure 4 shows remarkable features.E He values becomes less negative following a trend which changes smoothly at N = 14, N = 26 and N = 38.Interestingly, for these clusters, E evap exhibits steps which lead the energy to plateau regions of less negative values.For a larger number of He atoms the E evap values remain below approximately −8 meV and seem to behave slightly more erratic.
Classical configurations for the equilibrium geometries of He 14 -C 24 H 12 and He 38 -C 24 H 12 are shown in Figure 5.The case N = 14 corresponds to a (7, 7) configuration in which all hollow positions for the carbon hexagons are occupied.Analogously, the cluster with N = 38 He atoms displays a characteristic configuration (19, 19) with an equal division of the total number of He atoms in both sides of the coronene.Figure 5  The apparently erratic behaviour of the E evap values as a function of N beyond N = 38 may be then due to the fact that additional He atoms have not the option to occupy specific locations at hollow, bridge, or top positions directly over the coronene structure, now screened by the presence of other solvating atoms.They start on the contrary to explore along the borders of the molecular plane thus initiating very preliminary stages for a complete solvation.An example is shown in Figure 6 for He 68 -coronene.The calculations with the CMC approach at T = 2 K provides a first insight about the thermal effects on the energy of the clusters when compared with the energies obtained with the classical minimization methods.The corresponding energies E CMC , also included in Table II, revealed differences with respect to the E cl values which seem to increase slightly with N but which do not exceed ∼2.5 meV for N = 4.

B. Quantum mechanical calculations
Much pronounced differences with the classical results are found for the DMC (T = 0 K) and PIMC energies at T = 2 K shown in Table I for the case of small clusters, thus suggesting the importance of the QM effects for these systems.The DMC method, reported as a preferred approach to calculate accurate binding energies for T = 0 K in Ref. 23, has also been used in the present investigation and results for both local and global minima are presented in Table I.A proper convergence was not found for the (4, 0) case, due in principle to the complications this technique has when various local minima are available.E DMC and E PIMC values are in reasonably good agreement being the only discrepancies due to thermal effects.The corresponding probability densities for the z coordinate showing the distance of the He atoms with respect to the plane of the coronene are also in extremely good accord.The example of such a comparison for the case of N = 2 is shown in Figure 7. ZPE estimates have been added to the BH results to compare with the present QM energies as shown in Figure 3 and in Table I.The procedure is applied only up to clusters with N = 30 given the exponential increase of the amount of data of local minima required for the analysis.The agreement with the PIMC energies is good although some deviations are observed as the number of He atoms increases with a tendency to have less negative energies for the BH+ZPE results, which become consistently upper bounds for all QM results.This could be due to anharmonicity effects introduced by the potential, already taken into account within the PIMC method but not properly well described by the purely harmonic ZPE-corrected to the BH energies.
The three QM calculations lead to some inversions with respect to the sequence of minima found with the minimization algorithms.Thus, for example, E PIMC for (2, 1) is almost 3 meV deeper than for (3, 0), the classical global minimum.Analogously for N = 4, (2, 2) becomes the most stable configuration instead of (3, 1), as predicted by the classical approaches.Similar situations are observed for clusters such as He 5 -coronene and He 10 -coronene, which seem to suggest a tendency for structures with an almost equal number of He atoms at both sides of the coronene substrate as the preferred configurations, as opposed to the classical predictions.
Values obtained with the PIMC calculations for E evap and E He are included in Figure 4 in comparison with the classical predictions.The range described with more detail corresponds to the smallest clusters N ≤ 15.PIMC E evap display a more oscillating trend as a function of N and there is no indication of the clear minima obtained with the classical approaches for He 14 -coronene.In fact, the inspection of the corresponding density probabilities D(x, y) and D( y, z) shown in Figure 8 suggests a different location of the He atoms.In particular, besides the occupation of the central hollow, the other 12 helium atoms seem to be located near those C atoms from the outer hexagons with no bond with the external H atoms in a (7, 7) configuration.
The specific size of the cluster with indications of He atoms exploring the z ∼ 0 molecular plane depends on the method employed for the calculations.Thus, the CMC approach requires N = 39 He atoms whereas the PIMC method finds that the first cluster with He on the coronene plane is He 24 -C 24 H 12 .In turn, the BH+ZPE simulations yield N = 29, a size in between, a behaviour which was also reported in investigations on doped fullerenes with a similar algorithm. 65The inspection of the corresponding distribution, shown in Figure 9, reveals how the density of the He atoms starts to approach the plane of the coronene for the case of N = 24.For this cluster all hollow positions and open dimples formed by the outer H-free C atoms between two C-H bridges (those with deeper wells in the external structure of the molecule as shown in Fig. 1), are occupied.
QM and classical treatments lead to different predictions regarding the minimum size of the He N -coronene clusters at which the He atoms start to fill up the second layer over the molecular substrate.Thus, according to the PIMC method, the first signatures of solvating atoms occupying distances above the first layer are observed for N = 45.The top panel of Figure 10 shows the D(z) distribution obtained with the PIMC method for He 45 -coronene (amplified 5 times) with some small features around | z | ∼ 6.5 Å which indicate that the last He atom is located about ∼3 Å above the first layer on the coronene.Nothing at such distances is suggested in the CMC counterpart (also shown in Figure 10), which predicts in turn that He 69 -coronene is the smallest cluster for which a second layer of He atoms over the coronene molecular plane is formed.N = 68 displays, according to a classical description, the structure of a full layer solvating the coronene (see Figure 6).As shown at the bottom of Figure 10, the classical D(z) density function displays: (i) peaks at both sides of the coronene at ∼1.5 and ∼3 Å from the plane; (ii) a maximum at z ∼ 0, indicating the presence of solvating He atoms at the edges of such a plane; and (iii) a peak at z ∼ 5.5 Å, which suggests the existence of a second layer of He atoms.The PIMC D(z) distribution for N = 69 has well developed broad peaks associated to this feature at slightly larger distances.In fact Figure 10 also reveals the differences regarding the localization of the He atoms depending on the framework employed in the analysis of the cluster.Peaks in the CMC distributions are narrow and clearly resolved with null probability areas in between, possibly indicating the different distances with respect to the coronene plane depending on the specific position (top, bridge, hollow, and border) occupied by the He atoms.The QM distributions obtained with the PIMC method, on the contrary, suggest a much larger delocalization of the He atoms around the surface of the coronene.QM results also provide evidences of localization for some specific sizes of the coronene-doped helium cluster under study.The inspection of the probability density D(x, y) for He 38 -coronene, shown in Figure 11, reveals some differences with respect to the corresponding functions for some other smaller clusters such as He 14 -coronene or He 24 -coronene (in Figures 8 and 9, respectively).Peaks for the positions of the He atoms on the hollow sites on the coronene surface look narrower for N = 38, thus suggesting a more restricted localization in comparison with the N = 14 and 24 cases.This interaction-induced localization effect has been previously observed for protons in hydrogen bonded systems. 85n fact, the analysis of the σ and δ coefficients defined in Eqs. ( 5) and ( 6), respectively, for the mentioned clusters reveals significant differences.As shown in Figure 12, the δ Lindemann coefficient for all He atoms in the N = 14 case calculated with the PIMC method is clearly beyond ∼0.15, whereas the corresponding values for N = 38 are closer to δ = 0.1, the critical point which separate the fully confined character. 23It is then possible to attribute the localized maximum peaks observed for the He 38 -coronene system to a more rigid behaviour, whereas N = 14 would display signatures of a more delocalized nature in consistence with the He atoms avoiding the hollow sites and wider peaks.In addition, Figure 12  with the occupancy of a position above the tightly packed first layer.The corresponding value of σ on the other hand indicates a pronounced spreading for this external atom.Analogously, the CMC calculation for N = 69 yields to a similar conclusion within the classical framework for the external He atom with a sufficiently different value of the Lindemann coefficient to suggest a distinct behaviour with respect to the rest.The classical estimates of these two coefficients also reveals noticeable differences in comparison with their QM counterparts.First, given that M = 1 for these calculations, the fluctuations σ are identically 0 and second, as suggested in the D(z) distributions (see Fig. 10), He atoms are more localized and therefore the δ coefficients are smaller than those obtained by means of the PIMC simulation.

V. DISCUSSION
On a recent investigation on PAH cations coated with helium clusters, 23 classical binding energies of benzene, pyrene, coronene, and circumcoronene displayed a similar trend as a function of N as those shown here in Figure 3: The addition of successive He atoms makes the energies of the He N -PAH clusters more negative.The neutral coronene considered here yields however to somehow smaller binding energies than in the case of its cationic form.The remarkable definition of E evap , results of Ref. 23 should be interpreted as indicative of a progressive instability of the PAH doped helium clusters as more He atoms are added, although the extrapolation to much larger clusters in both studies seems to approach to the expected energies in bulk helium.In particular, the binding energy per atom in bulk is about 0.86 meV 86 and evaporative cooling in free expanding helium droplets generates binding energies corresponding to the temperature of ≈0.4 K (0.034 meV).
The BH approach has been proved as a more robust method to explore the equilibrium configurations as revealed by the evaporation energies, a more detailed quantity than the binding energy, with less oscillations with respect to N than the corresponding values obtained by the EA.The corresponding QM-corrected estimates include harmonic contributions to the ZPE obtained by diagonalizing the dynamical matrix at the minimum configuration (see Section II D and Ref. 23).The comparison between CMC estimates and purely minimization values reveals that T = 2 K does not seem a high enough temperature to affect the energies for the coronene-doped helium clusters under study.In fact, the observed trend for the E BH+ZPE energies, less negative than the corresponding E PIMC values, suggests that thermal effects are not relevant in the present study.Thus a more likely explanation for the divergence between both set of energies may rely on the anharmonicity of the potential.The harmonic ZPE contributions added to the BH values for a QM estimate would then deviate from the correct energies as N increases if the potential starts to exhibit signatures of relevant anharmonic behaviour.Despite these considerations and based on the reasonable agreement found between the ZPE+BH predictions and PIMC results for the smallest clusters one could use the eigenvectors of the Hessian matrix in the BH method to obtain the corresponding density profiles and histograms within the harmonic approximation as proposed in Ref. 87.
One key issue in the study of helium clusters regards the correct treatment of the bosonic permutation symmetry in the PIMC calculation.As discussed above in Section II A, the exchange effects can be neglected in our case since traces of superfluid behaviour for these systems are usually visible: (i) for higher layers of helium over the molecular substrate and (ii) at much lower temperatures.In an attempt to have a firm confirmation of the validity of this assumption, we have performed extra PIMC calculations for specific He N -coronene clusters by means of the worm algorithm 88,89 where the exchange symmetry is included by allowing the formation of long cyclic paths from the reconnection of path integrals.Estimates of the superfluid fraction (SF) calculated as reported in Refs.90-92 for N = 14, 28 and 69 and for the pure He 64 cluster at T = 2 K, yield negligible values: ∼0.01 for the three doped systems and 0.12 for the pure He droplet, thus confirming that the temperature under consideration is too large to appreciate features of superfluid behaviour.In addition to this, the presence of the impurity is found to hinder the transition to superfluidity observed for the pure helium cluster as the temperature decreases, since the SF for the doped systems remain negligible for T ≥ 1 K and does not exceed ∼0.3 for any of the He N -coronene clusters under analysis at T = 0.2 K.
The relevance of theoretical investigations on the structure and positions of He atoms over PAH molecules in order to explain spectral signatures has been highlighted in numerous occasions before.Thus, for example, important information related to specific signatures on measured spectra of microsolvated molecules in helium droplets can be obtained by means of systematic studies of the van der Waals modes of the solvating layers depending of the size of the helium cluster. 93In the experimental study on tetracene in helium droplets of Ref. 10, the ZPL splitting was explained in terms of different configurations of localized He atoms on the molecular substrate.However, recent theoretical works by Whaley et al. 94,95 have shown that the anomalous splitting in tetracene is not due to different static metastable configurations of helium close to the molecule, but to a combination of high localization and the quasi-one-dimensional nature of highly coherent trios of helium atoms adsorbed along the long axis of tetracene.In turn, the single ZPL feature found in the electronic spectroscopy investigation of coronene inside helium droplets by Birer et al. 12 was explained assuming a non-superfluid single solvation layer as the most likely structure for the system.Our calculations are consistent with the experimental findings of Ref. 12: On one hand, they show that this first layer is composed by localized helium atoms with non-superfluid behaviour and, on the other hand, no indications of different stable structures forming this layer over the coronene on its ground electronic state that could suggest a spectral ZPL splitting have been found.
The evaporation energies as a function of N according to the classical methods exhibit a stepped structure with abrupt changes at particular sizes such as 14, 26, and 38, where plateau regions with almost constant values are defined.Such effect could suggest the existence of certain "magic" numbers for these particular N, although in principle such denomination seems to be more adequate, not for mere potential energy minima, but for proper bound or quasi-bound QM states supported for the clusters.In this sense, both PIMC and ZPE-corrected BH evaporation energies exhibit a somehow similar structure with successive steps, although the specific values of N are not strictly the same ones.In addition nothing relevant in the distributions associated to these specific cases provide concluding explanations for this behaviour.In terms of the BH+ZPE analysis, we find an increase in the ZPE contribution due to larger normal vibrational frequencies for those cases in which the E evap suffer abrupt variations.
QM effects proved to be relevant in order to find the most stable configurations as manifested by the already discussed inversions between global and local minima.DMC, PIMC, and BH+ZPE calculations evidence that there is more than a mere exploration in search for strict potential minima when the He atoms locate on the coronene structure.But besides the actual values of the energies, QM effects are also noticeable in the probability densities of the He atoms over the coronene molecule.In particular, a remarkable situation is found for He 14 -coronene with very different classical and QM probability densities.The equilibrium obtained by means of the minimization algorithms shown in Figure 5 corresponds, in principle, to a stable configuration in which the number of He atoms suffices to fill all carbon hexagons at both sides of the coronene plane.As discussed before, this feature could be the reason for the noticeable change in the tendency shown by the evaporation energies as a function of N. The He atoms locate close to the central hollow positions in a balance between the minima for the He-coronene and He-He interactions, a situation which changes according to the PIMC calculation (see Fig. 8).Beyond the strict potential minima analysis performed by the classical approaches, once QM (and possibly thermal effects) start to play a role, the He atoms seem to feel certain repulsion outwards.Instead of following a radial direction, they shift slightly ∼30 • towards some local minima supported by the coronene potential at the outer half hexagons formed with 3 C atoms and 2 H atoms (see Figure 1).

VI. CONCLUSIONS
He N -coronene clusters have been investigated by means of classical and QM approaches.Energies and equilibrium configurations are first calculated with EA and BH methods specially designed to explore the potential energy minima of the He-coronene system.The comparison of the CMC approach at T = 2 K with the minimization algorithms shows no significant differences thus suggesting a negligible thermal effect on the system.QM results suggest inversions between local and global minima for certain clusters.PIMC simulations at the same temperature and zero-point-energy corrected BH calculations yield however to less negative binding energies than classical results for the entire range of clusters under study with larger differences as the number of He atoms N increases.QM effects are also observed in the geometries of the droplets and a smaller value of N is required in the PIMC calculations to find He atoms either around the molecular plane z ∼ 0 or occupying positions in a second layer over the substrate.
FIG. 1. Contour plot of the He-C 24 H 12 potential with a distance from the He atom to the coronene surface of z = 3 Å.Cuts are shown at −13, −12, −11.5, −11, −10, −9, and −8.4 meV from red to green.The structure of the coronene is displayed in black solid line.

FIG. 2 .
FIG. 2. Equilibrium configurations of the He N -C 24 H 12 clusters, (with N ≤ 10) obtained by the minimization algorithms.
FIG. 3. Total energy of the He N -C 24 H 12 clusters in meV calculated by means of the BH (blue full circles) and PIMC (black full squares) approaches.EA (open circles) predictions are also included for comparison.The QM corrected BH results are included in red empty triangles.The region for the smallest clusters is amplified in the inset for clarity.

FIG. 7 .
FIG. 7. Probability density functions D(z) normalized to unity for He 2 -coronene obtained by means of the DMC and PIMC approaches.

FIG. 8 .
FIG. 8. Probability density D(x, y) with C 24 H 12 structure in black lines (top) and D(y, z) (bottom) for the PIMC calculation at T = 2 K for He 14 -coronene.
FIG. 10.Probability densities as a function of the z-coordinate obtained with the PIMC (enlarged ×5 in red dashed line) and CMC (in black solid line) approaches for the He 45 -coronene (top) and He 69 -coronene (bottom) clusters.Arrows show the onset of probability corresponding to the first He atoms occupying positions at the second layer over the coronene.

TABLE I .
Parameters of the He-(C-C) and He-(C-H) pair potentials for the He-coronene PES from Ref. 32 employed in this work.Units for distances R are Å, energies ϵ are measured in meV and κ is a dimensionless parameter.