Experimental, theoretical and computational investigation of the inelastic neutron scattering spectrum of a homonuclear diatomic molecule in a nearly spherical trap: H₂@C₆₀

Abstract

In this paper we report a methodology for calculating the inelastic neutron scattering spectrum of homonuclear diatomic molecules confined within nano-cavities of spherical symmetry. The method is based on the expansion of the confining potential into multipoles of the coupled rotational and translational angular variables. The Hamiltonian and the INS transition probabilities are evaluated analytically. The method affords a fast and computationally inexpensive way to simulate the inelastic neutron scattering spectrum of molecular hydrogen confined in fullerene cages. The potential energy surface is effectively parametrized in terms of few physical parameters comprising an harmonic term, anharmonic corrections and translation–rotation couplings. The parameters are refined by matching the simulations against the experiments and the excitation modes are identified for transfer energies up to 215 meV.

---

1 Introduction

The topic of confined particles lies at the foundations of quantum mechanics and textbook examples of a particle confined in a well or moving in a harmonic trap provides the starting point for a deeper understanding of many quantum systems. With the advent of “molecular surgery” techniques in organic synthesis, pioneered in the laboratories of Komatsu,¹ Murata² and more recently Whitby,³ molecular complexes have become available in which a variety of quantum rotors may be fully enclosed and confined inside a fullerene cage. These are the small molecule endofullerenes, of which H₂@C₆₀ is the primary example, representing a model system of its kind.⁴

With a hydrogen molecule entrapped inside C₆₀, the endofullerene H₂@C₆₀ affords an actual realization of a quantum rotor confined in a three dimensional cage. The dynamics of H₂ encapsulated within a fullerene is fundamentally quantum mechanical for several reasons. Firstly, owing to its low mass, free molecular hydrogen is characterized by a large rotational constant with sparse levels whose separation is in excess of 170 K in temperature units. Secondly once hydrogen is confined inside the C₆₀ nano-cavity its translational degrees of freedom becomes quantised. Significantly, since the internuclear bond length in hydrogen (0.74 Å) is in the same scale as the fullerene radius, the quantised rotational and translational degrees of freedom become coupled. Thirdly, since the nuclei comprising H₂ are indistinguishable fermions, the antisymmetry principle applies in determining the allowable eigenstates and there exists ortho- and para-nuclear spin isomers of H₂. At low temperature these remain largely decoupled from each other and the two co-exist as independent and identifiable species in the absence of time-dependent inhomogeneous magnetic fields.

The study of the quantum dynamics of H₂@C₆₀ is relevant for understanding the interaction between carbon based materials and hydrogen. Molecular endofullerenes are ideal compounds in which the guest molecules are relatively isolated from each other within the homogenous environment provided by the cages. As a result, the spectrum of translation–rotation energy levels is sparse and when probed by various spectroscopic techniques, narrow well-separated peaks are observed.^5–8 This reveals in fine detail the quantum dynamics of the entrapped translator–rotator and enables the potential energy surface (PES) experienced by the rotor to be probed meticulously along with the various interactions that characterise the Hamiltonian. With its unique ability to induce transitions between ortho- and para-spin-isomers, inelastic neutron scattering (INS) is a prime technique for studying confined hydrogen.⁹

Pioneering work by the group of Bačić using a rigorous computational approach to rigorously solve the Hamiltonian^10–14 has been instrumental in developing an understanding of confined H₂, including endofullerenes.^15–17 Their methodology has been extended to determine the INS spectrum of hydrogen confined in nanocavities^18,19 and in fullerenes.^20,21 There, the potential energy surface is modelled using modified Lennard-Jones functions to represent the atom–atom potentials, then the Hamiltonian is solved exactly from which the spectrum may be determined.

The scope of this paper is different. Starting from the experimental spectrum, the observed low energy transitions are described in terms of excitations of a rotator angularly coupled to a harmonic oscillator. The quantum dynamics of a homonuclear diatomic molecules in a weakly perturbed isotropic harmonic trap is investigated by expanding the confining potential in symmetry adapted multipoles of the coupled angular variables. The Hamiltonian and the INS transition probabilities are represented analytically in the basis of the coupled rotor and harmonic oscillator. The INS spectrum is simulated by refining few Hamiltonian parameters, including harmonic and anharmonic terms along with translation–rotation coupling constants. In this methodology, the PES is not built from specific microscopic interactions but it is determined from a direct fit of the Hamiltonian to the experimental spectrum with the advantage that the fitted parameters represent physical attributes, such as the harmonic potential, the anharmonicity and translation–rotation coupling constants.

In the context of the spectroscopy of H₂@C₆₀, it is well-established that the quantum states of the coupled harmonic translator–rotator can be described in terms of the quantum numbers (n,l,j,λ).²² The entrapped rotor possesses both rotational and translational energy. The rotational states are characterised by j, while the translational states, defined by displacements of the centre-of-mass of the particle within the 3D-cage, are characterised by n and l. The former, n, is the principle quantum number while the latter, l, determines the translational angular momentum. Rotational and translational angular momenta couple according to the usual vector addition rules and we identify the quantum number λ to characterise the coupled translation–rotation states. A more complete account of the rotational and translational states will be discussed in Section 4. Finally, applying the antisymmetry principle, the Pauli Exclusion Principle (PEP) demands that the nuclear spin isomer para-H₂ is a nuclear spin singlet with total spin I_N = 0 and characterised by even j, while ortho-H₂ is a nuclear spin triplet with I_N = 1 and odd values of j. In this way the spatial and spin degrees of freedom are entangled. This has consequences for experiments since ortho ↔ para transitions are spin-restricted and forbidden to photon spectroscopy, whereas neutron scattering interactions couple space and spin enabling ortho ↔ para transitions to be directly observed in the INS spectrum with high intensity.

The paper is organised as follows. First the experimental INS spectrum is presented in Sections 2 and 3. The analytical approach used to represent the Hamiltonian and to determine the INS transition probabilities is discussed in Sections 4 and 5, respectively. The method is then applied to simulate the INS spectrum of H₂@C₆₀ in Section 6. The parameters characterising the cage potential and the spectral assignments are discussed in Section 7.

2 Experimental details

