Experimental Determination of the Interaction Potential between a Helium Atom and the Interior Surface of a C60 Fullerene Molecule

The interactions between atoms and molecules may be described by a potential energy function of the nuclear coordinates. Non-bonded interactions are dominated by repulsive forces at short range and attractive dispersion forces at long range. Experimental data on the detailed interaction potentials for non-bonded interatomic and intermolecular forces is scarce. Here we use terahertz spectroscopy and inelastic neutron scattering to determine the potential energy function for the non-bonded interaction between single He atoms and encapsulating C60 fullerene cages, in the helium endofullerenes 3He and 4He, synthesised by molecular surgery techniques. The experimentally derived potential is compared to estimates from quantum chemistry calculations, and from sums of empirical two-body potentials.


I. INTRODUCTION
Non-bonded intermolecular interactions determine the structure and properties of most forms of matter. The potential energy function specifies the dependence of the potential energy on the nuclear coordinates of the interacting moieties, within the Born-Oppenheimer approximation 1 . The estimation of potential functions for non-bonded interactions remains an active research area of computational chemistry [2][3][4] . Ab initio methods are capable of high accuracy but are usually too computational expensive to be applied to anything but very small molecular systems. Computational techniques with good scaling properties such as density functional theory (DFT) are generally imprecise for non-bonded interactions, unless customised adjustments are made [3][4][5] . The accuracy of quantum chemistry algorithms is often assessed by seeking convergence with respect to the calculation level, or number of basis functions 2 .
Advances in all fields of science require comparison with experiment. Unfortunately, detailed experimental data on intermolecular potential energy surfaces is scarce. Some information may be gained by comparing crystal structures and energetics with those derived from model potentials 6 . The equilibrium structures, dissociation energies and vibrational frequencies of intermolecular complexes and clusters may be studied in the gas phase and molecular beams [7][8][9][10][11][12] . However these measurements encounter difficulties with control of the local sample temperature, and only provide information on potential minima, and their local properties close to potential minima. Atomic beam diffraction may also provide information [13][14][15] .
An ideal set of systems for the study of intermolecular interactions is provided by atomic and molecular endofullerenes, in which single atoms or small molecules are encapsulated in closed carbon cages [16][17][18][19] . A range of small-molecule endofullerenes is available in macroscopic quantities through the multi-step synthetic route known as "molecular surgery" 20 , including H 2 @C 60 18 , H 2 @C 70 21 , H 2 O@C 60 19 , HF@C 60 22 , CH 4 @C 60 23 , and their isotopologues. Endofullerenes containing noble gas atoms, and containing two encapsulated species, may also be produced 21,[24][25][26][27][28][29][30] . Endofullerenes are chemically very stable, may be prepared in a pure and homogeneous solid form, and may be studied at almost any desired temperature.
At low temperatures, the translational modes (and for non-monatomic species, the internal degrees of freedom) of the endohedral species are quantized. The quantum levels may be probed by a wide range of spectroscopic techniques 31 , including infrared spectroscopy 22,[32][33][34][35][36] , pulsed terahertz spectroscopy 37 , nuclear magnetic resonance 22,29,34,[38][39][40] , and inelastic neutron scattering 22,34,41,42 . When performed at cryogenic temperatures, these techniques reveal a rich energy level structure for the quantized modes of the encapsulated systems 22,32,34,41,43 .The quantum structure has been studied in detail using models of the confining potential, sometimes combined with cage-induced modifications of the rotational and vibrational characteristics of the guest molecule 32,33,35,36,43-56 . There are two main ways to describe the interaction potential between the encapsulated species and the cage. One approach describes the interaction potential as a sum over many two-body Lennard-Jones functions involving each endohedral atom and all 60 carbon atoms of the cage [44][45][46]49,50,52,53,55,56 , sometimes introducing "additional sites" on the endohedral species as well 46,52,53 . One disadvantage of this approach is that the summed potential has an undesirable dependence on the precise radius of the encapsulating fullerene cage. An alternative approach, which we call "model-free", describes the interaction potential as a sum of orthogonal spatial functions 32,33,35,36,43,47,48 . The latter approach makes no assumptions about the cage geometry and is better-suited for a comparison with computational chemistry methods.
In this report, we "go back to basics" by studying the simplest atomic endofullerene, He@C 60  . Terahertz and neutron scattering data is acquired and fitted by a simple quantum-mechanical model consisting of a particle confined by a three-dimensional potential well. This allows us to define a "model-free" atom-fullerene potential, with no assumptions about whether it may be expressed as the sum of many two-body interactions. Although He@C 60 was first made in trace amounts by gas-phase methods 16,17,57 , molecular surgery techniques now provide both isotopologues 3 He@C 60 and 4 He@C 60 in high purity and macroscopic quantities 25,30 . These synthetic advances have made it feasible to perform terahertz spectroscopy and inelastic neutron scattering experiments on solid polycrystalline samples of He@C 60 at low temperature, with good signal-to-noise ratio.
At first sight, He@C 60 is an unpromising object of study by both terahertz spectroscopy and neutron scattering. Since He atoms are neutral, their translational motion is not expected to interact with electromagnetic radiation. Furthermore, both 3 He and 4 He isotopes have small neutron scattering cross-sections, and 3 He is a strong neutron absorber. Fortunately, although these concerns are valid, they are not fatal. The He atoms in He@C 60 acquire a small induced electric dipole through their interactions with the encapsulating cage, and hence interact weakly with the THz irradiation, as in the case of H 2 @C 60 32 . The feeble neutron scattering of both He isotopes may be compensated by a sufficiently large sample quantity.
We compare the experimentally determined potential to estimates from empirical two-body interaction potentials, and from quantum chemistry calculations. Empirical two-body potentials give widely divergent results, even when those potentials are based on experimental helium-graphite scattering data. Møller-Plesset perturbation theory techniques and density functional theory (DFT) methods which explicitly in-clude, or are empirically corrected to account for, dispersive interactions, are shown to provide good estimates for the interaction potential.

