Exploring the parameter space of an endohedral atom in a cylindrical cavity

Endohedral fullerenes, or endofullerenes, are chemical systems of fullerene cages encapsulating single atoms or small molecules. These species provide an interesting challenge of Potential Energy Surface (PES) determination as examples of non-covalently bonded, bound systems. While the majority of studies focus on C$_{60}$ as the encapsulating cage, introducing some anisotropy by using a different fullerene, e.g., C$_{70}$ can unveil a double well potential along the unique axis. By approximating the potential as a pairwise Lennard-Jones (LJ) summation over the fixed C cage atoms, the parameter space of the Hamiltonian includes three tunable variables: $(M,\varepsilon,\sigma)$ representing the mass of the trapped species, the LJ energy, and length scales respectively. Fixing the mass and allowing the others to vary can imitate the potentials of endohedral species trapped in more elongated fullerenes. We choose to explore the LJ parameter space of an endohedral atom in C$_{70}$ with $\varepsilon\in$ [20cm$^{-1}$, 150cm$^{-1}$], and $\sigma\in$ [2.85\r{A} , 3.05\r{A}]. As the barrier height and positions of these wells vary between [1cm$^{-1}$, 264cm$^{-1}$] and [0.35\r{A}, 0.85\r{A}] respectively, using a 3D direct product basis of 1D harmonic oscillator (HO) wavefunctions centred at the origin where there is a local maximum is unphysical. Instead we propose the use of a non-orthogonal basis set, using 1D HO wavefunctions centred in each minimum and compare this to other choices. The ground state energy of the X@C$_{70}$ is tracked across the LJ parameter space, along with its corresponding nuclear translational wavefunctions. A classification of the wavefunction characteristics, namely the prolateness and ``peanut-likeness'' based on its statistical moments is also proposed.


I. INTRODUCTION
Endohedral fullerenes or endofullerenes are molecular complexes where a fullerene cage traps an atom or molecule(s).The development of the technique known as "molecular surgery" 1,2 has allowed for controlled synthesis of not just noble-gas or metallo-endofullerenes, but also light molecular endofullerenes including H 2 @C 60 3 , H 2 O@C 60 4 and CH 4 @C 60 5 .The fullerene cage almost completely shields the guest species from the environment, and these isolated species epitomise behaviour dominated by quantum phenomena.
These strong quantum effects arise due to the encapsulating potential quantising the translational motion, and the coupling of these modes with with the rotational states of the trapped molecule 6 .The spectroscopic fine structure associated with this translational-rotational coupling has given rise to numerous theoretical studies focusing on H 2 @C 60 7-13   and H 2 O@C 60 [14][15][16][17][18][19][20] .In order to understand the coupling, the translational-rotational eigenstates must be found which requires construction of a Potential Energy Surface (PES).For both of these major systems, the PES is usually generated by assuming a pairwise-additive Lennard-Jones (LJ) form 6,21 , with energy scale ε and length scale σ , as this is well suited for modelling dispersion interactions 22 , with the parameters fitted to match the spectroscopic data 7 .Calculating an accurate ab-initio surface is computationally intensive because of the large number of atoms and interactions dominated by dispersion effects, and so only the PES of HF@C 60 has been acquired using these methods 23 .
While the majority of studies have focused on C 60 as the en-  capsulating fullerene due to its high I h symmetry and wealth of experimental data, this is not the only available choice.The next most abundant fullerene, C 70 , has its symmetry reduced by elongating the buckyball along one axis and its cavity can contain larger or more molecules including (H 2 ) 2 @C 70 24 and (H 2 O) 2 @C 70 25 .Elongating the cage even further in a similar manner produces larger fullerenes or fullertubes in the D 5d/h series which have been recently synthesised [26][27][28][29] .An important feature which arises in the PES of these lower symmetry fullerenes/fullertubes is a symmetric double well potential along the unique axis.Traversing this series for the fixed LJ parameters of Ne@C 70 given in Table I towards larger cages pulls the minima further apart, and increases the barrier height.The first row of parameters have been used in the calculation of the endohedral vibrational levels Ne@C 70 21,22 , but the latter set were used in the calculation of the potential of noble gas endohedral and exohedral complexes with C 60 and C 70 30,31 and are more commonly used.
In order to garner some intuition of the behaviour of a species trapped in a cylindrical cavity created by a fullerene from this series, knowledge of the PES is required.While either the endohedral species or the encapsulating fullerene can be altered, the limited choice in cages is restrictive.The alternative we choose to adopt is to fix the fullerene cage and arXiv:2306.08491v2[physics.chem-ph]20 Nov 2023 vary both LJ parameters in order to more flexibly change the nature of the double well.This allows for more careful tuning of the barrier height and positions of the minima, similar to what is seen when elongating the fullerene along the D 5d/h series.The geometries of the fullerenes C 60 , C 70 , C 80 , C 90 , C 100 , C 110 , and C 120 28,32,33 compared in Figure 1 were oriented to make the unique axis the z direction and symmetrised using IQMol 34 , and their symmetries validated 35,36 using in-house Python software.Adding an extra belt of C atoms pulls the minima further apart, but also increases the barrier height.Given the form of the potential, these quantities are quite sensitive to the size of the cage, and so changing the C -C bond length can also alter them.Changing the LJ parameters (even though the mass is fixed) can also correspond to changing the endohedral species.The LJ parameters for both the noble gas 30,31,37,38 and metal [39][40][41] -fullerene interactions can take differing values depending on the environment where they were calculated.
The specific LJ parameters of the endohedral complex can be found by either considering the atomic polarisabilities and susceptibilites 22 , or by minimising differences to experimental data 7 .
When the parameters are found by minimising the difference to the experimental data, there can be multiple solutions for them which further justifies the choice of allowing the LJ parameters to take a range of values.Varying the LJ parameters instead of just the fullerene cage allows the barrier height B and minima location z min to change smoothly.This provides more control in tuning these quantities, along with the zeropoint energy (ZPE) as its relation to B partitions the parameter space into two main regions: B > ZPE and B < ZPE.
In the region where B < ZPE, which is the case for Ne@C 70 21 , the atom can "tunnel" between both minima and the simple harmonic oscillator (HO) basis set used within the discrete variable representation (DVR) calculations is per-fectly suited to this region.However, as the minima get pulled apart, and as the barrier height increases as seen in Fig 3, using these functions centred at the origin in the finite basis representation (FBR) becomes unphysical as the potential is a local maximum.Instead of this we propose and benchmark an alternative basis with HO functions centred at both minima.Although this basis is non-orthogonal, the added flexibility can recover results in the usual case of B < ZPE, and we show it is much better suited when the alternative scenario of B > ZPE arises.These two regimes will lead to different forms of ground state wavefunction: one that has a single maximum at the centre of the cage, or one that exhibits two maxima above the minima in the potential.This classification is discrete and assignment is not necessarily obvious, hence a system based on the statistical moments of the wavefunction density is proposed.This description is continuous and may correspond to some experimental observables.
Alongside information about the ground state of the system, diagonalising the Hamiltonian also converges the energies of excited states which can be used to calculate transition energies.These can be discerned spectroscopically, and when compared to the values calculated by this technique, the LJ parameters that best describe the X -C interaction can be ascertained.
In this paper we consider a generic system of X@C 70 where X represents a single particle, e.g. a noble gas atom or the centre of mass of a small molecular species, and vary the LJ parameters to explore a range of barrier heights and minima locations which imitate what could be encountered with a fullerene further down the D 5d/h series.The behaviour of the ground state is tracked across a region of this parameter space, and the differing shapes of the translational wavefunction are classified.Fundamental transition energies of X@C 70 are also calculated across the LJ parameter space, alongside excited state energies of X@C n , where n ∈ [70, 80, 90, 100] for a fixed set of LJ parameters.The relevant background theory on the diagonalisation of the Hamiltonian and wavefunction analysis is discussed in Section II, with the results and discussion in Section III.Section IV contains the conclusions drawn from the systems studied and outlines possible future work, both theoretical and experimental.