INS experiments were conducted using the IN1-Lagrange spectrometer at the Institut Laue-Langevin (ILL) in Grenoble.²³ This instrument is supplied from the “hot source” moderator of the ILL reactor, resulting in a high-flux beam of energetic neutrons. A choice of three different single crystal monochromators is available in a primary spectrometer, where a Bragg angle defines the incident energy of the monochromatic neutron beam arriving at the sample. The monochromator crystals of Si or Cu are mounted in a double focusing geometrical arrangement, enabling selection of the Si(311), Cu(220) or Cu(331) Bragg reflections. By combining the 3 monochromators, the incident neutron energy may be scanned in the range from 5 meV to 1 eV. The scattered neutrons enter a secondary spectrometer comprising a large area array of pyrolytic graphite analyzer crystals. The Bragg scattering geometry in the secondary spectrometer ensures only neutrons of a fixed energy (4.5 meV) are detected by a ³He detector, with neutrons being collected from a wide solid angle subtended at the sample. By scanning the incident neutron energy, INS spectra were recorded in the energy transfer range 12.7 ≤ ΔE ≤ 215 meV. Over this range of energy transfer, the momentum transfer of detected neutrons scattered by the sample varies systematically in the range 2 ≤ q ≤ 12 Å. For the Si monochromator the resolution is given by the resolution of the analyser, which is 0.8 meV. For the Cu monochromators the resolution function of the instrument is a constant percentage of the incident neutron energy. Identifying singlet transitions in the INS spectrum of H₂@C₆₀, the Gaussian resolution function was measured as follows (FWHM; full-width half maximum): Cu(220) 2.3%; Cu(331) 1.5%.

The sample of H₂@C₆₀ was prepared according to published methods.¹ In a final stage of purification, the presence of any intercalated solvent was substantially reduced by sublimation. The powdered sample with mass 101 mg was uniformly loaded inside a rectangular Al foil sachet, folded into a cylindrical annulus of approximate diameter 1 cm and mounted in the Lagrange spectrometer. The sample temperature was maintained at 2.5 K using a closed cycle He refrigerator. To subtract background and scattering from Al and from the C₆₀ cage, a “blank” sample of empty C₆₀ cages with matching mass was prepared in an identical Al foil sachet.

INS spectra of H₂@C₆₀ were recorded using the three monochromators and the corresponding energy transfer ranges were selected to provide the best available resolution as follows: Si(331) 12.7 < ΔE ≤ 27 meV; Cu(220) 27 < ΔE ≤ 63 meV and Cu(311) 63 < ΔE ≤ 215 meV. The intensities of the INS peaks recorded with the three monochromators were normalized to the same scale by comparing the integrated intensities of peaks recorded in overlapping energy transfer regions. The systematic error in energy transfer is estimated to be 0.1 meV below approximately 50 meV rising to 0.5 meV at the highest energy transfers recorded. Random errors in energy transfer were small by comparison, typically 0.02 meV for singlet peaks. Unless otherwise stated, these error bars will apply in the following discussions.

3 Experimental results

The INS spectrum of H₂@C₆₀ is presented in Fig. 1. A mass-matched “blank” has been subtracted, ensuring the spectrum arises only from scattering by H₂ inside the fullerene cage. Compared with previous studies using time-of-flight spectrometers, this new spectrum substantially extends the range of energy transfer recorded from 60 meV to beyond 250 meV. Furthermore, in the range 30–60 meV the resolution function is typically two times better than the previously published IN4 spectrum (λ_n = 1.2 Å),²⁴ enabling TR fine structure to be resolved for the first time in the high energy modes. The sample is at low temperature so that only the ground para-H₂ and ground ortho-H₂ states are occupied. At these low temperatures and in the absence of suitable time-dependent magnetic interactions, ortho–para conversion is extremely slow with no conversion being observed on the timescale of the experiments (several days). Therefore the two nuclear spin-isomers exist as distinct species and all INS transitions originate in the ground states, either (0,0,0,0) or (0,0,1,1), using the notation (n,l,j,λ) to identify the energy levels.

Fig. 1 The INS spectrum of H₂@C₆₀ at 2.5 K as recorded on the IN1-Lagrange spectrometer at ILL, Grenoble (France). The vertical axis represents neutron counts in arbitrary units.

At low temperature, and particularly in the lower energy transfer region, the spectrum is sparse with sharply resolved peaks revealing a remarkable quantity of detail. H₂@C₆₀ is clearly an extremely well-characterized and homogenous system, revealing the quantum nature of the H₂ molecule and its wavelike nature in exquisite detail.

Fig. 2 shows an extract from the experimental spectrum for transfer energy up to 73 meV. The first rotational transition (0,0,0,0) → (0,0,1,1) is clearly observed as the singlet peak with energy transfer 2B = 14.7 meV, interconverting the ground states of para-H₂ and ortho-H₂. Another singlet peak is the second rotational transition (0,0,1,1) → (0,0,2,2) which has energy transfer approximating 4B and is centered on 29.2 meV. A third singlet with width equal to the resolution function appears at 72.2 meV. This energy approximates to 10B so the peak assigned to the ortho–ortho pure rotational transition (0,0,1,1) → (0,0,3,3). These observations indicates that barriers impeding free rotations are small in H₂@C₆₀.

Fig. 2 Expansion of the INS spectrum of H₂@C₆₀ at 2.5 K of Fig. 1 in the range 10–75 meV. Each peak is assigned to an INS transitions for which the labels indicate the value of n, l, j in the final states. Transitions originating from the para-H₂ ground state are in blue and transitions originating from the ortho-H₂ ground state are in red, respectively. The vertical axis represents neutron counts in arbitrary units.

The prominent band centered on 22.5 meV received particular interest in our earlier papers,²⁴ since it represents the transition from the ground state of ortho-H₂ to the first excited translational state of ortho-H₂, (0,0,1,1) → (1,1,1,λ). This is a ubiquitous “particle-in-a-box” mode and in the Lagrange spectrum this is observed with similar resolution to the earlier IN4 spectrum (λ_n = 1.6 Å). This mode is a partially resolved triplet comprising the λ = 0, 1, 2 components of the excited state. This band was significant since it represented the first observation in an INS spectrum of a splitting arising from TR coupling of rotational and translational angular momentum. Moving to higher energy transfer, a doublet is observed with energies 37.4 and 38.5 meV. Considering the energy level diagram and consistent with the assignments in our earlier papers,^8,24 we infer this band is the transition (0,0,0,0) → (1,1,1,λ) involving the same final state triplet referred to in the previous paragraph, but this time originating in the ground para-H₂ state. Within the resolution of the spectrometer only two components of this band can be identified while the third one is apparently missing. We shall return to this point more widely in Section 7 in the context of the discussion of the INS transition probability for diatomic molecules in spherical potentials.