A. Sample Preparation
3 He@C 60 and 4 He@C 60 were synthesised using a solidstate process for the critical step, as described in reference 30. The initial filling factors were 30% to 50%. The samples were further purified by recirculating HPLC on Cosmosil Buckyprep columns to remove trace impurities of H 2 O@C 60 . Without this precaution, strong neutron scattering by the hydrogen nuclei interferes strongly with the INS measurements. For THz spectroscopy samples of high filling factor were required to get sufficient signal and were prepared by further extensive recirculating HPLC. All samples were sublimed under vacuum before spectroscopic measurements.

B. Terahertz Spectroscopy
THz absorption spectra were measured with an interferometer using a mercury arc light source and a 4 K bolometer as an intensity detector. The typical instrumental resolution was 0.3 cm −1 , which is below the width of the measured THz absorption lines. The 4 He@C 60 sample had a filling factor of f = 88.2 ± 0.5% while the 3 He@C 60 had a filling factor of f = 97.2 ± 0.5%, as determined by 13 C NMR. The sample pellets were pressed from fine powders of solid He@C 60 . The temperature dependence of the absorption spectra was measured by using a variable-temperature optical cryostat. More information is in the Supplementary Material.

C. Inelastic Neutron Scattering
INS experiments were conducted using the IN1-Lagrange spectrometer at the Institut Laue-Langevin (ILL) in Grenoble. Incident neutrons are provided by the "hot source" moderator of the reactor, resulting in a high flux neutron beam. A choice of three different single crystal monochromators, namely Si(111), Si(311) and Cu(220) are used to define the incident energy of the monochromatic neutron beam arriving at the sample using Bragg reflection. The neutrons scattered by the interaction with the sample enter a secondary spectrometer comprising a large area array of pyrolytic graphite analyzer crystals. The focussing geometry of the secondary spectrometer ensures that only neutrons with a fixed kinetic energy of 4.5 meV are detected by the 3 He detector. INS spectra were recorded in the energy transfer range [5,200] meV for the 3 He@C 60 sample, while it was reduced to [5,60] meV for 4 He@C 60 as the time allowed for performing the latter experiment was reduced.
The powdered samples, with respective mass of 1067 mg for 3 He@C 60 ( f = 45%) and 294 mg for 4 He@C 60 ( f = 40%) were loaded inside an Al foil and further inserted inside a cylindrical annulus before they were mounted at the tip of an orange cryostat and placed inside the IN1 spectrometer beam. The sample temperature was kept around 2.7 K. In order to subtract background and scattering from Al and from the C 60 cage, a blank mass matching sample of C 60 was measured using the same setup and an empty cell was also measured. In order to account for the strong absorption of 3 He@C 60 , a Cd sample was also measured enabling to correct from the incident energy dependent absorption of the sample. The neutron counts in figure 3 were normalized to the incident neutron flux.

A. Terahertz spectroscopy
Terahertz absorption spectra for 3 He@C 60 and 4 He@C 60 at two different temperatures are shown in figure 2. For both isotopologues, the high-temperature spectrum displays a comb of several clearly resolved THz peaks, with the 3 He peaks having higher frequencies than those of 4 He. As discussed below, the combs of THz peaks indicate that the potential energy function V (r) for the encapsulated He does not have a purely quadratic dependence on the displacement r of the He atom from the cage centre. This indicates that the He dynamics is not well-described as a purely harmonic three-dimensional oscillator.
The 5 K spectra in figure 2 display a single peak with partially-resolved fine structure, for both 3 He@C 60 and 4 He@C 60 . These fundamental peaks correspond to transitions from the quantum ground states of He in the two isotopologues. The fine structure requires further investigation, but may be associated with a small perturbation of the confining potential by the merohedral disorder in the crystal lattice. Similar effects have been identified for H 2 @C 60 42 .