A. Choosing a Translational Basis Set
Using the 6-12 pairwise LJ potential, the Hamiltonian for X@C 70 in atomic units can be written as where M is the effective two-body reduced mass of the endohedral complex, (ε, σ ) are the energy and length scale of the X -C interaction 21 .The sum is taken over all fixed positions of the carbon atoms, with r i being the distance between the endohedral atom and the ith carbon atom.Due to the anisotropy in the cage, this is treated using a 3D finite basis representation (FBR) constructed as a direct product of 1D HO wavefunctions expanded from the origin where q r = √ αr is a scaled dimensionless coordinate, and r represents a single Cartesian coordinate.Instead of working in the FBR, the Hamiltonian can be diagonalised within a DVR 42 which requires choosing both an underlying basis set, as well as the set of points where these functions are centred.Choosing the same basis functions for the DVR, this is isomorphic to the FBR as they can be interconverted using a unitary transformation, with the advantage being that the potential energy matrix in the DVR is diagonal.Traditionally, the DVR points are chosen by selecting an interval for each coordinate alongside an energy cutoff and only the points where the potential lies under this are kept.Alternatively, a Potential-Optimised DVR 43 (PODVR) scheme can be considered where the scaling factor for the dimensionless coordinate q r is determined from the potential.This can be chosen by either taking the second derivative of the potential at the minimum and using the relation α = √ k eff M, or by leveraging the variational principle and using the ground state harmonic oscillator as a trial wavefunction in 1D By minimising E(α), a more appropriate effective force constant is recovered which corresponds to a higher zeroth order frequency.This is because the procedure accounts for the higher degree nature of the potential.This strategy is required for the unique axis even when B is small as the origin is a local maximum which would result in unbound basis states.