The second excited translational state of ortho-H₂, n = 2, j = 1, is revealed in a doublet with energy transfers 45.8 and 47.0 meV and in the small peak at 49.8 meV. These peaks were observed with the Cu(220) monochromator and have the resolution line width. The former doublet has a partner observed at 61.3 meV while the latter peak is closely related to the singlet at 64.2 meV observed with the Cu(331) monochromator, which also has the resolution line width. Allowing for systematic errors that arise from two different monochromators, these two set of peaks differ in energy transfer by 2B, so we are able to make the assignment to the same final state within the manifold n = 2, j = 1; the lower energy of the two originates in ortho-H₂, while the higher energy peak originates in para-H₂.

The peak in the region 51.5 meV has been observed with the Cu(220) monochromator. It is slightly broader than the resolution line width and slightly asymmetric so we can infer the feature has multiple components. As with the examples above, applying an offset equal to a multiple of 2B provides an energy that coincides with a known splitting in the spectrum. In this case subtracting 4B coincides with the translational splitting so we can infer these two peaks correspond to the ortho–para transition from (0,0,1,1) → (1,1,2,λ).

Finally, two low intensity peaks occur at 32.3 and 34.3 meV. These identify the states n = 2, j = 0, corresponding to the transitions (0,0,1,1) → (2,0,0,0) and (0,0,1,1) → (2,2,0,2).

As the above few examples show, with good instrument resolution and a spectrum which is sufficiently sparse, substantial progress can be made in the assignment of n, l, j quantum numbers based solely on the experimental data. This enables the energy level diagram to be determined for all low lying states up to approximately 73 meV. As both para to ortho and ortho to ortho transitions are allowed, and the same final states are represented by multiple peaks, the assignments may be confirmed by correlating transition energies. However, with increasing energy transfer, the number of transitions significantly increases while anharmonicities lead to non-integer relationships between peak energies. It becomes clear that a definitive assignment of the INS peaks, including the TR multiplet components, increasingly demands the assistance of a more sophisticated theoretical framework. Therefore we shall return to consider more detailed assignments in the spectrum in Section 7 after describing the computational simulation of the H₂@C₆₀ spectrum.

4 Quantum dynamics

In this section we discuss the quantum dynamics of a H₂ molecule confined in a C₆₀ cage. The C₆₀ cage is assumed to be rigid, fixed and non-rotating. The laboratory reference frame is taken with origin in the centre of the cage. The geometric configuration of the cage is completely identified by the orientation of a frame fixed with the cage C with respect to the laboratory frame L. Mathematically the orientation of the cage frame with respect to the laboratory frame is described by the set of three Euler angles Ω_LC. For example, in the case of C₆₀, the cage fixed frame can be chosen so that the Z axis is directed along a C₅ axis with the X axis along one of the 5 orthogonal C₂ axes and the Y axis determined by the right-handed condition. The configuration of the diatomic molecule is determined by the centre of mass vector R and the internuclear vector r.

The molecule–cage interaction is described by the potential V(R,r;Ω_LC). The function V can be thought as the PES of an H₂ molecule confined inside the fullerene cage. It is worth noting that any symmetry operation Ô of the icosahedral point group I_h acting on all the position coordinates leave the potential unaffected:

V(R,r;Ω_LC) = V(ÔR,Ôr;Ô[Ω_LC]) = V(ÔR,Ôr;Ω_LC) (1)

and the last equality follows from the fact that the Euler angles are unchanged under Ô. As a consequence the potential is invariant under operation of the symmetry group acting only on the variables R and r, see also ref. 25 and 26. In the following the dependence of V on Ω_LC will be understood. To discuss the effect of symmetry on the potential at fixed cage orientation, it is convenient to expand out the angular dependence in multipoleswhere the bipolar spherical harmonics are a complete set of functions for the angular variables defined by²⁷and where C is a Clebsch–Gordan coefficient and {θ_R,ϕ_R} and {θ_r,ϕ_r} are polar angles of the centre of mass and internuclear vector, respectively. Λ is the rank and M_Λ is the projection of the bipolar harmonic. For sake of simplicity, the dependence on the ordered set of angular variables will be understood henceforth.

The symmetry of the confinement is reflected in the number and types of multipoles present in the potential expansion, eqn (2). In icosahedral symmetry it is possible to prove that the rank of the multipoles that form invariant functions is restricted to Λ = 0, 6, 10, 12, 16.²⁸ Since the centre of the potential is a point of inversion symmetry, only multipoles with L + J = even are allowed. For multipoles with rank Λ = 0, only terms with Λ = 0, M_Λ = 0 and L = J are allowed. The next lowest non-zero rank terms for an icosahedral potential is:²⁹

In treating the quantum dynamics of H₂@C₆₀, the effect of multipoles with Λ > 0 can be neglected in first instance and eventually reconsidered perturbatively.²⁵ Retaining only terms with Λ = 0 is equivalent to treating C₆₀ as spherical. The internal symmetry of the homonuclear molecule determines the invariance of the potential under the exchange of r in −r which in turn implies that only even-J bipolar harmonics are present in eqn (2).

The vibrationally averaged effective five-dimensional Hamiltonian for a diatomic molecules moving in a spherical potential can be written as

where P̂ is the linear momentum operator associated with the centre of mass vector, Ĵ is the angular momentum operator associated with end to end rotations, B(ν) is the effective rotational constant in the vibrational state ν and ^νV^LL;k₀₀ are coefficients of the expansion of the averaged molecule–cage interaction potential in bipolar harmonics of the angular variables. For a homonuclear molecule, only even k terms enter in the expansion of the potential in eqn (5) because of the symmetry of the potential under inversion. In this work the attention is restricted to H₂@C₆₀ in the ground vibrational state ν = 0 and the label ν will be dropped henceforth.