B. Inelastic neutron scattering
Inelastic neutron scattering spectra for 3 He@C 60 and 4 He@C 60 are shown in figure 3. The figure shows the difference between the INS of the He endofullerenes and that of pure C 60 . The INS spectra before subtraction are shown in the Supplementary Material. Since C 60 has no vibrational modes below ∼ 250 cm −1 , and the low-energy phonon spectrum cancels precisely for the empty and filled fullerenes, the peaks below this energy threshold are clearly attributable to FIG. 3. Inelastic neutron scattering of He endofullerenes. (a) Inelastic neutron scattering spectra of 3 He@C 60 at a temperature of 2.7 K (blue). (b) Inelastic neutron scattering spectra of 4 He@C 60 at a temperature of 2.7 K (green). In both cases, a weighted difference between the scattering of He@C 60 and pure C 60 is shown, with the weighting factors adjusted for best subtraction of the C 60 background. The short vertical bars indicate the predicted positions of the INS peaks for the quantized He motion under the radial potential energy function specified in table I. The INS peaks are labelled according to the transition assignments in figure 4(b). The peaks above ∼ 250 cm −1 and marked by asterisks are due to scattering from the C 60 cages, whose modes are slightly modified in frequency by the presence of endohedral He. the quantized modes of the confined He atoms. As in the case of THz spectroscopy, the 3 He INS peaks are at higher energies than for 4 He.
The strong features above ∼ 250 cm −1 are attributed to the known vibrational modes of C 60 molecules 58 . Raman studies have shown that the radial vibrational modes of the C 60 cages are slightly blue-shifted by the presence of an endohedral noble gas atom 59 . These shifts lead to imperfect cancellation in the INS difference spectra, causing the dispersion-like features in figure 3 which are marked by asterisks. These subtraction artefacts are much stronger for 4 He than for 3 He, for two reasons: (i) the C 60 vibrational modes are slightly more shifted for 4 He than for 3 He, due to its larger mass; (ii) 4 He has a much lower scattering cross-section than 3 He.

A. Energy levels and transitions
The Schrödinger equation for the confined atom (within the Born-Oppenheimer approximation), is given bŷ where q describes a set of quantum numbers, q = {q 1 , q 2 , . . .}, and E q is the energy of the stationary quantum state. The Hamiltonian operatorĤ is given bŷ wherep is the momentum operator and M is the atomic mass.
In general, the energy levels E q and stationary state wavefunctions ψ q depend strongly on the potential energy function V (r), where r represents the nuclear coordinates of the encapsulated atom (figure 1b). The potential energy of the He atom inside the cage may be described by a potential function V (r, θ , φ ), where r is the displacement of the He nucleus from the cage centre, and (θ , φ ) are polar angles. The C 60 cage has icosahedral symmetry, but may be treated as spherical to a good approximation, at low excitation energies of the endohedral atom. The angular dependence may be dropped by assuming approximate spherical symmetry, V (r, θ , φ ) V (r). We assume a radial potential energy function of the form The energy eigenvalues and eigenstates are given by E n m and ψ n m (r, θ , φ ) respectively. The principal quantum number n takes values n ∈ {0, 1, . . .} with the angular momentum quantum number given by ∈ {0, 2, . . . n} (for even n) and ∈ {1, 3, . . . n} (for odd n) 60 . The azimuthal quantum number takes values m ∈ {− , − + 1, . . . + }. For spherical symmetry, the energies are independent of m, so the energy level E n is (2 + 1)-fold degenerate. The stationary quantum states ψ n m (r, θ , φ ) are given by products of radial functions R n (r) and spherical harmonics Y m (θ , φ ), just as for the electronic orbitals of a hydrogen atom 60 .
The eigenvalues and eigenstates depend on the potential coefficients {V 2 ,V 4 ,V 6 } and the mass of the He atom. The electric-dipole-allowed transitions, which are observed in THz spectroscopy and described by the induced dipole moment coefficient A 1 , have the selection rule ∆ = ±1, see Supplementary Material. There are no relevant selection rules for the neutron scattering peaks.