B. Tackling the Double Well
The C 70 cage is oriented such that the z-axis is the unique axis, with the x-axis aligned to be one of the C ′ 2 rotation axes.As B rises, the aforementioned PODVR method starts to break down as the depth of the well weighs down the effective parabola decreasing α.This necessitates the use of more basis functions (or equivalently, quadrature points), more of which are placed towards the edge of the cage, and even outside it.The potential in this region is of the order of 1,000,000 cm −1 and this leads to unstable solutions.This can be tackled in two different ways, either from a DVR perspective or from the FBR.
From the DVR viewpoint, accounting for the origin being a local maximum, a possible choice of basis set is to use sinc functions 44 .The advantage is that this is an orthogonal basis set, and the kinetic energy matrix elements are known analytically.However, these functions take no account of the nature of the potential and have equidistant sampling points.The further the minima are located from the origin, or as B rises, the more quadrature points that will be needed.An alternative, ignoring the fact that the origin is a local maximum, is to revert to the traditional DVR standpoint and impose the range of sampling points, instead of optimising the scaling parameter.Using the 1D HO basis functions as in Equation ( 2), this dictates an effective harmonic potential and fixes the scaling parameter.This maintains the traditional advantages of the DVR, but (ignoring truncation effects) the double well may imply that the energy doesn't always decrease with increasing number of basis functions.
From the FBR perspective, both minima in the double well can be acknowledged with basis functions expanded around them instead of the local maximum origin.Using 1D HO wavefunctions along the unique direction, the scaled dimensionless coordinate and basis functions can be written as where a is the location of the minimum of the potential.A similar optimisation choice for α as in Section II A exists, but here it is calculated by taking the second derivative at the minimum.As these functions are confined in a single minimum, generating delocalised wavefunctions that span both minima is achieved by taking in-phase and out-of-phase combinations of equivalent basis functions where the kets are defined as in Equation ( 6) and the phase factor chosen to maintain the expected nodal structure with increasing quantum number.This symmetrised double minimum basis set is somewhat reminiscent of the orbitals and configuration basis used when considering the interaction of two H 2 molecules within C 70 45 .As this basis is adapted from the potential, fewer points are expected to be needed to converge the eigenstates.
On the other hand the Hamiltonian cannot be simply manipulated such that these functions form an eigenbasis and so matrix elements have to be calculated explicitly.This can be done using Gauss-Hermite quadrature, further justifying this choice of basis set.The Hamiltonian in Equation ( 8) is split up into four terms: ĥ0 x with eigenfunctions |n x ⟩ and frequency ω x , ĥ0 y with eigenfunctions |n y ⟩ and frequency ω y .These chosen to be 1D HO operators, frequencies and wavefunctions of the form given in Equation (2).kz is the kinetic energy operator along the z direction and ∆V is the potential not accounted for in the definitions of ĥ0 x and ĥ0 y .The Hamiltonian and overlap matrix elements have the following forms Ĥ = ĥ0 x + ĥ0 y + kz + ∆V (8) Depending on the separation of the two minima, the overlap integral ⟨m − |n + ⟩ could be non-zero making the basis set non-orthogonal.Carelessly increasing the basis set size will lead to some linear dependence within the basis set.This overcompleteness arises when states with zeroth-order energies in the z direction larger than B start being included within the basis.The energies are given by (n z + 1 2 )ω z , where ω z is the frequency within the well.These functions are very diffuse and taking linear combinations of these eventually leads to overcompleteness.This is alleviated by projecting the Hamiltonian and overlap matrices into a linearly independent subspace before diagonalising.

C. Wavefunction Classification
In the regime B < ZPE, the ground state wavefunction is ellipsoidal.This shape can be described by the length of its semi-axes, but this requires selecting an isosurface to measure from.An alternative that can be calculated is the standard deviation in each direction.Leveraging the isotropy in the x and y directions, the prolateness of the spheroid, ς 2 can be calculated by considering the ratio of variances in the z and x directions where the last equality follows as µ z ≡ ⟨0|z|0⟩ and µ x ≡ ⟨0|x|0⟩ are zero by symmetry.As B and ZPE become more comparable in size, and eventually when B > ZPE, the wavefunction along the z direction moves away from being singlypeaked at the origin to having two equivalent maxima.This moves the wavefunction away from an ellipsoid, to a more peanut-like shape.The simplest way to classify this would be to just consider σ 2 z , which would be large when there is more density away from the mean.Although this will be the case for a double-maxima wavefunction, it could be generated fictitiously with a shallow and very diffuse single peaked ground state.Instead, the kurtosis of the normalised state can be calculated as This accounts for both the mean and variance of the distribution along the unique direction and measures the density at µ z ± σ z .It can be thought of as the relative weight of the shoulder of the distribution compared to the mean and tails.For a pure Gaussian distribution, the kurtosis is strictly equal to 3. A bimodal distribution which corresponds to a double peaked wavefunction has a kurtosis in the interval [1,3).Together, these two statistics describe the deviation of the shape of the ground state from an idealised spherical wavefunction which is the shape for an atom inside C 60 46 .The variance ratio measures the ellipsoidal nature and the kurtosis gives a measure of the peanut-like shape.