It is a general statement of quantum mechanics that the total angular momentum of a system with spherical symmetry commutes with the Hamiltonian, so providing “good” quantum numbers for labelling the energy levels. For a diatomic molecule in a spherical trap, the total angular momentum [capital Lambda, Greek, circumflex] = L̂ + Ĵ is given by the sum of the orbital angular momentum of the centre of mass L̂ and of the end to end rotational angular momentum Ĵ. The bipolar spherical harmonics are eigenfunctions of the angular momentum operators [capital Lambda, Greek, circumflex]² and [capital Lambda, Greek, circumflex]_z with eigenvalues λ(λ + 1) and m_λ, respectively. They provide a complete basis set for the coupled angular motion.

A complete basis set for the radial motion is provided by the set of radial eigenfunctions for the Schrödinger equation of a particle in a central potential. Early ab initio determinations of the potential energy surface of H₂@C₆₀ suggest to use the eigenfunctions of a three-dimensional isotropic harmonic trap. These are discussed in the appendix, Section A.1. The coupled harmonic translation–rotation wave-function is written in ket notation as

using the notation of Section A.1 and the definition of bipolar harmonic of eqn (3). Here is the scale length associated to the harmonic oscillator with mass M and frequency ω. The ket in eqn (6) is an eigenstate of the spherical three-dimensional oscillator represented by the Hamiltonian in eqn (5) when terms with L = 0 and up to k = 2 are considered. In this case . The zero-order energy of the undecoupled translator and rotator is

E⁰_n,j = ħω(n + 3/2) + Bj(j + 1) (7)

Anharmonic terms with L = 0, k > 2 mix states with different n while translation–rotation couplings are characterized by L = J > 0 and mix states with different n, l, j. In presence of anharmonicities and translation–rotation couplings the eigenfunctions of the Hamiltonian given in eqn (5) are characterized by definite angular momentum quantum number λ and m_λ and can be written aswhere the coefficients c do not depend on m_λ. In the following we will use the notation E_λ to denote the energy level comprising the set of 2λ + 1 degenerate states of eqn (8). The symmetry under inversion implies that parity is a good quantum number. Since the parity under inversion for the bipolar harmonics states is (−1)^l+j, it follows that either even- or odd-sum l, j pairs appear in the expansion (8).

5 INS transition probability

In a typical neutron scattering “spectrum” the intensity of a transition between an initial energy level i and a final energy level f is determined by the double differential cross section:³⁰where P_i is the population of the energy level i, represents the transition probability for the scattering of neutrons and the delta function δ enforces the energy conservation in the scattering process. As shown with greater detail in Section A.4, the INS transition probability for homonuclear molecules in unpolarised experiments can be written as a product of a spin s_fi and a space part S_fi:Specifically in the case of molecular hydrogen, I = 1/2 and I_N = 0 for para-hydrogen and I_N = 1 for ortho-hydrogen, so that the spin dependent part of the INS cross section becomeswhere o and p stands for ortho and para, respectively and where b_c and b_i are the bound coherent and incoherent scattering lengths of bound ¹H, respectively. The space dependent part S_fi for transitions between an initial energy level i = E_λ and a final energy level is obtained by tilting the reduced matrix element of the INS transition probability from the spherical basis to the appropriate eigen system basis:where the INS reduced matrix element is given by:

Interestingly, the derivation of the transition probability via the Wigner–Eckart theorem, as discussed in Section A.4, allows us to recover the selection rule for INS transitions in diatomic molecules confined in potentials of spherical symmetry, originally found by Xu et al.^20,31 and later extended and explained by Poirier in the context of group theory.²⁶ Indeed in the evaluation of the reduced matrix element of eqn (13) the terms with

are exactly 0 since both the Clebsch–Gordan term C^γ0_α0β0 and the reduced matrix elements of bipolar harmonics involve 3-j symbols with 0 projections, see eqn (29) and (30). These 3-j symbols are zero if the sum of the arguments are odd, see eqn (20). When the initial level is characterized by λ = 0 total angular momentum, l = j and γ = λ′ because of the triangular condition for the rows and column of the 9-j symbols implied by the reduced matrix element of the bipolar harmonic, eqn (29). By subtracting the last two identities out of the first one in 14, it follows that the reduced matrix element 〈n′,l′,j′,λ′‖h^±_γ‖n,l,j,λ〉 is zero when the initial level has λ = 0 and the final level angular numbers satisfy λ′ − l′ − j′ = odd. The selection rule for the transition probability follows from the observation that l′ + j′ is either even or odd for all the coefficients in eqn (12) as a consequence of the inversion symmetry for diatomic molecules in a spherical trap, see Section 4. The equivalent selection rule for final states with λ′ = 0 follows exchanging the role of primed and non-primed labels in the discussion above.

6 Calculations

The methodology described in Sections 4 and 5 can be applied to determine the INS spectrum of H₂@C₆₀ according to the following prescription. Firstly a finite coupled basis set is chosen by fixing the maximum number of translational and rotational excitations to n_max and j_max, respectively. The para (even j) and the ortho (odd j) coupled states |n,l,j,λ,m_λ〉 are gathered separately in blocks with the same angular momentum λ. For each spin isomer and each λ-block, the finite matrix basis representation of the Hamiltonian (5) is built in the basis |n,j,l,λ,m_λ = 0〉 using eqn (7), (22) and (23). The rotational constant in the ground vibrational state is written as B = B₀ − D₀j(j + 1) where the coefficient D₀ accounts for centrifugal corrections to the rotational energy.³² The potential expansion is truncated to include the harmonic term V^00;2₀₀, two anharmonic terms, V^00;4₀₀ and V^00;6₀₀ and one translational–rotational coupling term V^22;2₀₀. The blocks with different angular momentum λ are diagonalized separately. The eigenvalues are used to determine the thermal populations, the levels and the transition energies that fall within the observation range. Assuming no spin-isomer interconversion, the ortho to para ratio can be fixed to 3 and the thermal equilibrium populations can be determined for the two spin isomer independently. The eigenvector coefficients, eqn (8), are used to calculate the total intensity of the INS transitions as a function of q, by evaluating the double differential cross section, eqn (9), viaeqn (11) and (12). For numerical simulation the values b_c = −3.74 × 10⁻¹⁵ m and b_i = 25.27 × 10⁻¹⁵ m have been used as scattering lengths for the coherent and incoherent scattering length of bound ¹H, respectively.³³