B. Fitting of the Potential
We treat the V 4 and V 6 terms as perturbations of the quadratic V 2 term, which corresponds to an isotropic threedimensional harmonic oscillator.
The solutions of the Schrödinger equation for the isotropic 3D harmonic oscillator are well-known 60,61 , and are given by: with the fundamental vibrational frequency ω 0 = (2V 2 /µ) 1/2 , where µ is the reduced mass (assumed here to be equal to the mass of the 3 He or 4 He atom, since each C 60 molecule is more than two orders of magnitude more massive than the encapsulated atom, and is also coupled to the lattice). The Schrödinger equation was solved approximately for finite V 4 and V 6 by numerically diagonalizing a matrix with elements given by n m|V 4 r 4 + V 6 r 6 |n m . Since the assumed Hamiltonian retains isotropic symmetry, all matrix elements are independent of the quantum number m and vanish for = and m = m . In practice the matrix was bounded by quantum numbers n ≤ 18, after checking for convergence. The THz peak intensities and peak positions were fitted, as described in the Supplementary Material, allowing numerical estimation of the potential parameters V 2 (or ω 0 ), V 4 and V 6 , and the induced dipole moment amplitude A 1 . The derived eigenvalues were used to estimate the INS peak positions.
The fitting of the potential was performed independently TABLE I. Best fit polynomial coefficients and confidence limits for the radial potential function V (r) = V 2 r 2 +V 4 r 4 +V 6 r 6 and induced dipole function d 1q = 4π/3 A 1 rY 1q (θ , φ ) experienced by the confined He isotopes, see SI.
An energy level diagram for the confined He atoms, marked with the observed transitions, is shown in figure 4 The potentials used in (a) and (d) were used for the fitting of He · · · C scattering data 14 . The functional forms of the potentials and their associated parameters are given in the Supplementary Material. In all cases the He atom was displaced from the cage centre towards the nucleus of a carbon atom. The confidence limits in the structural data for C 60 66 lead to error margins on the empirical curves which are smaller than the plotted linewidths. derived potential curves for 3 He and 4 He, despite the different masses of the isotopes and the very different observed frequencies, attests to the validity of the determination of V (r).

C. Comparison with Empirical Potentials
There have been numerous attempts to model the nonbonded interactions between atoms using empirical two-body potential functions such as the Lennard-Jones (LJ) 6-12 potential, or by more complex functional forms. Suitable functions and parameters have been proposed for the He · · · C interaction 14,15,64,65,67 . Some of the proposed two-body potentials were developed for modelling the scattering of He atoms from a graphite surface 14,15 . Figure 5 compares the experimental V (r) curve with predictions from published He · · · C two-body interaction functions. In each case, the total potential energy V (r) was estimated by locating the He atom a distance r along a line from the centre of the cage towards a C atom, and summing the contributions from all 60 two-body He · · · C potentials. The direction of the He displacement has a negligible effect on the calculated potential curves over the relevant energy range (see Supplementary Material). The derived potentials are very sensitive to the geometry of the C 60 cage, especially its radius R. We fixed the locations of all C nuclei to the best current estimates from neutron diffraction 66 , as follows: Bond lengths h = 138.14 ± 0.27 pm for C-C bonds shared by two hexagons, p = 145.97 ± 0.18 pm for C-C bonds shared by a hexagon and a pentagon, and distance of all carbon atoms from the cage centre R = 354.7 ± 0.5 pm. The width of the curves in figure 5 is greater than their confidence limits, which are dominated by the uncertainties in the structural parameters. Explicit functional forms and parameters for the empirical two-body potentials are given in the Supplementary Material.
The most striking feature of Figure 5 is the wide variation of derived potentials for different two-body interaction models. Of all the proposed two-body potentials, the Lennard-Jones 6-12 potential with parameters given by Pang and Brisse 65 (curve a) provides the best agreement with experiment. The isotropic two-body potentials derived by fitting experimental He/graphite scattering data 14,15 (curves c and d) give poor fits to the experimental He@C 60 potential.