A. Computational Details
Along the x and y directions, 12 basis functions of the form shown in Equation (2) were used, with the scale factor determined by minimising Equation ( 4) with the corresponding one-dimensional Hamiltonian.For the unique direction, the choice between sinc DVR functions, HO DVR functions, or basis functions of the form shown in Equation ( 7) was determined by the convergence of the ground state energy of a fixed Hamiltonian.
A one-dimensional potential was constructed by considering V (0, 0, z) for the potential given in Equation ( 1) for two different sets of LJ parameters, one which corresponds to the zero-point energy being similar to the barrier height and one where the barrier height is larger than the zero-point energy.
The ground state energy is tracked against N representing the number of quadrature points used.For both DVR techniques this is equivalent to the number of basis functions used, whereas the symmetrised double minimum basis set used N 2 functions.Both DVR basis sets were constructed by using a fixed interval of points, with |z| < 2Å.
In both cases (ZPE ≈ B and ZPE < B), the sinc DVR basis set fails to converge to the ground state energy.In the former, it is much tighter and it appears that it will converge to the same energy but this is not as apparent in the latter.While both the HO DVR functions and symmetrised double minimum basis set converge the energy to sub milli-wavenumber accuracy, they achieve this in vastly different ways.When ZPE ≈ B, the symmetrised double minimum basis set is flat, indicating that it had already converged in the smallest basis set size, but when ZPE < B, there is a small increase in energy moving from 24 points to 28 points, signifying the integrals had not fully converged.The HO DVR functions on the other hand (alongside the sinc DVR functions) have a fluctuating behaviour during the convergence.In the simpler polynomial oscillator, which lacks any energy barrier, one would expect the energy to decrease monotonically with increasing basis set size.Due to the presence of the double well, this is not the case and there is oscillatory convergence.
As well as the ground state energy, another important quantity is the energy gap between the first excited state and ground state as it can be measured spectroscopically, corresponding to a tunnelling splitting.Considering that frequency of transitions is on the order 21 of 10 1 cm −1 , it is prudent to converge the energies and frequencies tightly.In both regions the sinc DVR basis frequency is larger than what is seen for the HO DVR functions and symmetrised double minimum basis set.The HO DVR transition frequency also suffers from oscillatory convergence in both regions, as it does for the ground state energy.The symmetrised double minimum basis set transition frequency also shows an initial rise as more quadrature points are needed to converge the higher energy state.
We choose to use the symmetrised double minimum basis set for our calculations, as this takes into account features of the potential, does not require arbitrarily choosing the interval of sampling points and does not have oscillatory convergence of energy levels.We place 12 functions in each minimum for the z direction, giving an overall basis size of 3456 functions.The Hamiltonian and overlap matrices were constructed using Equations ( 9) and (10).The basis set was then canonically orthogonalised 47 by discarding linear combination of basis functions whose overlap eigenvalues were negligible.

B. Potential Energy Surface
Of the three variable parameters in Equation ( 1), M was fixed to be the two-particle reduced mass of Ne@C 70 .Given the form of the Hamiltonian, the quantity Mε = const can be interpreted as an effective energy scale and therefore only one needs to vary.The LJ parameter space, (ε, σ ) chosen to explore was ε ∈ [20cm −1 , 150cm −1 ] and σ ∈ [2.85Å, 3.05Å].These parameters were chosen as they ensured that the points for Ne@C 70 given in Table I were within it, and that the barrier height, B, along the z direction varied over a substantial range: [1cm −1 , 265cm −1 ].This allows for regions where ZPE B < 1 and ZPE B >> 1.The location of minima also spanned a wide range: [0.35Å, 0.85Å] which leads to varying amounts of linear dependence within the basis set.
Allowing for a larger LJ parameter space can also introduce a double well in the x and y directions.These can be tackled with the same methodology outlined in Sections II B position of minima, in the unique direction as defined in Figure 1 of the PES tracked across the LJ parameter space.and III A, but are not considered for the scope of this article.For each point in the LJ parameter space where the Hamiltonian was constructed and diagonalised, the PES was linearly shifted such that its minimum was set to zero.
The location of the minima shows only σ dependence, and as this increases, the minima slowly converge towards the origin.For the 1D LJ potential r eq ∝ σ and inside this fullerene cage this corresponds to pushing the endohedral species further away from the C atoms and closer towards the centre.Once the minima merge together, increasing σ further has no more effect as the centre of the cage is the furthest the species can be from the cage atoms.On the other hand, the barrier height shows dependence on both LJ parameters.While the relation with ε is intuitive, increasing with larger ε as raising the energy scale deepens the wells; increasing σ reduces the barrier height.By pushing the minima closer together, they must eventually coalesce, eradicating the barrier completely.Due to this effective repulsive nature of increasing σ , once the minima merge the double well doesn't reappear.
The difference in 1D potentials throughout the LJ parame- x Cartesian directions.Solid line indicates V (0, 0, z) whereas dashed is V (x, 0, 0).V (0, y, 0) is not shown as for this region it is almost identical to V (x, 0, 0), differing only when the endohedral species is closer to the cage edge.The red curves have LJ parameters (ε, σ ) as (60, 2.98) and the blue curves have (140, 2.88).
ter space can be seen by considering two sets of parameters.
In the x (and y) directions, apart from the minimum having a different value, that of the barrier height, they both resemble potentials of an even polynomial oscillator as seen in Fig 6 .However, along the z direction, not only are the minima found at different positions, namely 0.57Å and 0.79Å for the red and blue curves the barrier height also differs significantly changing from 30cm −1 to 200cm −1 respectively.

