{'fraction_non_alphanumeric': 0.0313968584855633, 'fraction_numerical': 0.012086967652275678, 'mean_word_length': 4.502918092171463, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Localized helium excitations in 4He_N-benzene clusters

We compute ground and excited state properties of small helium clusters 4He_N containing a single benzene impurity molecule. Ground-state structures and energies are obtained for N=1,2,3,14 from importance-sampled, rigid-body diffusion Monte Carlo (DMC). Excited state energies due to helium vibrational motion near the molecule surface are evaluated using the projection operator, imaginary time spectral evolution (POITSE) method. We find excitation energies of up to ~23 K above the ground state. These states all possess vibrational character of helium atoms in a highly anisotropic potential due to the aromatic molecule, and can be categorized in terms of localized and collective vibrational modes. These results appear to provide precursors for a transition from localized to collective helium excitations at molecular nanosubstrates of increasing size. We discuss the implications of these results for analysis of anomalous spectral features in recent spectroscopic studies of large aromatic molecules in helium clusters.


I. INTRODUCTION
Helium droplets provide a unique, ultra-cold nanolaboratory for investigation of a variety of physical and chemical phenomena. 1 This has been increasingly used in recent years to analyze the behavior of a wide variety of atomic and molecular species in a quantum liquid environment, using spectroscopic techniques to probe both the molecular and helium dynamics. 2 Electronic spectroscopy in particular allows one to probe microscopic details of the helium dynamics. In the large droplet regime (N ≥ 10 3 ), the laser-induced fluorescence (LIF) spectra of molecules in cold helium droplets are usually characterized by a sharp zero-phonon line (ZPL) due to the transition for the electronic origin, accompanied by a broad, phonon wing sideband. The zero-phonon line can contain fine structure due to internal molecular transitions, while the phonon wing structure reflects the coupling to collective excitations of the surrounding helium. Thus, electronic excitation spectrum of glyoxal, a small 6atom molecule (C 2 H 2 O 2 ), in 4 He droplets at T = 0.38 K, exhibits a distinct phonon wing feature that is separated by a gap from the zero-phonon transition. 3 In contrast, in 3 He droplets the corresponding spectrum shows a phonon wing feature but no gap. 4 The phonon wing structure for glyoxal in has been successfully interpreted in terms of excitation of the collective phonon-roton modes in large 4 He droplets, 3 while the lack of a gap between zero phonon and phonon wing features in 3 He droplets has been interpreted as consistent with the presence of particle-hole excitations in normal 3 He. 4 Recent spectroscopic experiments involving larger organic molecules in 4 He N clusters have revealed interesting additional features beyond these basic elements. 5 For a number of the larger organic molecules studied so far, both the ZPL and phonon wing portions of the LIF spectra exhibit sharp peaks superimposed on the underlying features. These additional peaks are not compatible with this basic picture of molecular electronic excitation coupling to either internal molecular transitions or bulk compressional modes of helium. 6,7,8,9 In addition to these experiments in large helium droplets, a new class of size-selective experiments involving small numbers (N ≤ 20) of 4 He atoms attached to large planar aromatic molecules has also recently emerged. 10,11 These small cluster studies allow one to directly observe the size evolution of excited states involving helium motion, at sizes less than a full solvation shell around the molecule. For these small clusters, which are better described as weakly bound complexes with helium than as quantum solvated molecules, a very different situation pertains. Starting with N = 1, the experimental spectra show discrete lines that generally increase in complexity with increasing N , with the number of observed lines reaching a maximum at N ∼ 4 − 5. After this, many of the discrete features observed for smaller N disappear, until only a few lines persist at N ∼ 10.
The larger organic molecules that have been studied experimentally vary in complexity from planar aromatics such as tetracene (a fused conjugated system of four six-membered carbon rings connected by common bonds) 6,7 to more complex heterostructures such as phthalocyanine 9 and indole derivatives. 8 The presence of aromatic character due to π-electron conjugation provides a common feature in these systems. Because of their geometry and π-electron character, planar aromatics such as tetracene are particularly interesting for study in helium droplets because they can be considered as nanoscale precursors to a bulk graphite surface and their local quantum solvation structure concomitantly as a nanoscale precursor of the adsorption behavior of thin helium films on graphite. A considerable body of literature has been accumulated for helium films on graphite, 12,13,14,15,16,17,18,19,20,21,22 so that analysis of the solvation structure and helium excitations as a function of increasing size of polyaromatic molecule offers the possibility of developing a microscopic under-standing of the evolution from quantum solvation of an isolated molecule, to adsorption and film formation at a bulk surface in superfluid helium.
The unusual features recorded in experimental spectra for these organic molecules consist of unidentified peaks in the phonon wings of vibronic bands 9 and anomalous splittings of the zero-phonon lines. 6,7,8 It has been speculated that both of these types of features may be due to some type of excitation of helium atoms that are localized at the molecular surface, 7,8,9 but the true origin of both of these types of features is unclear. Supporting evidence for excitations of localized helium atoms derives from the similarity of the energies for phonon wing peaks with low-frequency modes of thin films on graphite (4 − 11 K) observed via neutron scattering 23,24 as well as from theoretical predictions of spatially localized helium atoms at an aromatic ring in the benzene molecule. 25 Path integral calculations reveal this to be a true localization of helium atoms that are effectively completely removed from participation in the superfluid, thereby constituting a single "dead" atom in the surrounding solvation layer. The characteristics of excitations of such localized helium atoms around aromatic substrates are expected to be very different from the collective helium modes found in large droplets, 26,27,28 but to show increasing similarity with surface localized modes in helium films 23,24 as the size of the molecular nanosubstrate increases.
Theoretical analysis of the phenomenology of these helium excitations near a molecular nanosubstrate is rendered difficult due both to the lack of accurate interaction potentials for these larger molecules with helium, and to the need for accurate calculation of excited states in a very inhomogeneous quantum liquid environment. Only a few direct calculations to elucidate the nature of these helium excitations in the presence of a large molecular impurity have been reported so far, and these have been restricted to very small numbers of helium atoms (N = 1, 2). 29,30,31 Calculations involving doped clusters require accurate helium-impurity potential energy surfaces, of which there is a definite lack for the larger organic molecules. Thus, most calculations with large molecules resort to the use of simple atom-atom pair potential models. Given a potential surface, excited state calculations employing basis set methods are still limited to small N , where N is typically no greater than two. 29,30,31 For larger numbers of helium atoms that are however not yet in the regime of bulk films, a practical approach is provided by calculations for excited states that are based on quantum Monte Carlo methods.
In this work, we present a series of such quantum Monte Carlo studies of the ground and excited states of small helium clusters (N = 1, 2, 3, 14) containing a benzene impurity. Benzene is a simple planar aromatic molecule, and it can be viewed as a primitive sp 3hybridized unit of a bulk graphite surface. Thus, it is the logical starting point for the kind of systematic study of helium adsorption on molecular nanosubstrates mentioned above. We note that this notion has been pursued extensively for the heavier noble gases, 32 but it is only recently that experimental data for helium has become available. We employ here both the well-known groundstate diffusion Monte Carlo method and the projection operator imaginary time spectral evolution (POITSE) method for excited states. 33 The latter approach allows exact excitation energies to be obtained, provided a highquality trial function approximating the ground state is available. We calculate a range of helium excitation energies at T = 0 K and analyze these in terms of their spatial nature and extent. The previous path integral calculations for the 4 He 39 -benzene system at a temperature of T = 0.625 K have shown that the effect of the strength and strong π-anisotropy of the He-benzene interaction serves to localize a single helium atom on each side of the benzene surface. 25 That is, in the Feynman path integral representation, a single helium atom attached to the benzene surface is completely removed from permutation exchanges with the surrounding helium environment. In this work, we find a set of collective helium vibrations which have energies of up to ∼ 23 K above the ground state, and can be characterized in terms of their transformation properties under the symmetry group of the He-benzene interaction potential. We show that even for sizes as small as N = 14 helium atoms, one can distinguish a subset of higher-energy excitations localized near the global minimum of the helium-benzene interaction potential, from a subset of lower-energy collective excitations which are delocalized around the periphery of the molecule surface. The former appear to be associated with the localized helium density identified from path integral calculations in Ref. 25. Sec. II begins with a discussion of the model Hamiltonian and potential surface. Technical details of the quantum Monte Carlo methodology are presented in Sec. III. Our results for ground and excited states of 4 He N -benzene (N = 1, 2, 3, 14) are presented in Sec. V, where we analyze the nature of these molecule-induced localized states around the benzene impurity and, and in Sec. VI we discuss the implications for helium excitations on larger aromatic molecules.

II. THE MODEL HAMILTONIAN
We treat the 4 He N -benzene cluster as a collection of N indistinguishable helium atoms, and a rigid benzene molecule which is free to translate and rotate in space. Thus the benzene intramolecular degrees of freedom are held fixed, which implicitly assumes a separation between the helium motion and the benzene vibrational modes. The positions of the N helium atoms in a space-fixed (SF) frame of reference are denoted as {R 1 , R 2 , . . . , R N }, and the SF position of the benzene center-of-mass is denoted as R I . The HamiltonianĤ is (3N + 6)-dimensional: there are 3N helium translational degrees of freedom, plus an additional six dimensions due to translation and rigid-body rotation of the benzene molecule. It takes the formĤ =T (imp) +T (He) +V , whereT (imp) ,T (He) , andV are the impurity rigid-body kinetic energy, helium kinetic energy, and total potential energy, respectively. The impurity kinetic energy is most easily expressed by introducing a moving bodyfixed (BF) frame (x,ŷ,ẑ), whose origin relative to the SF frame is fixed at the benzene center-of-mass R I . The Euler angles (ϕ, ϑ, χ) specify the orientation of the BF axes, which are set parallel to the benzene principal axes. The benzene kinetic energy is thuŝ with the prefactors The first term ofT (imp) is the center-of-mass translational kinetic energy for a benzene molecule, where we use a mass m I = 78.114 amu. The remaining terms give the benzene rigid-body rotational kinetic energy, where the prefactors d a = 0.188 cm −1 and d c = 0.0938 cm −1 are the (oblate) symmetric top rotational constants. The Laplacian ∇ 2 I is taken with respect to the translations of the benzene center-of-mass. The angular second derivatives are taken with respect to rotations of the benzene about the BF axes (x,ŷ,ẑ). Similarly, the next term of the Hamiltonian is the helium kinetic energy, where m 4 is the 4 He mass, and ∇ 2 j denotes the Laplacian taken with respect to the displacement of helium atom j.
The final term of Eq. (1) is the model potential energy. We use an additive, two-body potential, The helium-helium interaction V He−He is given by the semi-empirical HFD-B potential of Aziz et al., 34 and depends only on the distance r ij between helium atoms i, j. The helium-benzene interaction V He−I is most conveniently expressed in terms of helium BF frame coordinates {r 1 , r 2 , . . . , r N }. It is an analytical fit to ab initio MP2 data of Hobza et al., 35 and consists of a sum of atom-atom terms: The indices α and β run over the carbon and hydrogen atoms of the benzene molecule, respectively. The heliumhydrogen interaction V H is a standard Lennard-Jones 6-12 form which depends on the distance r βj between hydrogen atom β and helium atom j, The helium-carbon interaction V C has a modified angledependent Lennard-Jones 8-14 form, where ǫ ′ C = ǫ C cos 2 θ, σ ′ C = σ C cos θ, and θ is the spherical polar angle with respect to a coordinate frame centered on the carbon atom α and parallel to the BF frame. The potential fit parameters are listed in Table I. This helium-benzene potential possesses two equivalent global minima of −94.97 K at z = ±3.27Å along the benzene C 6 -axis. There are also 12 equivalent secondary minima of −62.54 K located between neighboring hydrogen atoms, six above and six below the benzene plane. The global minima are connected to the secondary minima by saddle points of −50.77 K, situated approximately above and below the C-H bonds. A cut of the potential at z = 3.27Å is shown in the top panel of Fig. 1, where the (x, y) positions of the carbon and hydrogen atoms on the z = 0 plane are also marked. The lower panel of Fig. 1 shows a cut of the potential along the y = 0 plane, which reveals one of the secondary minima at x = 2.84Å, z = 2.70Å, and one of the saddle points at x = 1.97Å, z = 3.54Å.
For N = 1, the model Hamiltonian of Eq. (1) belongs to the molecular symmetry group D 6h (M ), which is the subgroup consisting of only the feasible elements of the complete nuclear permutation and inversion group (CNPI). 36 This group is isomorphic to the D 6h point group, which is valid in the limit of small amplitude helium motions. For N > 1, the symmetry group is identically D 6h (M ) due to Bose symmetry, since the only allowed rovibrational states are those obtained from the tensor product of irreducible representations of D 6h (M ), with irreducible representations of the symmetric group S N which remain unchanged under permutation of iden- Three representative sites s γ ′ , s γ ′′ , s γ ′′′ are marked, which correspond to centers of the helium-benzene localization factors centered near the global, secondary, and tertiary potential minima, respectively (Secs. IV and V A). tical 4 He nuclei, i.e. the totally symmetric representation of S N .

III. COMPUTATIONAL METHODOLOGY
In this section we present a summary of the numerical methods employed in our study of the 4 He N -benzene system. Variational estimates of ground-state energies are obtained from variational Monte Carlo (VMC), and exact ground-state properties are computed from diffusion Monte Carlo (DMC) with importance sampling. 37,38 The projection operator, imaginary time spectral evolution (POITSE) approach is used to obtain excited-state energies. 33 Both the DMC and POITSE methods can provide results that are exact, to within a systematic time step error.

A. Ground states: variational Monte Carlo (VMC)
The starting point for our study of the 4 He N -benzene system begins with the Rayleigh-Ritz variational theorem, where the trial energy E T with respect to some parameterized trial function Ψ T represents an upper bound to the exact ground-state energy E 0 . In the coordinate representation, this becomes where R is a (3N + 6)-dimensional coordinate denoting the positions and orientations of all bodies governed by the HamiltonianĤ. The quantity E L is a local energy, defined as An optimized variational upper bound is obtained by varying the parameters of Ψ T to minimize E T . For a realistic N -particle helium system, Monte Carlo methods offer a practical means to compute the multidimensional integral of Eq. (10). In variational Monte Carlo, expectation values of observables are evaluated as averages over the normalized distribution which is numerically represented as a discrete ensemble of Monte Carlo walkers {R k }. Typically, we use an ensemble size of n w = 2000 walkers. This distribution can be generated from a Metropolis walk, in which a move R → R ′ is proposed with a transition probability consisting of factors having the form The vector quantity F is commonly referred to as the "quantum force", and is chosen to be G D can be viewed as the Green's function associated with a diffusion/drift process in the presence of an external field F, over a time step of ∆τ . All VMCbased computations reported here use a time step of ∆τ = 75 − 100 Hartree −1 . For the 4 He N -benzene system the full (3N + 6)-dimensional transition probability is comprised of N such factors with diffusion constant D = D 4 corresponding to the N helium atoms, a factor with D = D I corresponding to the benzene impurity, and three analogous one-dimensional factors corresponding to individual rotations of the benzene about each of its three principal axes. A proposed move R → R ′ is subsequently accepted with probability p = min(1, q), where the acceptance ratio q is With the choice of transition and acceptance probabilities given by Eqs. (13)-(15), the Metropolis walk converges to the asymptotic distribution Ψ 2 T (R). At this point, expectation values for quantities such as the variational trial energy E T are sampled from this distribution as where the notation . . . walk denotes a statistical average over the course of the equilibrated Metropolis walk. This is performed by sampling the quantity given in . . . walk every M time steps apart, with M chosen to be longer than the autocorrelation length in order to minimize correlation biases in the error estimates; typically M = 50 for our VMC computations here.

B. Ground states: diffusion Monte Carlo (DMC)
With diffusion Monte Carlo (DMC), one can move beyond the variational level of theory to obtain exact ground-state expectation values for the energy and quantities which commute with the HamiltonianĤ. This method derives from the imaginary time (τ = it/h) Schrödinger equation, where Ψ(R, τ ) is the many-body wave function, and E r ef is an arbitrary constant shift in the energy spectrum. Importance sampling is introduced by multiplying both sides of Eq. (17) by a trial function Ψ T (R), and rewriting in terms of f (R, τ ) = Ψ T (R)Ψ(R, τ ) to obtain a set of equations having the form The stationary solution is now the normalized "mixed" distribution where φ 0 is the exact many-body ground-state wave function. In DMC, the approximate short-time Green's function G which generates this distribution is Here, G D is the diffusion/drift Green's function of Eq. (13), and G B has the form All DMC-based computations here use a time step of ∆τ = 25 Hartree −1 . A proposed move R → R ′ is accepted or rejected according to Eq. (15), to ensure that the exact ground state is sampled in the limit of perfect importance sampling, i.e. when Ψ T (R) = φ 0 (R). For a non-exact trial function however, the approximate Green's function of Eq. (20) results in a systematic time step error. Operationally, the importance-sampled DMC procedure is similar to that described previously for VMC, except now associated with each walker R k is a cumulative weight w k due to the action of G B : The efficiency of the DMC method can be significantly improved by replicating walkers with large w k , and destroying walkers with small w k . At every 10 − 50 time steps, the ensemble is checked for walkers whose weight exceeds the empirically set bounds w min and w max . A walker R k with weight w k > w max is replicated into n k = int(w k + u) walkers, where u is a uniform random number on [0, 1). These new walkers are then assigned the weight w k /n k . A walker R k with weight w k < w min is destroyed with probability 1 − w k ; otherwise its weight is set to unity. This is to ensure that the branching scheme, on average, does not artificially alter the ensemble sum of weights W tot = w k . The parameters w min and w max are chosen empirically to maintain a stable DMC walk with respect to W tot and the ensemble size; here we use w min = 0.1 − 0.3 and w max = 2.0 − 2.2. Additionally, we also vary the reference energy E r ef during the course of the walk according to . (23) Here, the imaginary-time dependence of the E r ef and W tot is explicitly written. We set E r ef initially to be the starting ensemble average for the local energy E L . The empirical update parameter η is chosen to be as small as possible, typically η/∆τ = 0.0004 − 0.004, which gives fluctuations in W tot of roughly 5%.
Once the ensemble converges to the mixed distribution Ψ T (R)φ 0 (R), mixed expectation values for a Hermitian observableÔ can be obtained as In our ground-state DMC calculations, this average is computed by samplingÔ every M = 150 time steps apart, which is larger than that used for VMC because successive DMC iterations are more strongly correlated due to the smaller DMC time steps that we use here.
In the case where φ 0 is an eigenstate ofÔ, i.e. whenÔ commutes with the HamiltonianĤ, this procedure yields exact expectation values. All ground-state energies reported in this work are derived from the mixed estimator for the local energy, and are thus exact. For quantities which do not commute withĤ, such as the local helium density operatorρ(r) = j δ(r − r j ), the mixed estimator is biased by the trial function Ψ T . But as pointed out by Liu et al., 39 the asymptotic weight w k (τ → ∞) of a walker R k is proportional to φ 0 (R k )/Ψ T (R k ). Thus, in principle this trial function bias can be eliminated by computing instead the reweighted average where M ′ is taken to be as large as possible, typically M ′ = 1500. When branching processes are incorporated in the DMC, one needs to take care to keep track of walkers that descended from a walker R k at time τ . All helium densities reported in this work (N ≤ 14) are computed by reweighting walkers R k by their descendant weights, i.e. using Eq. (27). However, the statistical noise seen in the densities computed in this manner is much greater than from densities obtained using the mixed estimator, and thus this approach becomes problematic for 4 He N -benzene at still larger sizes.

C. Excited states: projection operator imaginary time spectral evolution (POITSE)
The POITSE approach is a DMC-based method which can provide exact excited state energies, subject to a systematic time step bias. 33 It begins with the DMC evaluation of the imaginary-time correlation functioñ Here, {φ n } and {E n } represent a complete set of energy eigenstates and eigenvalues of the HamiltonianĤ, respectively, andÂ † is an operator chosen to connect Ψ T to the excited state(s) of interest φ n . Ψ T should be a good approximation to the exact ground state φ 0 . To second order in x,κ(τ ) is a superposition of exponential decays whose decay rates correspond to the energy differences E n − E 0 . 40,41 The VMC procedure of Sec. III A is used to generate an initial ensemble of n w = 2000 walkers distributed according to f (R, 0) ∝ Ψ 2 T (R). This distribution is then propagated in imaginary time using the DMC procedure outlined in Sec. III B, during which the correlation functionκ(τ ) is sampled from the DMC walk as 41κ The index k ′ denotes walkers at time τ which descended from a parent walker k at time τ = 0, and here we also explicitly emphasize the imaginary-time dependence of the DMC quantities W tot , n w , and w k . An inverse Laplace transform of Eq. (31) yields the spectral function whose peak positions give the excited state energies E n − E 0 . The additive contributions O(x 2 ) are neglected in Eq. (34), since they do not affect the positions of the peaks associated with E n − E 0 , and in practice, their spectral weights are negligible for a reasonable choice of Ψ T . 33,40 The numerical inverse Laplace transform is performed using the Maximum Entropy Method (MEM), in which the Laplace inversion is formulated as a data analysis problem in terms of Bayesian probability theory. 42 We adopt here the vector notation κ ≡ {κ(E i )∆E} to denote the discrete spectral function specified at intervals of width ∆E, andκ ≡ {κ(τ i )∆τ } to denote the correlation function sampled at intervals of ∆τ . The Bayesian approach begins with the posterior distribution function p(κ|κ, α, m, I) ∝ e αS−L .
The posterior p represents the probability of obtaining the image κ, given the Monte Carlo dataκ, the parameter α, an initial guess for the image m (also referred to as the default model), and any other relevant background information I. We use a constant value for the default model m. The quantity S is the Shannon-Jaynes information-theoretic entropy, and L = χ 2 /2 is chosen to be the usual Gaussian likelihood. In matrix notation, χ 2 takes the form whereL is the Laplace operator, andC is the covariance matrix for the dataκ. In this work,C is obtained from averaging N d = 256 decays, which is sufficient to give well-converged spectra for N ≤ 3. A search is then made for the imageκ which maximizes the posterior p, using a search algorithm due to Bryan. 43 In the limit α = 0, this reduces to a standard least-squares data-fitting procedure involving the minimization of χ 2 , which is numerically unstable when one seeks to infer a free-form solution for the image from data with non-negligible Monte Carlo noise. For α > 0, the search forκ requires that the entropy S be simultaneously maximized while χ 2 is minimized. Thus, α can be viewed as a regularization parameter, which stabilizes the least-squares fit by constraining the minimization of χ 2 .
In the ideal situation where the spectrum consists of a single, sharp peak, we take the excitation energy E n − E 0 to be the first moment of the spectrum, The corresponding estimate for the variance of the mean σ 2 then follows from the usual procedure for the propagation of uncertainties, which requires knowledge of the covariance matrix C for the imageκ. By approximating p as a sharply-peaked Gaussian in the neighborhood ofκ, the covariance of the image becomes 44 When the spectrum is composed of multiple, wellseparated peaks, the mean value of the peak and the variance in the mean is obtained in the same manner as given by Eqs. (38)- (39), except that the summation is taken only over the spectral feature of interest. Note that the approach taken here in evaluating error estimates differs somewhat from the Gaussian approach advocated in Ref. 45. In the Gaussian approach the variance in E n −E 0 is associated with the peak width, and thus does not necessarily scale as 1/N d , where N d is the number of DMC decays used to compute statistics for the covariance ma-trixC of Eq. (37). Here, we obtain E n − E 0 from the mean peak position, and then evaluate the variance of the mean according to Eqs. (38)- (39), which by construction scales as 1/N d .

IV. TRIAL FUNCTIONS AND EXCITATION OPERATORS
While in principle a fully converged DMC result should be independent of the trial function Ψ T , a good choice of Ψ T which closely approximates the exact ground state φ 0 can significantly improve computational efficiency. On the other hand, a poor choice of Ψ T can produce misleading results. In this study, we use analytical forms motivated by basic physical considerations. The trial function is a product of two-body factors, Here, χ describes helium-helium correlations, and has the McMillan form 46 The function ξ describes helium-benzene correlations, and is a product of helium-benzene factors defined in the BF frame, The first of these factors ξ 0 controls the behavior of the trial function at short and long helium-benzene separations: The terms involving the parameters c C , c H ensure that Ψ T goes to zero as a helium atom approaches a carbon or hydrogen atom, where r αj and r βj denote helium-carbon and helium-hydrogen distances, respectively. The r −6 and r −5 forms are chosen to cancel the leading singularity in the local energy 47 due to the helium-benzene pair potential of Eq. (6). With ξ = ξ 0 alone, it is known that the resulting trial function describes a ground state with "liquid"-like characteristics. If the true ground state has "solid"-like qualities, i.e. individual helium atoms are strongly localized around the molecular impurity, such a state will not solidify out of a liquid-like trial function at the variational level of theory, within reasonable VMC simulation times. One simple way to incorporate solid-like features a priori is to introduce a set of sites {s γ } fixed in the BF frame. The remaining factors ξ γ (γ > 0) serve to localize helium atoms near these sites. In this work we use the exponential of a Gaussian to effect this localization: 48 The specification of the sites {s γ } is given later in Sec. V A. A more sophisticated variational approach would involve the extension of this trial function to shadow functions, in which the {s γ } are treated as subsidiary variables, thus obviating the need for an a priori specification of {s γ }. 49 While a shadow trial function is expected to yield a significantly better description of the ground state, we do not employ this approach here because it introduces additional complexity due to the need to sample the subsidiary variables {s γ } in VMC and DMC. In the POITSE formulation, an initial ansatz for the excited state is generated by the action of an excitation operatorÂ † on Ψ T , and exact excited-state energies are then extracted from the DMC imaginary-time evolution of this initial state. Due to the inherent numerical difficulties in the MEM inversion when the target image κ(E) contains multiple, closely-spaced peaks of comparable spectral weight, we seek operatorsÂ † which connect Ψ T to only one or a few energetically well-separated states. Thus, one consideration is to chooseÂ † to transform as an irreducible representation of the D 6h (M ) group, such that the spectral weight φ n |Â † |Ψ T is only non-zero for states φ n which transform identically. For systems with a relatively high degree of symmetry, such a choice of A † would significantly reduce the number of individual decays contributing toκ(τ ).
In the BF representation, a convenient primitive for A † are provided by the regular spherical harmonics: Similar operators have been used previously to study rotationally excited states of small 4 He N and (H 2 ) N clusters (N = 7, 40). 50,51 For a pure cluster, this yields trial excited states which are simultaneous eigenstates of the square of the total helium angular momentumL 2 , and its SF z-componentL z . Addition of a benzene impurity lowers the symmetry to D 6h (M ). In this case the set of regular spherical harmonics can be symmetrized by taking linear combinations of R lm and R l,−m to obtain real operators which transform as irreducible representations of D 6h (M ): The subscripts c and s indicate proportionality to cos(mφ) and sin(mφ), respectively. The various R lm are listed in Table II in terms of BF Cartesian coordinates for up to l = 4, along with their corresponding symmetry labels under D 6h (M ). Note that for up to l = 4 under D 6h (M ), R lmc and R lms are doubly degenerate except for m = 3, and thus from this point on we drop the subscripts c and s for the degenerate functions except when the distinction is necessary. For N > 1, these one-particle operators are Bose-symmetrized over the N identical bosons to give single-excitation operators: Application ofÂ † lm on Ψ T then gives the trial excited state: = N j=1 ξ(r 1 )ξ(r 2 ) . . . R lm (r j )ξ(r j ) . . . ξ(r N ) N i<j χ(r ij ). x, y E1u R20 x 2 + y 2 − 2z 2 A1g R21c, R21s xz, yz E1g R22c, R22s Thus, the action ofÂ † lm is to impose a trial nodal structure on the helium-benzene factor ξ, generating a trial excited state corresponding to a Bose-symmetrized singleexcitation excited state. As will be discussed further in Sec. V B, these operators can be further refined to take into account the localization of the helium density near the molecule surface. In principle this approach can also be readily extended to multi-excitation operators.

A. Ground states
Ground-state energies and structures for the 4 He Nbenzene cluster (N = 1, 2, 3, 14) were computed using the methodologies described in Secs. III A and III B. Our general strategy for ground states is to first begin with a preliminary unbiased DMC calculation with Ψ T = 1, which gives an asymptotic DMC distribution proportional to the exact Bose ground state φ 0 . However, such an unbiased DMC approach for helium clusters has a number of convergence issues associated with it. In the absence of a guiding function, the DMC walk will spend too much time sampling unimportant regions of configuration space, and thus unbiased DMC can give an energy which is not converged and is higher than the true value of E 0 . 38 Nevertheless, it can be useful in situations where no other information is available, in particular to act as a guide for constructing a suitable starting trial function Ψ T . We first compute a reduced ground-state wave functionφ 0 (r) in the BF frame by binning the unbiased DMC distribution φ 0 (R) onto a three-dimensional grid. Trial function parameters for the helium-benzene factor  (45), in atomic units. The position of a site sγ is listed in terms of its (x, y, z) Cartesian coordinates in the BF frame. Only one center from a set of sites which are equivalent by symmetry is given, and the remaining can be generated by application of the elements of the D 6h (M ) group.    Table III. Unbiased DMC refers to DMC with ΨT = 1. IS-DMC is the best estimate for the exact ground-state energy, using the trial function whose corresponding VMC energy is given in the second column. The value in parenthesis represents the standard deviation in the last significant figure. ξ, Eqs. (43)- (45), are then obtained by fitting to this φ 0 (r). In some cases, these parameters are further varied to lower the variational energy, and the resulting trial function is then used as input for an importance-sampled DMC calculation. The final optimized trial function parameters for the various sizes are listed in Table III. The VMC and DMC results for the ground-state energy per particle E 0 /N are summarized in Table IV. For N = 1 − 2, we find that a simple trial function involving ξ = ξ 0 only, Eq. (44), is sufficient to give good DMC results. With the two global minima occupied at N = 2, additional helium atoms begin to occupy the twelve secondary potential minima (Fig. 1). To describe this situation for N = 3, we incorporate localization factors ξ γ , Eq. (45), centered at sites {s γ ′ } and {s γ ′′ }. The set of sites {s γ ′ } are centered near the two equivalent global minima of the helium-benzene potential, the set {s γ ′ } are centered near the twelve equivalent secondary potential minima. The N = 14 trial function is determined in the same manner as N = 3, except that we now also incorporate additional localization factors centered at sites {s γ ′′′ } which are situated near the six equivalent tertiary potential minima. This allows for a description of helium binding in these regions in the larger clusters. Three representatives of the sites s γ ′ , s γ ′′ , s γ ′′′ are marked on the potential contour plot of Fig. 1. The positions of all corresponding sites are readily generated from application of the elements of the D 6h (M ) group.
Cuts of the helium density along the y = 0 plane are shown in Fig. 2. These cuts are taken along a plane perpendicular to the benzene molecular plane, bisecting two parallel carbon-carbon bonds. All structural quantities given here are obtained from DMC using the descendant weighting procedure described in Sec. III B, and are thus sampled from the exact ground-state density distribution φ 2 0 . For N = 1, 2, the helium density along the benzene C 6 -axis has two maxima at z = ±3.7Å, which is in good agreement with the corresponding value of z = ±3.5Å inferred from spectroscopic measurements. 52 The positions of these two density maxima remain unchanged as N increases to 14. As evident in Fig. 2, the local density near the benzene impurity is highly structured, which is reflective of the strong anisotropy in the helium-benzene interaction potential. We note that neglecting the rotational terms in the benzene kinetic energy, Eq. (2), leads to a considerably more strongly peaked helium density distribution in the BF frame. 53

B. Excited states
We have calculated helium vibrational excited-state energies for the 4 He N -benzene system, using the POITSE methodology described in Sec. III C. For the N = 1 dimer, each individual excitation operatorÂ † lm gives a well-defined excitation, by which we mean a spectral function κ(E) consisting of a single, sharp peak. This indicates that the trial excited state given byÂ † lm Ψ T is a good approximation to a true eigenstate of the system, having non-negligible overlap with only a single φ n . The various spectral functions are superimposed and shown in the upper panel of Fig. 3. For N = 1, these helium vibrational modes begin at ∼ 9 K above the ground state. Below this onset, we would expect to see excitations involving rotational motion of the much heavier benzene impurity, which lie in the spectral range studied experimentally by Beck et al. 52 The excitations shown in the upper panel of Fig. 3 can be grouped into pairs split by ∼ 0.6 − 2.3 K, corresponding to states generated by the application of pairs of projectors which are symmetric and antisymmetric with respect to reflection about the benzene molecular plane (Table II). Since the N = 1 projectors used here all give well-defined excitations, we conclude that these states represent symmetric and antisymmetric tunneling pairs. We are not able to extract the lowest tunneling excitation given by the application of the R 10 = z operator. This is likely due to its en-ergy being too close to the ground state E 0 , and thus its DMC imaginary-time evolution is too slow relative to the DMC propagation times used for the excited-state calculations here. These tunneling energetics are comparable to those obtained from basis set calculations of the 2,3-dimethylnaphthalene·He complex by Bach et al. 29 These authors found tunneling excitations associated with the motion of the single complexed helium atom from one side of the planar aromatic moiety to the other side, with splittings ranging from < 0.0014 K for strongly localized states, and up to 4.6 K for highly delocalized states.
The Bose-symmetrized versions of the excitation operators used for N = 1 give a similar set of well-defined excitations for N = 2, which are shown in the lower panel of Fig. 3. Thus, we conclude that these represent singleparticle-like excitations, which is reasonable since the two helium atoms reside in the two equivalent global minima along the benzene C 6 -axis, and are well-separated by the benzene plane. In general, the N = 2 excitation energies E − E 0 tend to be slightly lower than those for N = 1, with the onset of helium vibrational excitations beginning at ∼ 7.5 K for N = 2, as compared to ∼ 9 K for N = 1. Unlike the situation for N = 1, the addition of a second particle now allows for a well-defined excitation of A 2u symmetry at 7.63(4) K, obtained from the R 10 = z operator. We note that neglecting the rotational terms in the benzene kinetic energy, Eq. (2), alters the energy spectrum significantly. 41 With the two global minima occupied at N = 2, an additional third helium atom must be delocalized over the twelve secondary minima, due to the effect of heliumhelium repulsions. Thus, the character of the excited states is expected to change dramatically as the number of helium atoms increase from N = 2 to N = 3. This is evident in the POITSE calculations, where the projectors that give well-defined excitations for N = 1, 2 ( Fig. 3) now give multiple peaks for N = 3, shown in the upper panel of Fig. 4. This indicates that the starting excited state ansatz generated byÂ † lm Ψ T now has appreciable overlap with more than one eigenstate. Since the POITSE methodology does not provide information on excited-state wave functions, interpretation of the N = 3 excitations is now less straightforward. Additional insight can be gained by making modifications to the projectors. The most noticeable new feature of the N = 3 spectrum here is the appearance of lower-energy states at < 9 K. Since the two atoms situated near the global potential minima experience a very different local environment from the third atom that is delocalized over the twelve secondary minima, the question then arises as to whether one can ascribe any of the excitations to motion localized near the global potential minima alone. We explore this here by defining a set of weighted "local" excitation operatorŝ where the product in the denominator is taken over all localization factors that are not centered near the global minima. Application of theseB † -type projectors on Ψ T gives the following local excited-state ansatz: By substituting Eq. (43) for ξ(r j ) in Eq. (54) above, it can be seen that the action ofB † lm is to place an atom in an excited single-particle state that is spatially localized near the set of sites {s γ ′ } only, i.e. near the two global potential minima. These operators are local in the sense that they act primarily on helium density near these two minima, while still maintaining spatial and Bose permutation symmetry. In contrast, by comparison with Eq. (51), we see that theÂ † operators are global opera-tors, acting on all sites.
The spectrum for N = 3 generated by the set ofB †type projectors is shown in the lower panel of Fig. 4. Compared to the spectrum derived fromÂ † -type projectors (upper panel of Fig. 4), the spectral weight of lowerenergy states are now reduced relative to the higherenergy excitations. A few low-energy features at < 9 K persist, in particular the R 10 and R 11 excitations, albeit with reduced weights. Thus we conclude that the higherenergy excitations above ∼ 9 K are associated with states that are spatially localized near the global minimum of the interaction potential, while the excitations below this range are associated with collective helium states of a more delocalized nature. For each set of states generated by a projector of given l, corresponding states of different m are clustered more closely together than in the N = 1, 2 spectra, with each cluster of different l spaced ∼ 3 K apart.
As N increases from 3 to 14, the ground-state DMC calculations reveal that the two global minima and twelve   (Table II). secondary minima are completely occupied (see Sec. V A and Fig. 2). An unambiguous determination of excited state energies at this size becomes more difficult due to the increase in Monte Carlo noise. However, certain qualitative features remain apparent. The upper panel of Fig. 5 shows the POITSE excitation spectrum for N = 14, derived from theÂ † -type global projectors. These projectors generate trial excited states which overlap with multiple eigenstates, and the resulting spectrum for a given projector has multiple peaks at both low and high energies. This problem, along with the accompanying increase in statistical noise in going to N = 14, gives a spectrum which is not fully converged in the sense that the peak widths are broader and the peak positions shift (by up to ∼ 3 K) with additional sampling. Nevertheless, we can again achieve qualitative insight by comparing with spectra derived from the localizedB † -type projectors. The lower panel of Fig. 5 shows the spectrum derived from theB † -type projectors. There, the spectral weights of low-energy excitations (< 9 K) are again significantly reduced just as the N = 3 case, and each projectorB lm gives a κ(E) consisting of a single dominant peak. While the higher-energy l = 4 excitations are   (Table II).
probably still not fully converged, the l = 1 − 3 peaks do not shift much (less than 0.5 K) with additional sampling, and appear to be converged. Thus the qualitative trends can be clearly established. Similar to the N = 3 case, the localized projectors give rise to a set of higher-energy excitations, except red-shifted by ∼ 3 K in comparison to the N = 3 energies. For each l, states of different m are now clustered together as was also the case for N = 3 (Fig. 4). These states represent the localized vibrational motion of a single helium atom near the global potential minima, and appear to be associated with the localized helium density found in the path integral calculations of Ref. 25 (Table II).
highly anisotropic helium-benzene interaction potential gives rise to a very structured helium density distribution in the BF frame. In particular, a single helium atom is well-localized at each of the two equivalent global potential minima, above and below the benzene surface. We find a set of collective helium excitations with energies of up to ∼ 23 K above the ground state. Among these excitations, the higher-energy states (> 9 K) can be characterized as a localized excitation deriving from the helium density near the global potential minima, i.e. adsorbed on the molecular nanosurface. The existence of these localized modes is consistent with the localization of a single "dead" helium atom seen in Ref. 25. Helium excitations of lower energies were also obtained, which correspond to collective vibrations of helium atoms delocalized over equivalent sites of lower binding energy, situated near the edges of the benzene surface. Both localized and delocalized excitations can also be further classified by their symmetry with respect to the D 6h (M ) group. The energetic range of these helium excitations is similar to that observed experimentally as vibronic structure in the mass-selective excitation spectra of planar aro-matic molecules in small helium clusters (N ≤ 16). 10,11 It is also similar to the energetic range of the peaks observed in the phonon wings of vibronic spectra of planar aromatic molecules in large helium droplets. 6,7,8,9 Thus, these calculations provide support for the picture of helium atoms adsorbed and vibrating on the molecule surface, with the specific details of the helium motions being determined by the geometry of the surface. In particular, the close relationship between the anisotropy of the molecule-helium interaction and the nature and energetic range of the excitations seen here for benzene implies that these molecule-induced vibrational modes will be very sensitive to the identity of the molecule, possibly even to the extent of providing a spectral "fingerprint" of complex polyaromatic species. Since these calculations have less than one solvation shell of helium surrounding the molecule, they are directly relevant to the recent experimental observations for small cluster sizes. 10,11 Moreover, the analysis of the excitations has been made in terms of single-particle type excitation operators acting on the many-body ground state, and is thus applicable to any number of solvating helium atoms. The fact that the localized excitation associated with the most strongly bound helium density at the global minimum persists in the largest cluster studied here (N = 14), together with the previous identification of localization at that site in larger helium clusters, 25 indicates that these localized excitations will also be present in much larger helium droplets. For the 4 He N -benzene system here we have examined the limit of a single helium adsorbed atom on a single site given by the benzene nanosurface. But given the energetics that we find here, in a more general sense this class of surface-adsorbed vibrations are very likely to be responsible for the structure seen in the phonon wing sideband in experiments using helium droplets. 9 A detailed analysis of the experimental phonon wing data clearly requires molecule-specific calculations with accurate potential surfaces. For example, the region of spatial confinement for an adsorbed helium atom near the global potential minimum is smaller in benzene than in the larger polyaromatic molecules studied experimentally with small numbers of helium atoms in Ref. 11. In a general context, the present results are very encouraging in that they show that the diffusion Monte-Carlo based methodologies can be used to systematically study these excitations as a function of the number of helium atoms, provided that realistic molecule-helium interaction potentials are available. We emphasize here that a quantitative analysis, both for the small cluster vibronic excitations and for the phonon wing structure in electronic absorption spectra, will require detailed knowledge of the molecule-helium interaction potential in both ground and electronically excited states.
Less obvious than the relation to phonon wing structure is whether the excitations of the type studied here are responsible for the splittings on the order of ∼ 1 K in the zero-phonon line, which were observed for some large organic impurity molecules in large helium droplets. This zero-phonon line splitting has been experimentally studied in detail for tetracene, 7,8 and for indole derivatives. 8 In both experimental studies rotational fine structure has been ruled out as being responsible and several different interpretations have been advanced. For tetracene, a molecule possessing a plane of mirror symmetry, it has been suggested that either two inequivalent helium binding sites exist, possibly due to an inhomogeneous solvation environment, or that some kind of two-particle tunneling is manifested. 7 For indole derivatives, molecules possessing aromatic rings but no mirror symmetry and not usually completely planar, the theoretical evidence of localized helium atoms at aromatic rings in larger, superfluid, helium clusters 25 has been used to propose models of helium atoms similarly localized on either side of the aromatic portion of the molecule. 8 The 4 He N -benzene calculations reported here show low-lying excited states (≤ 9 K) for the larger clusters, N = 3, 14. Furthermore, the spectral weights of these low-lying states do decrease significantly relative to the higher-energy states when the localB † -type operators are applied. This shows that the low-energy spectral features we observe the our calculations are due to vibrational motion of helium located near the periphery of the benzene surface, since theB † -type operators specifically de-emphasize states that are not strongly localized near the global potential minima. In contrast to the higherlying localized states, it is more difficult to say whether these low-energy states would remain low-energy with increasing N , since addition of more helium atoms would be expected to significantly change the local details of the helium wave function in these edge regions around the molecule. Thus, it is also more difficult to conclusively claim that these low-energy delocalized modes are responsible for the experimentally observed splittings in the zero-phonon lines.
However, with aromatic molecules of lower overall symmetry like the indole derivatives, it appears reasonable that the effective potential on opposite sides of the aromatic ring may differ, due to the effect of non-symmetric, three-dimensional side chains, resulting in slightly different localized helium densities on either side of the aromatic ring, as proposed in Ref. 8. This suggests that an alternative approach to interpret the splitting of zerophonon lines in experiments would be in the context of impurity spectra in solids, where the impurity molecule can be trapped in structurally inequivalent trapping sites, giving rise to small differences in the spectral shifts of the electronic origin. 54 Quantitative analysis of this kind would also require determination of accurate helium solvation densities around the molecule.
In extrapolating conclusions made from small cluster studies, whether theoretical or experimental, to explain experimental observations made in large droplets, the key point to consider is whether the localization induced by the helium-aromatic molecule interaction is sufficiently strong such that the localized helium excitations characterized here persist with increasing N . As the number of helium atoms increases, collective compressional and surface modes of the droplet are expected to develop, and the coupling to these modes are manifest in the phonon-wing sideband in electronic spectra measured in large droplets. 3 The path integral results for benzene in 4 He 39 provide evidence that the single localized helium atom located on either side of the molecule near the global potential minimum does indeed retain its localized identity in large clusters. 25 This provides a strong argument for the persistence of the single-particle-like localized states whose excitations have been characterized here for smaller clusters (N ≤ 14). An important task for future work is to understand how the discrete spectra for larger polyaromatics vary with the number of helium atoms when this is still less than a solvation shell, 11 as well as what happens to these excitations in much larger helium clusters. A theoretical prerequisite for study of all these excitations is now the development of accurate trial functions for these very anisotropic aromatic systems at larger numbers of helium atoms (N > 14). This will then allow analysis of the transition from localized to collective helium vibrations, as well as identification of any persistent localized modes that may exist at an aromatic molecular nanosubstrate.