D. Comparison with Quantum Chemistry
The He@C 60 system is too large to be treated at the full ab initio level of quantum chemistry. The availability of an experimental radial potential function V (r) allows the direct evaluation of approximate computational chemistry techniques -not only at the equilibrium geometry, but also for displacements of the He atom from the centre of the C 60 cage.
The radial potential V (r) was evaluated by estimating the energy of a He@C 60 system using a range of computational chemistry algorithms, with the He atom displaced by r from the centre of the C 60 cage. In all cases the locations of the carbon atoms were fixed to the C 60 geometry as determined by neutron diffraction 66 , with the same parameters as used for the evaluation of the empirical potentials. The He was moved on the line connecting the cage centre to a carbon nucleus. The direction of the He displacement has a negligible effect on the predicted potential curves over the relevant energy range (see Supplementary Material). The potentials were calculated using the Psi4 program 68 . The functionals used for DFT were: (i) the B3LYP functional, which is one of the most popular semi-empirical hybrid functionals [69][70][71][72][73] ; (ii) the B3LYP functional including the Grimme D3 empirical dispersion correction with Beck-Johnson damping 5,74 ; (iii) the ωB97X-V functional, which includes a contribution from the non-local VV10 correlation functional and is designed to handle noncovalent interactions 69 . The potential was also calculated using second-order Møller-Plesset perturbation (MP2) theory 2 including empirical spin-component-scaling factors (SCS) 75 . All potential calculations employed a counterpoise basis-setsuperposition-error correction, and converged to a good approximation with the correlation-consistent cc-pVXZ (X=D, T, Q, 5) basis sets 76,77 . More details on the quantum chemistry calculations are given in the Supplementary Material. Some comparisons are shown in Figure 6. Density functional theory with the popular B3LYP functional 69-73 overestimates the steepness of the confining potential, although the correspondence with experiment is improved by including the empirical D3 correction with Beck-Johnson damping 5,74 . DFT with the ωB97X-V functional 69  (SCS) 75 , both give an acceptable correspondence between the calculated and experimentally determined potentials.

V. DISCUSSION
We have showed that the quantized energy levels of helium atoms encapsulated in C 60 cages may be probed by THz spectroscopy and INS, despite the weak interactions of the He atoms with the electromagnetic field and with neutrons. The spectroscopic features were analysed to obtain a detailed potential energy function for the interaction between the encapsulated species and the surrounding cage -an interaction dominated by non-bonded dispersion forces which are hard to estimate experimentally. An excellent correspondence was obtained between the interaction potentials derived from independent 3 He@C 60 and 4 He@C 60 measurements, despite the different peak positions for the two samples.
The experimental V (r) curve was compared with sums of published two-body He · · · C interactions. With few exceptions the summed two-body potentials have a poor correspondence with the experimental result. It is not a great surprise that the interaction of a He atom with a highly delocalized electronic structure such as C 60 is hard to model as the sum of individual atom-atom interactions.
We also compared the experimentally derived interaction potential with those derived by quantum chemistry techniques. This allowed the validation of DFT methods which have been developed to deal with dispersive interactions, including the popular B3LYP functional with the D3 empirical dispersion correction 5,74 , and the ωB97X-V functional which incorporates the non-local VV10 correlation functional and has been parameterised using a training set rich in nonbonding interactions 69 . Møller-Plesset perturbation theory with spin-component-scaling factors 75 also provides a good description of the confining potential of the encapsulated He atoms.
There are small discrepancies between the calculated and observed potentials. However it is not yet known whether the remaining discrepancies reflect the limitations in the quantum chemistry algorithms, or the limitations in the assumptions made when interpreting the experimental data -for example, the neglect of the influence exerted by the encapsulated He atoms on the cage radius. Precise measurements of the He@C 60 cage geometry by neutron scattering or X-ray diffraction are planned.
He atoms are small, have no static dipole moment, and a low polarizability. This makes He@C 60 a relatively easy case for computational chemistry. A stiffer challenge for computational chemistry is likely to be presented by compounds in which the endohedral species is polar, such as H 2 O@C 60 24 and HF@C 60 22 , and by endofullerenes such as CH 4 @C 60 23 , where the fit with the cage is much tighter. Furthermore, the study of systems with multiple atoms or molecules encapsulated in the same fullerene cage 24,26,27 should allow the study of non-bonded molecule-molecule and molecule-atom interactions.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