Energy
The translational ground state of X@C 70 was tracked across the aforementioned LJ parameter space and can be seen in Fig 7 .The zero-point energy shows a strong dependence on ε, with the ZPE getting larger as ε increases.Increasing ε corresponds to a more steeply growing potential as the endohedral atom moves further away from the origin.This increase in the effective force constant correlates with a larger frequency and The ZPE also increases with increasing σ , but this is a much weaker dependence than seen for ε.Increasing σ pushes the two minima closer together and lowers B, but this is compensated by the faster growing nature of the potential in the region |z| > z min .As the size of this region increases, this acts as an increase in the effective force constant, pushing the ZPE to be larger.
The red contour indicates where ZPE = B and partitions the LJ parameter space into two regions: above it, where ZPE > B, and below, where B > ZPE.Moving further away from this contour in this space, the more dominant ZPE is above the contour and B is below it.As B dominates, the particle is more likely to be localised at the minima and the less likely it is to tunnel between them.As the ZPE dominates, the more likely the particle is to be found in the centre as the double well is not a significant effect.

Moments
Across the LJ parameter space, the standard deviation σ z , the kurtosis κ z , and standard deviation ratio ς as defined in Equations ( 11) and (12) in Section II C were calculated.The standard deviation and kurtosis were tracked throughout the parameter space and shown in Fig 8 .The standard deviation lies in the interval [0.24Å, 0.81Å], and this value is always smaller than the position of the minima seen in Fig 5b .Acknowledging that ⟨z⟩ = 0 by symmetry this indicates that the endohedral atom is more likely to be located between the two minima than closer to the cage.While the overall skewness of the wavefunction must also be zero due to symmetry, when the wavefunction exhibits two maxima each of those is individually skewed towards the centre of the cage.This is due to the potential experiencing a local maximum at the origin, seen in Fig 6, instead of growing steeply towards the cage.
The standard deviation increases with decreasing σ and with increasing ε.Both of these effects change the nature of the minima in the potential, by either pulling them apart and further away from the centre or by accentuating their depth.The influence of the minima is to drag wavefunction density towards themselves and away from the origin which leads to an increase in the standard deviation.
As mentioned in Section II C, the standard deviation cannot distinguish between a wavefunction that exhibits a single maximum at the origin or two maxima near the minima in the potential.This is determined by the kurtosis, defined in (12).For this system, the kurtosis lies in the interval [1.07, 2.58] implying all these wavefunctions represent a platykurtic dis-  tribution, one where the kurtosis is lower than a Gaussian.The kurtosis decreases with increasing ε and decreasing σ .Varying the LJ parameters in this way pulls the minima apart and emphasises their depth.This drags more wavefunction density away from the origin, in the region |z| > σ z which puts more density in the tails and this implies a decrease in kurtosis.
To fully classify the 3D wavefunction, instead of just considering σ z , the quantity ς can be used instead as this is the ratio of standard deviations in the z and x directions which gives a measure of the prolateness of the ellipsoidal wavefunction, as defined in Equation (11).This quantity increases with increasing ε, as the weight of the minima in the potential increases, pulling density towards themselves, and it decreases with increasing σ , as the minima get closer together.The prolateness is considered alongside κ z , defined in Equation ( 12), which gives an indication of how much the idealised singlepeaked Gaussian wavefunction has been squashed, pushing density further away from the origin.
Four paths were chosen to exhibit these statistics on the ground state wavefunction for the X@C 70 system.The first two were for fixed ε and fixed σ , which isolate the effect of a single LJ parameter on the shape of the wavefunction.A third path is one where both LJ parameters vary, along the line σ = 3.08 − 0.0015ε in order to check whether the effects of a single LJ parameter compound or cancel each other.The last path is along the ZPE = B contour, as the range of these two statistics can determine which of them is a good measure of double-peakedness.
Keeping either ε or σ fixed and varying the other LJ parameter, ς lies in the interval [2.45,6.05].With increasing σ and decreasing ε, the prolateness ς decreases indicating a move to a more spherical wavefunction.This effect is compounded when both quantities are varied as ς spans a larger range.Along B = ZPE, the dependence of ς on both LJ parameters remains the same but the range is tightened.
The kurtosis on the other hand has the opposite dependence on the LJ parameters than the standard deviation ratio.The range spanned by varying ε is smaller than when varying σ as changing the latter has a stronger effect on whether the two maxima in the wavefunction eventually coalesce.This is the same behaviour as the two minima in the potential; as they get closer together the wavefunction moves towards a single peak.Once again both these effects stack up when both parameters are varied as κ z approaches its largest value of 2.58 which is only slightly smaller than for a pure Gaussian.Along the ZPE = B contour, as for ς , the range for κ z is very narrow.
The reason to use κ z as a more informative measure on the type of wavefunction becomes apparent when looking along the ZPE = B contour where ς is in the interval [3.76,5.56]but κ z is in the interval [1.17,1.31].Only when both LJ parameters are varied does the kurtosis reach values above 2.2 illustrating the move to a single peaked wavefunction.These statistics vary smoothly as the LJ parameters are varied, indicating the wavefunction shape also smoothly changes throughout the parameter space.The quantity ς changes to show how stretched the ellipsoidal shape is along the unique axis, whereas κ z illustrates the deformation of the idealised Gaussian, single maximum wavefunction into two equivalent peaks.This occurs by  effectively squashing the single peak of a standard Gaussian and stretching it outwards such that the maxima occur near the potential minima.This changeover occurs around a kurtosis of 2.2.As these statistics completely define the wavefunction shape, this can then be used to re-construct the PES and the LJ parameters of the X -C interaction can be discerned.