To achieve a more direct visual comparison between experiments and simulations, the transition probability, eqn (10), is integrated over the q-range appropriate to the instrumental apparatus at energy ΔE_fi = E_f − E_i. In an inverse geometry neutron spectrometer the final energy of the scattered neutron is fixed to E_n, while the incident neutron energy is scanned through. A detector at angle θ with respect the direction of the incident neutron beam records scattered neutrons at transferred momentum

so that integrating over q across the polar angular opening of the instrumentThe integration can be effectively performed averaging over an appropriate angular distribution describing the geometry of the apparatus. This procedure gives rise to a stick peak for each possible transition. To emulate finite width transitions the peaks are finally convoluted with the instrumental resolution, assuming negligible intrinsic line-width in H₂@C₆₀.

In the simulation of the IN1-Lagrange spectrum of H₂@C₆₀ the number of scattered neutrons was determined by averaging over a uniform distribution of 20 polar angles in the nominal angular ranges 34° ≤ θ ≤ 70° at fixed E_n = 4.5 meV and using the instrumental resolution discussed in Section 2.

7 Discussion

The INS spectrum of H₂@C₆₀ was simulated using a set of purpose-built packages in Mathematica.³⁴ In the actual calculations, the basis was defined by n_max = 11 and j_max = 7, giving rise to 10 192 para states and 13 104 ortho states, of which only 800 distinct reduced para states and 972 reduced ortho states were necessary for the evaluation of the Hamiltonian and of the INS transition probability. The initial rotational parameters were fixed to the values of free H₂: B₀ = 7.357 meV and D₀ = 0.0056 meV.^35,36 Starting from the uncoupled oscillator with translational excitation ħω = 22 meV, the potential parameters , Ṽ^00;4₀₀ = ρ⁴V^00;4₀₀, Ṽ^00;6₀₀ = ρ⁶V^00;6₀₀ and Ṽ^22;2₀₀ = ρ²V^22;2₀₀ were varied until satisfactory agreement with the experimental spectrum was found. The optimized values are reported in Table 1.

Table 1 Summary of the optimized values of the Hamiltonian parameters used for calculating the INS spectrum of H₂@C₆₀ in Fig. 3

Hamiltonian parameter	Optimized value
B 0/meV	7.375
D 0/meV	0.012
ρ/Å	0.311
Ṽ ^00;400/meV	−3.90
Ṽ ^00;600/meV	1.52
Ṽ ^22;200/meV	8.9

---

Fig. 3 shows the comparison between the experimental and calculated spectrum of H₂@C₆₀ at 2.5 K.

Fig. 3 Comparison between the experimental and calculated spectrum of H₂@C₆₀. For clarity the plot is split into two equally wide panels spanning the transfer energy ranges from 10 meV to 215 meV. The black line with error bars represent the spectrum of H₂@C₆₀ recorded at 2.5 K on IN1-Lagrange at ILL. The bars represent the calculated stick spectrum of H₂@C₆₀ at 2.5 K, where the colors blue and red are used to distinguish transitions originating in the para and ortho ground states, respectively. The orange line overlayed to the experimental spectrum represent the calculated spectrum obtained by convolving the stick spectrum with the instrumental resolution of IN1-Lagrange. For clarity the peaks in the stick spectrum are scaled by a factor 0.5 with respect to the numerical values and vertically shifted.

The calculated spectra confirm the peak assignment up to 73 meV given in Section 3, that were mainly based on the approximation of uncoupled rotation and translation, see eqn (7). However even at low energy the presence of couplings has a major impact on the appearance of the spectrum. By providing intensities and transition energies as well, the calculations establish a firm framework for the analysis of the spectrum.

The experimental peaks at 14.7, 29.2 and 72.2 meV are singlet rotational peaks, corresponding to the para–ortho (0,0,0,0) → (0,0,1,1) and ortho–para (0,0,1,1) → (0,0,2,2) transitions, respectively. The presence of anharmonic and translation–rotation coupling has a minor impact on pure rotational transition between states with n = 0, l = 0. However other rotational transitions are more difficult to identify because of the reduced intensity and their use as milestones for spectral assignment is less poignant.

The structure underlying multi-component peaks is clearly resolved in the stick spectrum. The experimental peak at 22 meV is resolved into three main components corresponding to the intra-ortho translational transitions (0,0,1,1) → (1,1,1,λ) with λ = 1, 2, 0 in order of increasing energy. As shown in ref. 37 using group theoretical methods and successively in ref. 25 by using perturbation theory, the order of the λ levels in the (1,1,1,λ) multiplet is dictated by the positive sign of the translational–rotational coupling term Ṽ^22;2₀₀. As observed previously the λ = 0 component is much weaker in intensity than the other two and, importantly, these observations are in agreement with the simulation of the H₂@C₆₀ spectrum by Xu et al.^15,20 In a similar manner, the experimental peaks at 37.4 and 38.5 meV are assigned to the λ = 2 and λ = 0 components of the para–ortho translational transition (0,0,0,0) → (1,1,1,λ), respectively. There is no evidence of the λ = 1 component of the triplet. This highlights a striking prediction from the earlier papers,^20,21 namely the existence of a selection rule for INS transitions in diatomic molecules in a spherical confinement.^26,31

Transitions towards levels with two quanta of translational excitations have been identified in the experimental discussion. The peaks at 32.3 and 34.3 meV are assigned to the transitions (0,0,1,1) → (2,2,0,2) and (0,0,1,1) → (2,0,0,0), respectively. Transitions towards states with n = 2 and j = 1 reveal a larger degree of mixing between different translational–rotational states. The calculations identifies the transitions towards (2,2,1,2) and (2,2,1,3) with the peaks at 45.8 and 47.0 meV, respectively. Concealed under the latter peak, there is a transition to an energy level comprising a superposition of the (2,0,1,1) and (2,2,1,1) states with a ratio close to 1 : 2. Similarly the peak at 49.8 meV corresponds to a transition to an energy level comprising a superposition of the (2,0,1,1) and (2,2,1,1) states in a ratio close to 2 : 1. The partner transitions to these superposition states originating from the para ground state are found at 61.3 and 64.2 meV, respectively, while the para transition towards (2,2,1,2) is forbidden and the para transition towards (2,2,1,3) has negligible intensity.