THz spectroscopy
The THz transmission spectra were obtained using the Bruker interferometer Vertex 80v and an optical coldfinger type continuous flow cryostat with two thin-film polypropylene windows. For the measurements the powder of resublimed He@C 60 samples was pressed under vacuum into a 3 mm diameter hole of a brass frame. The mass and thickness of 3 He@C 60 sample were 28 mg and 2.16 mm and of 4 He@C 60 sample 21 mg and 1.72 mm. The brass frame with the pellet was inserted into a sample chamber with two thin-film polypropylene windows and with a vacuum line for pumping and filling with Helium heat exchange gas. The sample chamber was in a thermal contact with the cold finger of the cryostat. The cryostat was placed inside the Vertex 80v sample compartment. The cold finger was moved up and down by letting the beam through the sample chamber or through a reference hole with 3 mm diameter. The mercury arc lamp, 6 µm Mylar beamsplitter and 4 K bolometer were used to record the transmission spectra below 300 cm −1 . The apodized resolution was typically 0.3 cm −1 or better.
The transmission T r (ω) was measured as the light intensity transmitted by the sample chamber divided by the light intensity transmitted by the reference hole. The absorption coefficient α(ω) was calculated from the transmission T r (ω) through α(ω) = −d −1 ln T r (ω)R −1 corr , with R corr as the amount of radiation reflected from the sample pellet surface and from the windows of the sample chamber. R corr adds to the background absorption but does not affect the intensities of He@C 60 absorption lines.
Spectra presented in the paper were measured at low temperature, 5 K, and at high temperature, 100 K for 4 He@C 60 and 125 K for 3 He@C 60 . The low and high T absorption spectra were treated differently before fitting the peaks with Gaussian functions. The baseline correction by subtracting the slowly changing background was performed on the low temperature absorbance spectra and on the lowest frequency line in the high T spectra. The baseline of the high T spectrum above the lowest frequency resonance was corrected by subtracting the 5 K spectrum from the high T spectrum. The position of a broad line observed at 140 cm −1 in the 5 K spectra is independent of the mass of He and independent of T as it subtracts out in the high T spectra if using the 5 K spectrum as a baseline. Therefore, the 140 cm −1 resonance is not caused by the presence of endohedral He and is likely a C 60 lattice mode.
He@C 60 filling factors from 13 C NMR The 13 C solution NMR spectra of 3 He@C 60 and 4 He@C 60 , in ODCB-d 4 (Sigma-Aldrich), were acquired in order to measure the filling factors of the endofullerenes. The measurements were performed at 298 K and a field of 16.45 T on a Bruker Ascend 700 NB magnet fitted with a Bruker TCI prodigy 5 mm liquids cryoprobe and a Bruker AVANCE NEO console.
The 13 C solution NMR spectrum of 3 He@C 60 is shown in fig. S1 (a), which results in a filling factor f = 97.2 ± 0.5%. The 13 C solution NMR spectrum of 4 He@C 60 is shown in fig. S1 (b), which results in a filling factor f = 88.2 ± 0.5%.

Inelastic neutron scattering 3 He@C 60 INS
The low temperature INS spectra of 1067 mg 3 He@C 60 (f=45%) and 1067 mg C 60 are shown in fig. S2 (a), in red and black respectively. The 3 He@C 60 spectrum is corrected for the neutron absorption by the encapsulated 3 He. The C 60 spectrum is scaled to best match the 3 He@C 60 in the low energy transfer region of the spectrum, in order to obtain optimum background subtraction. In fig. S2 (b) the difference between 3 He@C 60 and C 60 is shown in blue.

He@C 60 INS
The low temperature INS spectra of 294 mg 4 He@C 60 (f=40%) and 293 mg C 60 are shown in fig. S3 (a), in green and black respectively.
The C 60 spectrum is scaled to best match the 4 He@C 60 spectrum, in order to obtain optimum background subtraction. In fig. S3 (b) the difference between 4 He@C 60 and C 60 is shown in dark-green. An experimental artefact, present in both 4 He@C 60 and C 60 spectra, is marked with a dagger ( †).

Quantum mechanics for spherically symmetric potentials
The time-independent Schrödinger equation is given by: where the vectors r 1 , r 2 . . . are the coordinates of the particles, the quantum state is described by the eigenfunction ψ q (r 1 , r 2 . . .), q represents any set of quantum numbers (q = {q 1 , q 2 . . .}) and E q is the eigenvalue (energy) of the eigenfunction ψ q . When applying the Hamiltonian operatorĤ (eq. 2) to the eigenfunction ψ q , one obtains the eigenfunction back and the eigenvalue E q associated with it. The Hamiltonian operatorĤ is given by:Ĥ wherep i and M i are the momentum operator and mass for particle i. In the absence of a potential V , eq. 2 is the Hamiltonian for a free particle. For the cases studied here the potential is not zero, and for endofullerenes V is the confining potential which keeps the endohedral moiety enclosed. The meaning of V is such that a quantum particle i at coordinates r 1 has potential energy given by V (r 1 ).
In the cases of noble gas endofullerenes studied here, there is only one confined particle so the eigenfunction depends only on three spatial variables, ψ q (r 1 ) = ψ q (r, θ, φ) in spherical polar coordinates. The variables are separated into a radial part and angular part: The Schrödinger equation is factorised into a radial (eq. 4) equation and an angular equation (eq. 6), 2 where M is the mass of the particle: u q (r) = rR q (r) (5) The wavefunctions which solve the angular equation are spherical harmonics 2 and are given below: Here P m (x) are the associated Legendre functions and P (x) are the Legendre polynomials. Thus, for the angular part of the Schrödinger equation the solutions are known. is the angular momentum quantum number and m is azimuthal quantum number.