Nuclear Orbitals
In order to give more tangible meaning to the wavefunction statistics described in Section III C 2, the wavefunctions themselves can be visualised.This gives more of an intuition as to what exactly ς and κ z correspond to.Along the ZPE = B path, there is not too much discernible change in the wavefunction contours seen in Fig 10 .As ε and σ increase, the wavefunc- ), and (d) (150cm −1 , 2.951Å).The contours are taken at 1%, 10%, 25%, 50%, 75%, 90%, and 99% of the maximum amplitude of the wavefunction.
tion becomes more contracted.This is weakly due to σ as this pushes the minima closer to the centre, from 0.78Å down to 0.66Å.The stronger dependence is on ε, as this dictates the depth of the well and steepness of the potential, and the larger it is the more tightly bound the wavefunction will be.
The stark difference in wavefunction shapes can be seen along the line which varies both LJ parameters, as illustrated in Figs 11 and 12.As ε increases, and σ decreases on the line σ = 3.08 − 0.0015ε, there is a smooth deformation of the 1D wavefunction from a more typical Gaussian, which becomes squashed and stretched along the line into two shoulder peaks and eventually into two very distinct, localised maxima.Con-sidering the 2D contour plots, when (ε, σ ) are (40, 3.02), the wavefunction exhibits a typical, simple ellipsoidal character.Moving to (70, 2.975), while maintaining a mostly overall ellipsoidal shape, there are clearly two maxima present in this wavefunction.When the LJ parameters are (110, 2.915), the two maxima are very prominent and the wavefunction shows a significant peanut-like shape.Moving to (150, 2.855), this is the extreme case showing two isolated lobes of wavefunction density, indicating two prominent maxima in the minima of the potential, with no discernible wavefunction amplitude connecting these regions.
functions.Between the first two wavefunctions ς is fairly similar, only increasing slightly but the kurtosis drops quickly from 2.5 down to 2.0.For the latter two, the kurtosis is nearly the same, but ς increases from 5.5 to 6.8.The decrease in kurtosis and increase in standard deviation ratio both correspond to moving wavefunction density away from the origin.The latter is more sensitive initially on the move from a single peak wavefunction to having two maxima.Once these maxima are prominent, putting more density there has less significant effect on the kurtosis as this is lower bounded by 1, but σ z keeps on increasing.

D. Excited States
In addition to converging the ground state, excited state energies were also calculated.Here, instead of varying the LJ parameters which were fixed to be (43.79cm−1 , 3.03Å), the parameters for Ne@C 70 as shown in Table I, the fullerene cage was varied along the D 5d/h series from C 70 to C 100 .The lowest 50 eigenstate energies were calculated and assigned quantum numbers (L, m, n z ), analogously to Refs 7 and 21, where L and m are the angular momentum quantum numbers of the 2D isotropic harmonic oscillator in the xy plane and n z the  quantum number of the oscillator in the z direction.These quantum numbers are assigned by analysing the nodal structure of the wavefunctions, with n z being the number of nodes along the z direction.By considering the L−m 2 radial and m angular nodes, the angular momentum quantum numbers L and m can also be assigned.In Figure 13, the eigenstates of X@C 70 are in black, X@C 80 in red, X@C 90 in green and X@C 100 in blue.For fixed (L, m), going up in energy corresponds to increasing n z .For a fixed L and n z , states increase in energy with decreasing |m|, with the restriction that |m| ≤ L. Once the fullertube reaches C 90 , the minima are far enough apart that for all these eigenstates the wavefunctions are not interacting and all the states appear as (L, m, n z ) and (L, m, n z + 1) degenerate pairs.
The lifting of degeneracy of states with differing |m| for a fixed L indicates the anharmonic nature of the potential.Across the set of fullerenes, this anharmonicity is negative, as seen in H 2 @C 70 7 due to the higher degree nature of the nearly radially symmetric potential.Along the z direction, while for X@C 70 this also shows negative anharmonicity, the longer fullerenes show a slightly different pattern.Because of the larger barrier height, these states are more influenced by the double well, and so in X@C 80 , the negative anharmonicity only starts from the (0, 0, 2) state, with an unusually large energy difference between the (0, 0, 1) and (0, 0, 2) states seen in Fig 13, which could be rationalised as the endohedral atom almost "breaking-free" of the double well influence.
For X@C 90 and X@C 100 , after accounting for the (0, 0, n z ) and (0, 0, n z + 1) states arising in degenerate pairs for n z even, the states seem to show a small positive anharmonicity, below an energy of 250cm −1 .This is because as the states lie below the barrier height, the potential almost mimics what would be seen for a diatomic stretch.In that limit, the diatomic dissociates forcing positive anharmonicity on the potential.But in this case, this is the case only until the atom is in a high enough energy state as to not feel the effect of this double well, and then the negative anharmonicity returns.
Traversing the D 5d/h series, the transition energies between corresponding states decreases.This can be rationalised by considering Figure 1.As the fullerene is elongated, the double wells become deeper hence the particle is more attracted in the local minima.Spectroscopically, the transition energies between states can be measured.If this is conducted at a low enough temperature, only the ground state will be sufficiently populated for transitions to be visible, and so the fundamental transitions (0, 0, 0) → (0, 0, 1) and (0, 0, 0) → (1, 1, 0) as seen in Figure 14 are of importance.These correspond to the fundamental frequency of the 1D oscillator in the z direction and 2D isotropic oscillator in the xy plane for X@C 70 respectively.
The xy fundamental transition shows a stronger dependence on ε, as this is directly related to the "stiffness" of the isotropic oscillator and a larger stiffness implies a larger frequency.The dependence on σ is less intuitive, but as this is the length scale of the X -C interaction, the larger the value, the "stronger" the interaction leading to a larger frequency.As can be see in Figure 13, the number of (0, 0, n z ) states lying below this first degenerate state differs depending on the position of the minima and barrier height which can also suggest a larger frequency.
The fundamental transition, which in this case is always the transition in the z direction, shows very different characteristics to the xy fundamental.The dependence on ε is very weak as this modulates the depth of the double wells, seen in Figure 5a.There is a stronger dependence on σ , which characterises the distance between the minima as seen in Figure 5b and as σ increases, the minima get closer together increasing the interaction of the wavefunctions in both wells leading to a larger tunnelling splitting.The red contour indicates where the barrier height is equal to the zero point energy, whereas the yel- low contour shows a frequency of 1cm −1 , which is near the resolution of what could be seen experimentally.Spectroscopically assigning the transitions, the LJ parameters for the X -C interaction can be determined.
To compare the varying nature of these fundamental transitions across the LJ parameter space to the D 5d/h fullertube series, a quantity to consider would be the ratio of the fundamentals.However, due to the double well and small tunnelling splitting, the ratio of xy fundamental and z fundamental is a not a well behaved quantity as the z fundamental approaches 0. Instead, the z fundamental could be replaced by the (0, 0, 1) → (0, 0, 2) transition, the z hot band.In Fig fig : Fundamental Transitions, this is the transition what looks as if its the fundamental, and is better behaved.For the cases where ZPE > B, as for X@C 70 and X@C 80 , this ratio is approximately 3. For the longer fullerenes where B > ZPE, this ratio decreases to approximately 1.5 Due to the highly regular nature of the wavefunctions of these lowest 50 eigenstates, there is no discernible coupling between the z and xy modes for the low lying energy states, which was also observed in Ne@C 70 21 and H 2 @C 70 7 .However, a weak coupling between higher energy states of the appropriate symmetry may become apparent, as the energy gap between them decreases and the anharmonic nature of the potential allows them to mix, which has also been seen in Ne@C 70 21 .