The peak centered at 51.5 meV is resolved in three transitions starting from the ortho ground state towards (1,2,1,λ) state with λ = 2, 3, 1 in order of increasing energy.

The stick spectrum demonstrates the existence of para–para transitions that are easily overlooked in the experimental spectrum. Owing to the small ratio between the coherent and incoherent cross section for ¹H, para–para peaks generate a very low scattered neutron count and are difficult to observe. For example the fundamental translational transition in the para manifold (0,0,0,0) → (1,1,0,1) is shadowed by the ortho transitions (0,0,1,1) → (1,1,1,λ).

At higher energy transfer the absolute resolution becomes worse, the number of accessible transitions increases and the peaks become relatively less sparse and less intense. This is particularly evident in the region above 96 meV. The calculations identify 162 transitions originating from the ground para state and 191 transitions originating ground ortho state in the energy range between 96 and 215 meV. In contrast to the low energy case, where most of the levels are represented by the nominal (n,l,j,λ) state with at least 95% probability except where noted, states with different n, l, j become more and more mixed at high energies, so hampering their use as labels for the energy levels. A semi-empirical approach that tries to fit more prominent features to a minimal number of Gaussian lines may still work where few high intensity transitions are clustered.⁸ However in other regions where transitions are more spread in energy and there is no dominant component, it is practically impossible to produce any meaningful assignment of this type. Simulations remains the only method to make sense of the complexity of the spectrum at high energy transfer. Indeed the simulated spectrum provide a realistic representation of the experimental observations over all the accessed energy range. The simulated transitions follow closely but not exactly the experimental ones. Higher order spherical terms may improve the match at the expense of increasing the parameter space. However improved accuracy appears to be not worthwhile with the current experimental resolution also considering that shifts and splittings of the m_λ components in the energy levels of order of fractions of meV are introduced by terms of icosahedral symmetry³⁸ as well as intermolecular effects.³⁹

8 Conclusions

In this paper the low-temperature INS spectrum of H₂@C₆₀ covering excitations up to 215 meV was presented. The assignment of the experimental peaks was aided by developing a model which makes use of the nearly spherical symmetry of the confinement. The guest–host interaction was expanded in terms of coupled angular spherical bipolar harmonic functions. The effective potential was then described by few physical parameters comprising a harmonic term and terms describing anharmonic corrections and translational–rotational couplings. This approach affords a good flexibility in representing the overall non-bonding interaction between the hydrogen molecule and the fullerene cage when compared with other interaction models based on interactions between atomic units. The effective potential then represents a benchmark for any ab initio determination of the PES in H₂@C₆₀. The model does not involve any numerical integration: matrix elements and transition probabilities are shown to integrate into closed algebraic expressions in the coupled rotational translational oscillator basis. The analytical expressions may be useful in the analysis of the q-dependence of the INS transitions probability in similar systems. More importantly, the numerical evaluation of these expressions comes at a relatively inexpensive computational cost. Starting from the uncoupled harmonic system, the effect of each term on the spectrum can be varied independently until the spectrum is reproduced correctly. The agreement between experiments and calculations is overall very good. The INS transitions are well reproduced indicating that the icosahedral non spherical terms give rise to minor splittings and shifts within the instrumental resolution.

The theoretical and computational model developed in this paper may be useful for investigating the quantum dynamics and the INS spectrum of confined molecular hydrogen in systems such as clathrates hydrates,^40,41 zeolites^42,43 and MOFs.^9,44,45 We anticipate its extension to treat heteronuclear diatomic molecules as well as water molecules confined in C₆₀. In systems with high symmetry or when the effect of non spherical multipoles are small, such as H₂@C₆₀, realistic INS spectra can be generated promptly. For systems with lower symmetry the method can be used to provide a starting approximate description useful to identify the main features of the INS spectrum and it can be eventually augmented to include non-spherical terms perturbatively.

Appendix A: mathematical appendix

In this section we provide an extended presentation of the mathematical aspects that lie at the basis of the algebraic and numeric treatment developed for the analysis of the INS spectroscopy of H₂@C₆₀. The solution of the Schrödinger equation for the spherical oscillator provides the radial wave-functions for the coupled translational–rotational states. The matrix elements for local operators appearing in the Hamiltonian and in the INS probability are separated into a radial and an angular part. The radial matrix elements are integrated into closed algebraic forms. The angular matrix elements are evaluated using the Wigner–Eckart theorem. Finally the INS transition probability for diatomic molecules in spherically symmetric confinement is derived.

A.1 The isotropic three-dimensional harmonic oscillator

Following the treatment presented in ref. 46 and the relations between the hypergeometric confluent function ₁F₁ and Laguerre polynomials,⁴⁷ the normalized eigenfunctions of the Schrödinger equation for the three-dimensional spherical oscillator with mass M and angular frequency ω are written as the product of a radial and an angular term:in which is the scale length associated to the harmonic oscillator, L^(α)_n is a generalized Laguerre polynomial⁴⁷ andis a normalization constant. The principal quantum number n is a non-negative integer that determines the energy, E_n = ħω(n + 3/2), of the eigenfunction in eqn (17). The orbital angular momentum quantum number l is a non-negative integer. For a given n the angular momentum quantum number is limited to the range l = n, n − 2,…, 3, 1 for odd n and l = n, n − 2,…, 2, 0 for even n. The degeneracy of the level with energy E_n is g_n = n(n + 1)/2.

A.2 Clebsch–Gordan coefficients and 3-j symbols

Clebsch–Gordan coefficients appear naturally in quantum mechanics when adding the angular momentum of two subsystems as in formula (3). 3-j symbols are a more symmetrized version of the Clebsch–Gordan coefficients. The two symbols are related byThe 3-j symbols are invariant under even transpositions of the columns and take the sign (−1)^α+β+γ under odd transposition. 3-j symbols are 0 if m_γ ≠ m_α + m_β and if the triangular condition α − β ≤ γ ≤ α + β is not satisfied. It is worth to note the following parity property of the 3-j symbol

Similarly 6-j symbols and 9-j symbols appear when combining the angular momentum for three subsystems and four subsystems. 9-j symbols appear in the evaluation of the reduced matrix element of two coupled subsystem as well, see eqn (29). They are defined in terms of 6-j symbols as:

where max(|a − j, |b − f|, |d − h|) ≤ x ≤ min(a + j, b + f, d + h). 9-j symbols are zero if the triangular conditions is not satisfied by any of the rows or of the columns.

3-j symbols and 6-j symbols can be calculated from generating formulae²⁷ and are available as built-in objects in computational softwares such as Mathematica.³⁴

A.3 Evaluation of matrix elements

The evaluation of matrix elements for a function which is a direct product of a radial and a coupled angular part, such as the potential term in eqn (2) and the INS cross section in eqn (10), are evaluated as direct product of a radial and an angular part:

A.3.1 Radial matrix elements. The evaluation of radial matrix elements for non-negative power of local operator f(R) = R^k over the radial translational wave-functions, |n,l;ρ〉 of eqn (17), is found by integrationwhere Γ is the gamma function⁴⁷ andare the coefficients for the Laguerre expansion of L^(α)_n(x) in power of x^m.

Similarly the radial matrix elements of the spherical Bessel function⁴⁷j_β(qR) for integer β ≥ 0, which enter in the space part of the INS transition matrix, can be integrated in a closed form

For sake of completeness, we note that only matrix elements of spherical Bessel functions j_β with |l′ − l| ≤ β ≤ l′ + l and l′ + l + β = even needs to be evaluated. In such case the confluent hypergeometric function ₁F₁ reduces to:⁴⁷where

(a)_n = a(a + 1)(a + 2)…(a + n − 1), (a)₀ = a (27)

is the Pochammer symbol.⁴⁷

A.3.2 Angular matrix elements. The evaluation of the matrix elements of bipolar harmonics over the coupled spherical basis is easily achieved via the use of the Wigner–Eckart theorem⁴⁸where the round brackets represent a 3j symbol and the double bars represents the reduced matrix elements for the bipolar spherical harmonics. The reduced matrix elements in the coupled angular momentum scheme⁴⁸is given in terms of reduced matrix elements for spherical harmonics⁴⁸

The angular matrix elements in eqn (28) are diagonal in λ and m_λ for bipolar harmonics with Λ = m_Λ = 0:

A.4 INS transition probability

In neutron scattering experiments involving transitions between two energy levels in a homonuclear diatomic molecule, labelled here as i and f respectively, the main quantity of interest is the differential cross section in Born approximation³⁰where the sum is over the nuclear species ξ in the system at position x_ξ = R + (−1)^ξ(r/2), k and k′ are the wave-vectors labelling the incident and scattered neutron wave-function, q = k − k′ is the transferred momentum, b̂ is the scattering length operator andFor nuclei with non zero spin I, the scattering length operator is defined by³³where ŝ and Î are the neutron and nuclear spin operators, respectively and b_c and b_i are the bound coherent and incoherent scattering length.

When evaluating the neutron cross-section for homonuclear diatomic molecules, only the symmetric combination intervenes in transitions where the total molecular nuclear spin does not change, ΔI_N = 0, while the antisymmetric combination intervenes in transitions where the total molecular nuclear spin changes of one unit, |ΔI_N| = 1. In experiments with non polarised neutrons and non polarised molecular system, the cross section can be written as product of a spin dependent part s_fi and space dependent part S_fi:

when I^f_N ← Iⁱ_N andwhen I^f_N ± 1 ← Iⁱ_N.

The spatial dependent part of the INS cross section involves the evaluation of the matrix elements of the local operator g_±(q) where q = k − k′ is the transferred momentum. The plus and minus sign are alternatively required for direct (ortho–ortho and para–para) and cross spin transitions (ortho–para and para–ortho), respectively. For our scopes, this operator can be conveniently written as

where * stands for complex conjugate andhaving used the expansion of a plane wave into spherical harmonics and spherical Bessel functions j_α²⁷and the expansion of the product of two spherical harmonics into a single spherical harmonic (Clebsch–Gordan series)²⁷together with the definition of bipolar harmonics, eqn (3). The shorthand notation j_β,±(x) = [j_β(x) ± j_β(−x)]/2 has been introduced in eqn (38).

In systems with spherical symmetry, the space dependent part of the INS cross section for transitions between an initial energy level i = E_λ and a final energy level is given by the transition probability of the operator g_±(q) summed over the final states m_λ′ and averaged over the 2λ + 1 initial states m_λ:

wherehaving used the Wigner–Eckart theorem,⁴⁸ see eqn (28) together with the orthogonality properties of the 3j symbols²⁷and of the spherical harmonics²⁷The condition 20 together with the identity j_β,±(x) = [1 ± (−1)^β]j_β(x)/2 valid for integer β values allows to drop the ± sign from eqn (43) and (44).

Acknowledgements

This work was funded by the UK Engineering and Physical Sciences Research Council (EPSRC), grant number EP/M001970/1. The authors are indebted to Nicholas J. Turro for providing the original samples and Mark Denning for purifying it.

References