S5
The radial solutions of the Schrödinger equation depend strongly on the potential V (r).
Which means: = 0, 2, ..., n − 2, n; if n = even = 1, 3, ..., n − 2, n; if n = odd Since the potential is spherically symmetric, the 2 + 1 degeneracy of m substates of each state is not broken. The usual n quantum number of a harmonic oscillator is given by n = + 2n r ; where the quantum number n r represents the number of radial nodes that the wavefunction R n (r) posses. Furthermore, eigenstates of different but same n quantum numbers have accidental degeneracy because the potential is purely harmonic.

3D polynomial oscillator
The 3D HO quadratic potential (V 2 r 2 = k 2 r 2 ) is not sufficient to describe the confining potential of He@C 60 . We use the more general polynomial form: where {V 2 , V 4 , V 6 } are polynomial coefficients, assuming spherical symmetry. There are no known analytic solutions of the radial Schrödinger equation (eq. 4) for polynomial potentials. Numerical linear algebra diagonalisation is used to calculate the solutions, which converge given a large enough basis set. An element of the matrix representation of a given HamiltonianĤ (a) is defined as follows: 2

S6
Eq. 18 is simplified by the orthonormality of the spherical harmonics (eq. 19 below): 2 Eq. 19 dictates that the final Hamiltonian is block diagonal, with each different subspace making up the blocks. When using the 3D HO as a basis, with the wavefunctions from eq. 12 as basis states, the matrix representation ofĤ 0 (eq. 10) is diagonal: ψ n m |Ĥ 0 |ψ n m = E n δ nn δ δ mm (20) where E n are the 3D HO eigenvalues (eq. 16) and the Kronecker delta δ ab = 1 if a = b or = 0 if a = b. See fig. S4 for the form of the matrix representation ofĤ 0 . The left-over terms in the potential, V 4 r 4 + V 6 r 6 (eq. 17), are written as: By using eq. 18 all matrix elements ofĤ (4) andĤ (6) are computed analytically, using the Wolfram Mathematica software. Matrix representations ofĤ (4) andĤ (6) (in the 3D HO basis) are seen in fig. S4, for = 0 and n r = 0, 1, ..., 18. The Hamiltonian for the 3D polynomial oscillator discussed here isĤ: squared, χ 2 = [y − f (ω n , {κ})] 2 , the numerical values were substituted for symbols and the Hamiltonian diagonalized numerically. Here f (ω n , {κ}) is the theoretical spectrum with same linewidth and lineshape as the synthetic experimental spectrum; {κ} is the set of fit parameters: Hamiltonian and dipole operator parameters, {V 2 , V 4 , V 6 , A 10 }. As first partial derivatives are zero at the best fit, we used second derivatives to calculate the error margins ∆κ i of the fit parameters: In the spherical approximation the energy does not depend on m. Therefore it is practical to use a reduced basis and reduced matrix elements of spherical tensor operator T kq of rank k which are independent of m and q. 11 The rank of potential spherical operators is k = 0 and the rank of dipole operator is k = 1. In the reduced basis the number of states is smaller by factor 2 + 1 for each . In such reduced basis there are 100 states for n max = 18.
Forcing V 6 = 0 increased the χ 2 , as compared to the fit with all three potential parameters, three times and two times for 3 He and 4 He, respectively.
The calculated energy levels for the best-fit parameters (with n max = 18) are given in table S1 and S2.
S9 Table S1: Translational energy levels of 3 He@C 60 obtained from the fit of 125 K THz absorption spectrum.
Translational energy E, the angular momentum quantum number , and the amplitude squared |ξ| 2 of the main component of eigenstate with the principal quantum number n. The energies are given relative to the ground state. The potential parameters are given in Table 1

Empirical potentials
The confining potentials in this section are obtained by summing 60 two-body potentials between the enclosed Helium and the Carbon atoms constituting the C 60 cage.
To obtain the confining potential V (r), first a 3D structure of C 60 is generated with a given radius and HP & HH bond lengths. A two-body interaction potential U (ρ) between the carbon atoms of the cage and the endohedral Helium atom is chosen. The confining potential at any one point is obtained by summing all the 60 He · · · C two-body interactions. To get the confining potential along an axis, the position of the Helium atom is moved along that axis and then the potential is computed for each position.
where ρ is the interatomic distance between Helium and Carbon, is the potential well depth, and σ is the interatomic distance at which the potential energy is zero. The empirical He · · · C two-body Lennard-Jones (6-8-12) interaction potential, U LJ 6−8−12 (ρ), obtained from He/graphite scattering experiments (ref. 15), is defined in eq. 30 below, with parameters given in table S4.
where ρ is the interatomic distance between Helium and Carbon. The equilibrium position (minimum energy) is denoted by ρ m . The parameter s specifies the ratio of the two attractive terms (6 and 8) at ρ m . Lennard-Jones 6-12 potential (ref. 16) The He · · · C two-body Lennard-Jones interaction potential (U LJ (ρ)), from ref. 16, is defined in eq. 31 below, with parameters given in table S5.
where ρ is the interatomic distance between Helium and Carbon. is the potential well depth and σ is the interatomic distance at which the potential energy is zero. Another two-body He · · · C interaction potential is the modified-Buckingham (mB) potential, U mB (ρ), from ref. 17, which is obtained from the MM3 molecular mechanics program. 18 It is defined in eq. 32, with parameters given in table S6.