IV. CONCLUSION
We have explored the parameter space of an endohedral atom trapped inside C 70 , interacting via a summed Lennard-Jones potential with the Hamiltonian as given in Equation ( 1) and the mass fixed to be the two-particle reduced mass of Ne@C 70 .By varying the two LJ parameters, a wide range of symmetric double well potentials were produced, with a variety of different barrier heights B, and positions of minima z min .These quantities may also be accessible experimentally, with the latter being seen in diffraction or NMR experiments which could also give information about the localisation of the endohedral atom in the centre or at the minima.The barrier height could possibly be calculated from the tunnelling splitting, which if on the order of 1cm −1 or greater can be visible spectroscopically.
The Hamiltonians constructed from the differing (ε, σ ) values were diagonalised using a non-orthogonal basis set, built from 1D harmonic oscillator wavefunctions expanded from the minima.Attempting to use PODVR basis sets built from 1D HO wavefunctions failed when B grew too large as minimising α in Equation (4) leads to a very small effective harmonic potential as the double well weighs down the parabola.Sinc functions also had problems converging to a tight upper bound, especially in the region where B > ZPE, seen in Fig 4 .The HO DVR functions required fixing the interval of sampling points and had an oscillatory behaviour in the convergence.The non-orthogonal basis set seemed better behaved, requiring a similar number of quadrature points to converge as the HO DVR functions and was the basis set of choice.Its possible linear dependence was alleviated by canonically orthogonalising the functions before solving the eigenvalue problem in the subspace.This non-orthogonal basis set converged the ground state energy to milli-wavenumber accuracy.As the non-orthogonal basis set performs similarly to the HO DVR functions when there is a single trapped species, it may be better suited when considering multiple encapsulated species for use in a configuration basis 45 .
The results of the diagonalisation of the three-dimensional Hamiltonian for a variety of (ε, σ ) partitioned the space into two regions in two separate ways.One is where the partition lies along the B = ZPE contour, splitting the LJ parameter space into regions where B > ZPE or B < ZPE.As the ratio of ZPE B decreases, the relevance of the non-orthogonal basis set increases as the endohedral species sits deeper in the potential wells.The other partition is where the ground state wavefunction has a single peak, or two maxima where the changeover occurs at κ z = 2.2.These two differing partitions do not oc-cur in the same place, as there is a smooth deformation of the wavefunction from a single peak which is both squashed and stretched apart to show two maxima well before the ZPE drops down to the value of B.
The wavefunctions were classified based on their statistical moments, namely the standard deviation and kurtosis.The standard deviation gives a measure of the amount of wavefunction density away from the mean, in this case the origin but it cannot distinguish the peakedness of the wavefunction.The kurtosis on the other hand is a good measure of the tailedness of the distribution, describing the density around ±σ z .When the kurtosis drops below 2.2, this is a good indication of when the wavefunction moves away from have a single peak at the origin and starts exhibiting two maxima.These statistics are only part of the covariance and co-kurtosis tensors, which are linked to the quadrupole and hexadecapole moment of the system, which perhaps could be accessible experimentally.The shape of the wavefunction could also be probed using inelastic neutron scattering, which could allow for information of its moments to be calculated.
Excited states of X@C n for fullerenes in the D 5d/h series from C 70 to C 100 were also calculated, alongside the frequencies of the fundamental transitions of X@C 70 .These transitions could be seen spectroscopically and in combination with other diffraction and NMR data, the potential of the endohedral atom could be reconstructed and the best LJ parameters of the X -C interaction can be found.The excited z and xy modes tend to show negative anharmonicity, due to the high degree nature of the potential.Eigenstates lying below the barrier height show a small positive anharmonicity, which persists until the barrier height is exceeded and the negative anharmonicity reappears.
While the parameter space explored was small, the space can be expanded straightforwardly.In the region where ZPE < B, extending the LJ parameter space can introduce a radial double well and we envisage extending our methodology into this space.We aim to use a more accurate potential energy surface, based on some electronic structure calculations which can accurately describe the non covalent interactions 48 in order to converge the eigenstates to a higher accuracy, and compare how well this model LJ potential describes the X -C interaction.Varying the LJ parameters was the alternative to changing the encapsulating cage and traversing the D 5h/d series 26 , as B and z min were more sensitive to this.We intend to compare this non-orthogonal basis set to 1D HO functions when describing the centre of mass of a small molecular species e.g.H 2 by tracking its eigenstates 7 and ground state wavefunction behaviour along the D 5h/d fullerene/fullertube series.