1. K. Komatsu, M. Murata and Y. Murata, Science, 2005, 307, 238–240.
2. K. Kurotobi and Y. Murata, Science, 2011, 333, 613–616.
3. A. Krachmalnicoff, R. Bounds, S. Mamone, S. Alom, M. Concistrè, B. Meier, K. Kouřil, M. E. Light, M. R. Johnson, S. Rols, A. J. Horsewill, A. Shugai, U. Nagel, T. Rõõm, M. Carravetta, M. H. Levitt and R. J. Whitby, Nat. Chem., 2016, 1–5.
4. M. H. Levitt, Philos. Trans. R. Soc., A, 2013, 371, 20120429.
5. S. Mamone, M. Ge, D. Hüvonen, U. Nagel, A. Danquigny, F. Cuda, M. C. Grossel, Y. Murata, K. Komatsu, M. H. Levitt, T. Rõõm and M. Carravetta, J. Chem. Phys., 2009, 130, 081103.
6. M. Ge, U. Nagel, D. Hüvonen, T. Rõõm, S. Mamone, M. H. Levitt, M. Carravetta, Y. Murata, K. Komatsu, J. Y.-C. Chen and N. J. Turro, J. Chem. Phys., 2011, 134, 054507.
7. A. J. Horsewill, S. Rols, M. R. Johnson, Y. Murata, M. Murata, K. Komatsu, M. Carravetta, S. Mamone, M. H. Levitt, J. Y.-C. Chen, J. A. Johnson, X. Lei and N. J. Turro, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081410.
8. A. J. Horsewill, K. Goh, S. Rols, J. Ollivier, M. R. Johnson, M. H. Levitt, M. Carravetta, S. Mamone, Y. Murata, J. Y.-C. Chen, J. A. Johnson, X. Lei and N. J. Turro, Philos. Trans. R. Soc., A, 2013, 371, 20110627.
9. T. Pham, K. A. Forrest, B. Space and J. Eckert, Phys. Chem. Chem. Phys., 2016, 18, 17141–17158.
10. Z. Bačić and J. C. Light, Annu. Rev. Phys. Chem., 1989, 40, 469–498.
11. M. Xu, Y. S. Elmatad, F. Sebastianelli, J. W. Moskowitz and Z. Bačić, J. Phys. Chem. B, 2006, 110, 24806–24811.
12. F. Sebastianelli, M. Xu, Y. S. Elmatad, J. W. Moskowitz and Z. Bačić, J. Phys. Chem. C, 2007, 111, 2497–2504.
13. F. Sebastianelli, M. Xu, D. K. Kanan and Z. Bačić, J. Phys. Chem. A, 2007, 111, 6115–6121.
14. M. Xu, F. Sebastianelli and Z. Bačić, J. Phys. Chem. A, 2007, 111, 12763–12771.
15. M. Xu, F. Sebastianelli, Z. Bačić, R. Lawler and N. J. Turro, J. Chem. Phys., 2008, 129, 064313.
16. M. Xu, F. Sebastianelli, B. R. Gibbons, Z. Bačić, R. Lawler and N. J. Turro, J. Chem. Phys., 2009, 130, 224306.
17. S. Ye, M. Xu, Z. Bačić, R. Lawler and N. J. Turro, J. Phys. Chem. A, 2010, 114, 9936–9947.
18. M. Xu, L. Ulivi, M. Celli, D. Colognesi and Z. Bačić, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 241403.
19. M. Xu and Z. Bačić, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 1–10.
20. M. Xu, S. Ye, A. Powers, R. Lawler, N. J. Turro and Z. Bačić, J. Chem. Phys., 2013, 139, 064309.
21. M. Xu, M. Jiménez-Ruiz, M. R. Johnson, S. Rols, S. Ye, M. Carravetta, M. S. Denning, X. Lei, Z. Bačić and A. J. Horsewill, Phys. Rev. Lett., 2014, 113, 123001.
22. S. Mamone, J. Y.-C. Chen, R. Bhattacharyya, M. H. Levitt, R. G. Lawler, A. J. Horsewill, T. Rõõm, Z. Bačić and N. J. Turro, Coord. Chem. Rev., 2011, 255, 938–948.
23. A. Ivanov, M. Jimenéz-Ruiz and J. Kulda, J. Phys.: Conf. Ser., 2014, 554, 012001.
24. A. J. Horsewill, K. S. Panesar, S. Rols, J. Ollivier, M. R. Johnson, M. Carravetta, S. Mamone, M. H. Levitt, Y. Murata, K. Komatsu, J. Y.-C. Chen, J. A. Johnson, X. Lei and N. J. Turro, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 205440.
25. S. Mamone, PhD thesis, University of Southampton, UK, 2011.
26. B. Poirier, J. Chem. Phys., 2015, 143, 101104.
27. D. A. Varshalovich, A. Moskalev and V. Khersonskii, Quantum Theory of Angular Momentum, World Scientific, Singapore, 1988, p. 514.
28. S. Mamone, G. Pileio and M. H. Levitt, Symmetry, 2010, 2, 1423–1449.
29. W. G. Harter and D. E. Weeks, J. Chem. Phys., 1989, 90, 4727.
30. S. W. Lovesay, Theory of Neutron scattering from Condensed Matter, Oxford University Press, 1984, vol. 1, p. 329.
31. M. Xu, S. Ye and Z. Bačić, J. Phys. Lett., 2015, 6, 3721–3725.
32. G. Herzberg, Molecular Spectra and Molecular Structure – I. Spectra of Diatomic Molecules, Van Nostrand, 1950.
33. V. F. Sears, Neutron News, 1992, 3, 26–37.
34. Mathematica, Wolfram Research Inc., 2015.
35. G. Herzberg and L. L. Howe, Can. J. Phys., 1959, 37, 636–659.
36. U. Fink, T. Wiggins and D. Rank, J. Mol. Spectrosc., 1965, 18, 384–395.
37. T. Yildirim and A. B. Harris, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 214301.
38. P. M. Felker and Z. Bačić, J. Chem. Phys., 2016, 145, 084310.
39. S. Mamone, M. R. Johnson, J. Ollivier, S. Rols, M. H. Levitt and A. J. Horsewill, Phys. Chem. Chem. Phys., 2016, 18, 1998–2005.
40. L. Ulivi, M. Celli, A. Giannasi, A. J. Ramirez-Cuesta, D. J. Bull and M. Zoppi, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 161401.
41. K. T. Tait, F. Trouw, Y. Zhao, C. M. Brown and R. T. Downs, J. Chem. Phys., 2007, 127, 134505.
42. J. M. Nicol, J. Eckert and J. Howard, J. Phys. Chem., 1988, 92, 7117–7121.
43. J. Eckert, J. M. Nicol, J. Howard and F. R. Trouw, J. Phys. Chem., 1996, 100, 10646–10651.
44. M. H. Rosnes, M. Opitz, M. Frontzek, W. Lohstroh, J. P. Embs, P. A. Georgiev and P. D. C. Dietzel, J. Mater. Chem. A, 2015, 3, 4827–4839.
45. T. Pham, K. A. Forrest, E. H. L. Falcão, J. Eckert and B. Space, Phys. Chem. Chem. Phys., 2016, 18, 1786–1796.
46. S. Flügge, Practical Quantum Mechanics, Springer, Berlin, Heidelberg, 1971, p. 620.
47. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Issue 55, ed. M. Abramowitz and I. A. Stegun, U.S. Government Printing Office, 10th edn, 1964, p. 1046.
48. R. N. Zare, Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics, John Wiley & Sons, 1988, p. 368.

---