Comparison of two-body interaction potentials
All the two-body He · · · C interaction potentials U (ρ) used in this study are plotted in fig. S5. Figure S5: Helium · · · Carbon two-body interaction potentials U (ρ) against interatomic distance ρ. (a, green) Lennard-Jones 6-8-12 potential with parameters from reference 15; (b, blue) Modified Buckingham potential (as implemented in the MM3 program 17,18 ); (c, orange) Lennard-Jones 6-12 potential with parameters from reference 16; (d, red) Lennard-Jones 6-12 with parameters from reference 15. The potentials used in (a) and (d) have been used for the fitting of He/graphite scattering data. 15

Anisotropy of confining potentials
All the He@C 60 confining potentials V (r) computed by summing the 60 He · · · C two-body interaction potentials U (ρ) are slightly anisotropic, since C 60 is not a perfect sphere. However, in all cases the anisotropy was negligible S12 in the energy range probed by our measurements; changing the direction along which the He is moved makes a very small difference to the confining potential. Figure S6 shows an example for the Lennard-Jones potential in eq. 31 and ref. 16. Similar plots are obtained for all the other potentials. Figure S6: He@C 60 confining potentials V (r); obtained by summing the 60 He · · · C empirical Lennard-Jones two-body interaction potentials U (ρ) (eq. 31 and ref. 16

Computational chemistry
We calculated the He@C 60 radial potential V (r) using Psi4. 19 Psi4 uses density-fitting (DF) to achieve favorable scaling with respect to system size by casting expensive electron repulsion integrals into linearly scaling auxiliary basis. Coupled-cluster theory 20,21 [CCSD, CCSD(T)] is generally considered to be state-of-the-art method for interaction energy calculation but is practically not applicable for systems of this size. We calculated the potential with density functional theory (DFT) and Møller-Plesset perturbation (MP2) theory. 20 The DF procedure was used for both the self-consistent field (SCF) reference energy and the MP2 energy calculation. By default, MP2 calculates same-spin and opposite-spin contributions to the correlation energy with different accuracy. The accuracy of the calculated energy was improved by applying the empirical spin-component-scaling factors (SCS). 22 The DFT potential was calculated with the ωB97X-V, 23,24 and B3LYP, 25-28 hybrid functionals. The first functional is designed to handle non-covalent interactions with in-build contribution from the non-local VV10 correlation functional. 24 The latter functional is one of the most popular semi-empirical hybrid functional. 24 B3LYP cannot by default handle dispersion interaction and the Grimme D3 empirical dispersion correction with Beck-Johnson damping 29,30 was applied to the functional.
The potential was calculated as a function of the Helium distance from the centre of the cage, with C 60 structural parameters from neutron diffraction measurements 14 (see the main paper or section above). The Helium was moved along the axis between two inversion-related carbons. The basis set convergence was tested with the Helium at 50 pm and 100 pm from the center of the cage. The calculations converged to a good approximation at X= 5 with the correlation-consistent cc-pVXZ (X=D, T, Q, 5, 6) basis sets. 31,32 Figures S7  and S8 show the calculated potentials as a function of basis set size when the Helium atom is at 50 pm and 100 pm from the center of the cage respectively. The auxiliarity basis sets used by the DF algorithm were the defaults chosen by Psi4. The counterpoise basis-set-superposition-error correction was applied in all the calculations, 33 and all calculations were run in C1 symmetry. The dependence with respect to the direction the Helium is moved was investigated with MP2. The potential was calculated with the Helium moved along the three different symmetry axes of C 60 and an axis between two inversion-related carbons. The difference between the calculated potentials was found to be negligible, with the Helium in the range of 0 pm to 100 pm from the center of the cage. The potentials shown in the main paper were calculated in cc-pVQZ level with the Helium moved along the axis between two inversion-related carbons. S13 Figure S7: Quantum chemically calculated potentials (with respect to the minimum) as a function of the basis set size with Helium at 50 pm from the center of the cage. Figure S8: Quantum chemically calculated potentials (with respect to the minimum) as a function of the basis set size with Helium at 100 pm from the center of the cage.