FIG. 2 :
FIG. 2: Orientation of the fixed C 70 cage in the (a) yz and (b) xy planes.The red, green and blue arrows correspond to the positive x, y, and z Cartesian directions respectively.

FIG. 3 :
FIG. 3: A generic symmetric double well potential with barrier height B and z min indicated with red arrows.The dashed blue lines indicate two different regions of ZPE, smaller or larger than B.

FIG. 4 :
FIG. 4: Convergence of the ground state energy and energy gap to the first excited state of a 1D Hamiltonian as in Equation (1) with (ε, σ ) Lennard-Jones parameters in (a) (43.79cm −1 , 3.03Å) and (b) (43.79cm −1 , 2.88Å).Symmetrised double minimum basis energies are shown as black crosses, HO DVR energies are red stars, and sinc DVR energies are green pentagons.Open shapes correspond to the ground state energy (left), and filled in shapes correspond to the energy gap to the first excited state (right).

FIG. 5 :
FIG. 5: Features of the double well, (a) barrier height and (b)position of minima, in the unique direction as defined in Figure1of the PES tracked across the LJ parameter space.

1 FIG. 6 :
FIG.6: 1D slices of Lennard-Jones potential along the z and x Cartesian directions.Solid line indicates V (0, 0, z) whereas dashed is V (x, 0, 0).V (0, y, 0) is not shown as for this region it is almost identical to V (x, 0, 0), differing only when the endohedral species is closer to the cage edge.The red curves have LJ parameters (ε, σ ) as (60, 2.98) and the blue curves have (140, 2.88).

FIG. 7 :
FIG.7: Ground state energy of X@C 70 with Hamiltonian shown in(1) with M corresponding to20 Ne in wavenumber across LJ parameter space.The red contour indicates the path where the barrier height B and zero-point energy ZPE are equal.

FIG. 8 :
FIG. 8: Wavefunction moments, (a) standard deviation (b)kurtosis and (c) prolateness, of the ground state eigenfunction calculated along the unique z axis across the LJ parameter space.

FIG. 9 :
FIG. 9: Ground state wavefunction prolateness and kurtosis in (b) tracked along the contours in the Lennard-Jones parameter space in (a).The black line is at fixed σ = 2.93Å, the blue line is at fixed ε = 100cm −1 the green line is along the path σ = 3.08 − 0.0015ε and the red curve is along the ZPE = B contour shown in Fig 7.

1 FIG. 13 :
FIG. 13: Energy levels of lowest 50 eigenstates of X@C n in cm −1 , with fixed LJ parameters (ε, σ ) as (43.79cm −1 , 3.03Å).C 70 is in black, C 80 in red, C 90 in blue, and C 100 in green.States are indexed with quantum numbers (L, m, n z ).m = 0 shown with solid lines, |m| = 1 as dotted lines, |m| = 2 as dashed lines, and |m| = 3 as dot-dashed lines.States in C 90 and C 100 appear as near-degenerate pairs with quantum numbers (L, m, n z ) and (L, m, n z ) + 1 for n z even.
States with |m| > 0 are doubly degenerate with quantum numbers ±m.States with m = 0 are solid lines, |m| = 1 are dotted lines, |m| = 2 are dashed lines and |m| = 3 are dot-dashed lines.

1 FIG. 14 :
FIG. 14: Transition energies of the fundamental transition in the (a) anistropic z direction and (b) in the 2D isotropic xy plane.

TABLE I :
Lennard Jones parameters for Ne in a fullerene.