## Abstract

We report rigorous quantum calculations of the inelastic neutron scattering (INS) spectra of HD@C_{60}, over a range of temperatures from 0 to 240 K and for two incident neutron wavelengths used in recent experimental investigations. The computations were performed using our newly developed methodology, which incorporates the coupled five-dimensional translation–rotation (T–R) eigenstates of the guest molecule as the initial and final states of the INS transitions, and yields highly detailed spectra. Depending on the incident neutron wavelength, the number of computed INS transitions varies from almost 500 to over 2000. The low-temperature INS spectra display the fingerprints of the coupling between the translational and rotational motions of the entrapped HD molecule, which is responsible for the characteristic splitting patterns of the T–R energy levels. INS transitions from the ground T–R state of HD to certain sublevels of excited T–R multiplets have zero intensity and are absent from the spectra. This surprising finding is explained by the new INS selection rule introduced here. The calculated spectra exhibit strong temperature dependence. As the temperature increases, numerous new peaks appear, arising from the transitions originating in excited T–R states which become populated. Our calculations show that the higher temperature features typically comprise two or more transitions close in energy and with similar intensities, interspersed with numerous other transitions whose intensities are negligible. This implies that accurately calculated energies and intensities of INS transitions which our methodology provides will be indispensable for reliable interpretation and assignment of the experimental spectra of HD@C_{60} and related systems at higher temperatures.

## 1. Introduction

The endohedral complexes of C_{60} with H_{2} and HD have received a great deal of attention [1] because of the profoundly quantum character of the dynamics of the guest molecule, which is surprisingly complex despite the apparent simplicity of these systems [2,3]. Large quantum effects dominate the translational motions of the centre of mass of H_{2}/HD as well as the rotation of the molecule; they arise from the combination of the lightness of H_{2} or HD and the small size of the C_{60} cavity. The three translational degrees of freedom of the incarcerated H_{2}/HD are quantized and their eigenstates are well separated in energy. The same holds for the quantized rotational levels of H_{2}/HD owing to their exceptionally large rotational constants, 7.24 and 5.52 meV in the ground vibrational state, inside C_{60} [4]. This gives rise to a very sparse translation–rotation (T–R) energy level structure.

H_{2}@C_{60} and HD@C_{60} differ in one profound aspect. For the caged H_{2}, whose two protons with the nuclear spin are fermions, the Pauli principle requires the total molecular wave function to be antisymmetric with respect to the nuclear exchange. This gives rise to two nuclear spin isomers, *para*-H_{2} (*p*-H_{2}) and *ortho*-H_{2} (*o*-H_{2}), having the antisymmetric *I*=0 and the symmetric *I*=1 total nuclear spin states, respectively. For *p*-H_{2}, only even rotational quantum numbers *j* are allowed (*j*=0,2,…), whereas *o*-H_{2} can have exclusively odd rotational quantum numbers *j*=1,3,…. In the absence of a catalyst, the interconversion between the two nuclear spin isomers is extremely slow. This means that inside C_{60} *p*-H_{2} and *o*-H_{2} do not equilibrate, and the population of *o*-H_{2} in its *j*=1 ground state persists at very low temperature close to 0 K [1]. By contrast, the Pauli exclusion principle does not apply to the heteronuclear HD molecule, which therefore has both even- and odd-*j* rotational states. At temperatures below 20 K, only the ground rotational state (*j*=0) of HD is populated.

Our quantum five-dimensional calculations of the T–R eigenstates of H_{2} and HD in C_{60} [2,3,5] have characterized quantitatively the ‘rattling’ dynamics of the endohedral hydrogen. They have revealed strong coupling between the orbital angular momentum associated with the translational motions of H_{2}/HD and the rotational angular momentum of the molecule. This T–R coupling manifests conspicuously through the entire T–R energy level structure, partially lifting the degeneracies, and splitting the eigenstates, of H_{2} and HD in C_{60} with simultaneous translational and rotational excitation [2,3,5]. Shortly thereafter, the first infrared (IR) spectra measured for H_{2}@C_{60} [6] showed absorption lines split into distinct peaks owing to the T–R coupling, as predicted by our calculations [2,3]. Subsequent IR spectroscopic studies of H_{2}@C_{60} [4] as well as of endohedral HD and D_{2} in C_{60} [7] have provided a great deal of valuable information about the T–R energy levels of the guest molecules and their interaction potentials with C_{60}. The physical properties of hydrogen endofullerenes have also been probed by nuclear magnetic resonance [1,8,9]. What has enabled the experimental investigations is the synthetic organic breakthrough, known as ‘molecular surgery’, by Komatsu and co-workers [10,11] for producing macroscopic amounts of H_{2}@C_{60} and HD@C_{60}.

Inelastic neutron scattering (INS) spectroscopy has proved to be a powerful and versatile tool for elucidating with atomic-level resolution the behaviour of molecular hydrogen entrapped in a variety of nanoscale cavities [12,13]. Horsewill *et al.* have pioneered the use of INS to probe the quantum dynamics of H_{2} and HD inside an open-cage fullerene [14] and in C_{60} [15,16]. These studies have determined the energies of the majority of the lower lying T–R eigenstates of H_{2}@C_{60} [16]. In their interpretation and assignment of the measured INS spectra, Horsewill and co-workers have been guided by the general theoretical framework developed by us [2,3], in particular the T–R coupling scheme, and the accurately calculated T–R energy levels of H_{2}@C_{60} [5].

The analysis and assignment of the experimental INS spectra would be more complete and reliable if the theory could provide not only accurate T–R excitation energies but also the intensities of the INS transitions. Recently, we have developed the quantum methodology that achieves this goal by fully including the T–R coupling [17,18]. It allows rigorous quantum simulation of the INS spectra, i.e. the energies and the intensities of the INS transitions, of a hydrogen molecule confined inside a nanoscale cavity with an arbitrary shape. The distinguishing feature and the main strength of our approach is its incorporation of the coupled quantum five-dimensional T–R energy levels and wave functions of the entrapped hydrogen molecule as the initial and final states of the INS transitions [17,18]. As a result, the simulated INS spectra fully capture the intricacies of the T–R dynamics of the guest molecule in nanoconfinement and are exceptionally realistic, as demonstrated by the excellent agreement between the calculated and measured INS spectra of *p*- and *o*-H_{2} inside the small cage of the structure II clathrate hydrate [17].

In the present paper, this methodology is applied to the calculation of the INS spectra of HD@C_{60} and exploring their temperature dependence. The quantum T–R dynamics of HD@C_{60} has been shown to differ significantly from that of the caged H_{2}, on account of the mass anisotropy of HD and the significant rotational state mixing it causes [3]. The INS spectroscopy has so far provided only very limited experimental information about this interesting system [15,16]. Therefore, the predicted spectra which have emerged from this study should be very valuable for future INS investigations of HD@C_{60} [17].

## 2. Theory

### (a) Calculation of the coupled five-dimensional translation–rotation eigenstatesof HD@C_{60}

Our approach to the rigorous calculation of the coupled five-dimensional T–R energy levels and wave functions of a hydrogen molecule inside fullerenes, applied in this study to HD@C_{60}, has evolved in our group over a number of years, in the course of the theoretical investigations of the quantum T–R dynamics of H_{2} and its isotopologues entrapped in the cages of clathrate hydrates [19–21], fullerenes C_{60} and C_{70} [2,5], and an open-cage derivative of C_{60} (ATOCF) [22]. A detailed description of the quantum five-dimensional methodology is available in [18], and therefore only the salient features are given here. C_{60} is treated as rigid, and its geometry used in our calculations has been determined experimentally, from the gas-phase electron diffraction study [23]. The set of five coordinates (*x*,*y*,*z*,*θ*,*ϕ*) is used; *x*,*y* and *z* are the Cartesian coordinates of the centre of mass of H_{2}, whereas the two polar angles *θ* and *ϕ* specify its orientation. The coordinate system is aligned with the three principal axes of the fullerene, and its origin is at the centre of mass of the cage. The three rotational constants of C_{60}, a spherical top, are very small and equal to 2.803×10^{−3} cm^{−1}, which justifies treating the fullerene as non-rotating. In this case, the five-dimensional Hamiltonian for the T–R motions of the caged diatomic molecule is
2.1In equation (2.1), *μ* is the reduced mass of HD in C_{60} (3.0094 amu), and **j**^{2} is the angular momentum operator of the diatomic. *B* denotes the rotational constant of the endohedral molecule whose bond length is held fixed, and will be discussed shortly. *V* (*x*,*y*,*z*,*θ*,*ϕ*) in equation (2.1) is the five-dimensional intermolecular potential energy surface (PES) for the interaction between the entrapped HD and the interior of C_{60}, which is described below. The T–R energy levels and wave functions of the Hamiltonian in equation (2.1) are obtained utilizing the efficient computational methodology that originated in our group [19,24]. The final Hamiltonian matrix, its size drastically reduced by the sequential diagonalization and truncation procedure [25], is diagonalized, yielding the T–R eigenstates which are numerically exact for the five-dimensional PES used.

The parameters extracted from the IR spectra of HD@C_{60} [7], listed in table 1, give the rotational constant *B*_{0}=44.53 cm^{−1} (5.521 meV) for HD in the ground vibrational state *v*=0. This value of *B*_{0}, together with the centrifugal distortion parameter *D*_{e} in table 1, was used in the quantum five-dimensional calculations of the T–R eigenstates of HD (*v*=0) in C_{60} reported here. In these calculations, the dimension of the sinc-discrete variable representation (DVR) basis was 15 for each of the three Cartesian coordinates *x*, *y* and *z*, spanning the range −2.83 bohr ≤*λ*≤ 2.83 bohr (*λ*=*x*,*y*,*z*). The angular basis included the even- and odd-*j* rotational functions up to . The cut-off parameter for the size of the intermediate three-dimensional eigenvector basis was set to 468, resulting in the final five-dimensional Hamiltonian matrix of dimension 29 952. These basis set parameters were chosen following extensive testing, assuring that the T–R energy levels reported in this paper are converged to at least five significant figures.

### (b) Calculation of the inelastic neutron scattering spectra of HD@C_{60}

The methodology for accurate and efficient quantum calculation of the INS spectra, both the energies and intensities of the transitions, of a hydrogen molecule confined inside a nanoscale cavity has been introduced recently [17]. It uses the five-dimensional T–R eigenstates of the Hamiltonian in equation (2.1) as the initial and final states of the INS transitions. For the case of the entrapped homonuclear molecule such as H_{2}, the details of the formalism are given in [18]; its extension to the caged heteronuclear diatomics such as HD is presented in [26]. In the following, we elaborate on one aspect of the quantum simulation of the INS spectra which has not been addressed previously.

Let *E* and *E*^{′} be the energies of the incident and scattered neutron, respectively, and ** k** and

*k*^{′}their respective wavevectors. Then, the neutron energy transfer can be written as , and the neutron momentum transfer as . The double differential cross section for neutron scattering in the first Born approximation is [27,28] 2.2where 2.3and 2.4In equations (2.2)–(2.4), |

*i*〉 stands for the initial state of the scattering molecular system with the energy

*ϵ*

_{i},

*p*

_{i}is its statistical weight, |

*f*〉 is the final state with the energy

*ϵ*

_{f}, is the scattering length operator [18] and

*r*_{n}is the position of nucleus

*n*. In our case, |

*i*〉 and |

*f*〉 are represented by the T–R eigenfunctions of the Hamiltonian in equation (2.1). The evaluation of

*S*(

**,**

*κ**ω*) in equation (2.3) using the quantum five-dimensional T–R eigenstates is described in [18].

The experimental INS spectra of a hydrogen molecule in C_{60} are taken from powdered samples [15,16], in which the fullerene cages are randomly oriented with respect to the incoming neutron beam. Therefore, in order to achieve a more realistic comparison with the measured spectra, the computed INS spectra are averaged over all possible orientations of C_{60}, as described previously [18]. Moreover, with the IN4C spectrometer used by Horsewill and co-workers [16], the INS spectra of H_{2} and HD in C_{60} are recorded over a broad range of *κ*≡|** κ**| values. This means that

*S*(

*κ*,

*ω*), obtained by averaging

*S*(

**,**

*κ**ω*) in equation (2.3) over all directions of the wavevector

**[18], needs to be integrated also over**

*κ**κ*to obtain the desired

*S*(

*ω*), which is done below.

We denote the angle between ** k** and

*k*^{′}, the wavevectors of the incident and scattered neutrons, as

*θ*

_{0}≡

*θ*(

**,**

*k*

*k*^{′}). It is easy to derive the relationship between

*κ*, on the one hand, and the neutron energy transfer

*ΔE*=

*E*−

*E*

^{′}and the angle

*θ*

_{0}, on the other hand: 2.5From it, 2.6Given the

*κ*range , based on equation (2.6), the following expression is obtained for the angular range : 2.7where 2.8The integration of

*S*(

*κ*≡

*κ*

_{θ0},

*ω*) over

*κ*is replaced by the integration over

*θ*

_{0}from to in equations (2.7) and (2.8). From equation (2.5), one can obtain 2.9where

*k*≡|

**| and**

*k**k*

^{′}≡|

*k*^{′}|. Then, 2.10

The integration of the INS spectrum over both *κ* and the angles (*θ*_{κ},*ϕ*_{κ}), which define the direction of ** κ** (to account for the random orientations of the C

_{60}cages), is performed for each distinct INS transition in the range of excitation energies considered. For example, at the incident neutron wavelength

*λ*

_{n}=1.1 Å, there are 2623 such transitions. This procedure is computationally very time consuming, but it is worthwhile, because it results in an extremely realistic simulation of the experimental INS spectra. The integration in equation (2.10) is performed using the 20-point Gauss–Legendre quadrature. The numerical integration involved in the powder averaging is carried out using the 10-point Gauss–Legendre quadrature in

*θ*

_{κ}and the 20-point Gauss–Chebyshev quadrature in

*ϕ*

_{κ}.

### (c) Potential energy surface

As in all our previous theoretical studies of a hydrogen molecule inside fullerenes, the five-dimensional intermolecular PES *V* _{H2−fullerene} between the confined hydrogen (H_{2}/HD) molecule and the carbon atoms of C_{60} is assumed to be pairwise additive,
2.11where **q** are the coordinates (*x*,*y*,*z*,*θ*,*ϕ*) of the endohedral H_{2}/HD defined above, *V* _{H2−C} is the pair interaction specified below between H_{2}/HD and a carbon atom of the fullerene, and the index *k* runs over all C atoms, whose coordinates *Ξ*_{k} are fixed. *V* _{H2−C} is represented by the three-site pair potential introduced by us [5]. Its parameters used in this paper are given in table 1. They were determined by achieving the best match between the calculated T–R energy levels of H_{2} (*v*=0) in C_{60} and those obtained by fitting the IR spectroscopic data [7] for this system.

## 3. Results and discussion

### (a) Translation–rotation energy levels of HD@C_{60}

The T–R energy levels of HD (*v*=0) in C_{60} from the quantum five-dimensional calculations on the optimized PES are shown in table 2, together with their quantum number assignments (*n*,*j*,*λ*,*l*) and other useful information. They exhibit all the key features identified previously [2,3]: (i) the orbital angular momentum **l** of the centre of mass of HD and the rotational angular momentum **j** of HD couple to give the total angular momentum ** λ**=

**l**+

**j**, whose values are

*λ*=

*l*+

*j*,

*l*+

*j*−1,…,|

*l*−

*j*|, with the degeneracy of 2

*λ*+1; (ii) the integer values of

*l*are those allowed for the quantum number

*n*of the three-dimensional isotropic harmonic oscillator (

*l*=

*n*,

*n*−2,…,1 or 0, for odd or even

*n*, respectively), whereas the quantum number

*j*of the rigid rotor takes the values

*j*=0,1,2,…; (iii) the T–R states with the same quantum numbers

*n*and

*j*are split into as many distinct levels as there are different values of

*λ*, each with the degeneracy 2

*λ*+1; (iv) the translationally excited states are not harmonic, and consequently their energies depend not only on

*n*but also on

*l*; (v) the T–R levels with

*λ*=3,4 appear as closely spaced pairs of levels with the degeneracies 3 and 4 for

*λ*=3, and 4 and 5 for

*λ*=4, a manifestation of the ‘crystal field’ splitting by the icosahedral

*I*

_{h}environment of C

_{60}.

We highlight the T–R levels *i*= 3–5 in table 2, which have *n*=*j*=1, and owing to the T–R coupling [2,3] are split into three *λ* sublevels. This type of splitting has received considerable experimental scrutiny, in both INS [15,16] and IR spectroscopy [6]. Interestingly, in the case of HD@C_{60}, the ordering of the *λ* sublevels in ascending energy is *λ*=2,0,1, whereas, for H_{2}@C_{60}, the energy ordering is *λ*=1,2,0 [2,3,5]. Thus, the energy ordering of the *λ* sublevels is *not* the same for HD and H_{2}.

The T–R energy levels of HD@C_{60} in table 2 display another property predicted by the earlier calculations [3]—that at higher excitation energies, *j* is generally no longer a good quantum number for the incarcerated HD, in contrast to the homonuclear isotopologues H_{2} and D_{2} [2,3]. This is evident from the fact that for many of the computed T–R levels of HD, especially the higher lying ones, the contribution *c*(*j*) of the dominant rotational basis function *j* is often only slightly greater than 0.5, and sometimes even smaller. In other words, there is extensive mixing of two or more basis functions with different *j* values, which has been attributed to the pronounced mass anisotropy of HD [3]. The manifestations of the predicted strongly mixed character of the T–R states of HD in C_{60} have been observed in the INS spectra of this system [15].

At higher excitation energies, strong mixing of the rotational basis functions eventually precludes the assignment of the T–R eigenstates. This is evident in table 2, where starting with the level *i*=18, the quantum numbers *n*,*j* and *l* cannot be determined with confidence. The only remaining good quantum number is *λ* [3], which can be established from the degeneracies of the levels. The exception in this energy range is the closely spaced pair of states *i*=28,29, which are highly pure *j*=3 rotational levels of HD; their assignment is (0,3,3,0).

### (b) Low-temperature inelastic neutron scattering spectrum and a new inelastic neutron scattering selection rule

The INS spectrum of HD@C_{60} calculated for 0 K and the incident neutron wavelength *λ*_{n}=1.1 Åis displayed in figure 1. This wavelength was chosen because it has been used in the INS spectroscopy of H_{2} in C_{60} [16], and because it can probe a broad range of T–R excitations up to about 70 meV. Shown on the *x*-axis is the neutron energy (NE) transfer *ΔE*=*E*−*E*^{′}, where *E* and *E*^{′} are the energies of the incident and scattered neutron, respectively. *ΔE* is positive in NE loss and negative in NE gain; hence, the INS spectrum shown is in NE loss. The computed stick spectra are convolved with the instrumental resolution function. At this low temperature, all transitions originate solely in the ground T–R state (0,0,0,0) of HD because it is the only one populated. They are labelled with the (*n*,*j*,*λ*,*l*) quantum numbers of their respective final states. The calculated energies and intensities of the INS transitions which appear in figure 1 are listed in table 3.

The peak at 10.0 meV corresponds to the purely rotational excitation of HD, namely the transition , whereas the neighbouring peak at 18.7 meV arises from the fundamental translational excitation of HD. These two peaks have been observed in the measured INS spectra of HD@C_{60} at 10.0 and 19.1 meV, respectively [15,16].

The partially resolved band centred at 26 meV is particularly interesting. It is associated with the transitions from the ground T–R state of HD to its *n*=1,*j*=1 triplet of (1,1,*λ*,1) levels which, in the order of increasing energies, have the total angular momentum quantum numbers *λ*=2,0,1, respectively (table 2). Therefore, one expects to see three transitions in the calculated spectrum, corresponding to , *λ*=(2,0,1). But, instead, only two transitions appear in figure 1, at 25.2 meV and at 27.8 meV. The transition at 28.3 meV to the *λ*=1 sublevel is missing, and table 3 shows that its intensity is zero. This is not accidental; this transition is *strictly forbidden* according to the novel and unexpected *selection rule* that we present here:

*Consider a* *p*-H_{2} *or* HD *molecule confined in a near-spherical environment, such as that of* C_{60}, *where its orbital angular momentum* *l**and the rotational angular momentum* *j**couple to give the total angular momentum* ** λ**,

*with the eigenvalues*

*λ*=

*l*+

*j*,

*l*+

*j*−1,…,|

*l*−

*j*|.

*Then, the INS transition involving the ground T–R state*

**(0,0,0,0)**

*of*

*p*-H

_{2}

*or*HD

*and the excited state which is assignable as*(

**,**

*n***,**

*j***,**

*λ***) is**

*l***forbidden**

*for*

**=**

*λ***+**

*j***−**

*l***1**,

*in both NE loss and NE gain*.

Clearly, the state (1,1,1,1) satisfies the condition *λ*=*j*+*l*−1 of the selection rule. Rigorous derivation of this selection rule is rather lengthy and involved, and will be presented elsewhere in the near future.

The splittings among the three *λ* sublevels of *n*=1,*j*=1 are highly uneven; according to tables 2 and 3, *λ*=2 and *λ*=0 sublevels are split by 2.6 meV, whereas *λ*=0 and *λ*=1 sublevels are only 0.49 meV apart. Given that the transition to the *λ*=1 sublevel is forbidden by the above selection rule, we predict that the width of this band in the measured INS spectrum of HD@C_{60} will reflect the calculated splitting of the *λ*=2 and *λ*=0 sublevels. Clearly, the theoretical information regarding the relative intensities of the transitions to the *λ*=2,0,1 triplet will be crucial for the correct interpretation of the experimental INS spectrum.

The peak at 33.5 meV corresponds to the HD rotational transition, .

The asymmetric band peaking at 39.3 meV is due to the two two-quanta translational excitations, at 39.3 meV and at 41.1 meV; the former transition is about four times more intense. The fact that the two *n*=2 states have different energies for *l*=0 and *l*=2 is evidence of their anharmonicity [2].

Two small peaks at 44.5 and 47.5 meV arise from two weak transitions to the *n*=2,*j*=1 multiplet of HD split by the T–R coupling (table 2), and , respectively. The (also) weak transition to the *λ*=3 component of this multiplet, at 41.6 meV, is well inside the 39.3 meV band, close to the much stronger transition to the (2,0,0,0) state, and will be difficult to identify in the measured INS spectrum. Missing from the computed spectrum is the transition to the *λ*=2 component, , at 45.1 meV. Its intensity is zero (table 3), as this transition is also forbidden according to the selection rule above.

The peaks at 51.1 and 54.1 meV correspond to the transitions to the *λ*=3 and *λ*=1 sublevels of the *n*=1,*j*=2 multiplet of HD, and , respectively. Again, one transition, to the *λ*=2 component of this multiplet which should be at 51.8 meV (table 3), has zero intensity and does not appear in the spectrum. It is yet another example of a forbidden transition.

Two bands centred at 60 and 65 meV comprise several transitions to final states which cannot be assigned (apart from the quantum number *λ*; table 2), with the exception of the *j*=3 rotational state (0,3,3,0) at 65.9 meV.

### (c) Temperature dependence of the inelastic neutron scattering spectra: NE loss

The INS spectra exhibit a strong temperature dependence. The reason for this lies in the expressions for the INS scattering cross sections in equations (2.2)–(2.4), which involve summation over the initial states of the system weighted by their relative populations according to the Boltzmann distribution for the given temperature. While at low temperatures only the ground T–R state of HD@C_{60} is significantly populated, as the temperature increases a growing number of excited T–R states becomes populated as well. They can serve as the initial states of the INS transitions, thus greatly enlarging the number of transitions that can give rise to new peaks in the spectrum. The density of the INS transitions at higher temperature (10 or more per meV) is much greater than close to 0 K, and they blend into broad bands which are often only partly resolved, or not at all. As shown below, most features in the higher temperature INS spectra have contributions from two or more energetically close transitions with comparable intensities. In addition, numerous other transitions whose computed intensities are very small lie in the same energy range. All this makes the assignment of the INS spectra recorded at elevated temperature much more challenging than assigning the spectra taken near 0 K, and also much more reliant on theory.

Figure 2 shows the INS spectra of HD@C_{60} in NE loss calculated at 120 and 240 K for the incident neutron wavelength *λ*_{n}=1.1 Å, together with the 1.6 K spectrum identical to that at 0 K in figure 1. It is clear that raising the temperature causes large changes in the overall appearance of the INS spectra. The relative size of the peaks originating in the ground T–R state decreases as its population transfers to excited T–R levels, and between them new features arise from transitions originating in higher lying T–R eigenstates. The more prominent new features in the 240 K spectrum are labelled with numbers **1–12** in figure 2 and the INS transitions mostly responsible for them are listed in the figure caption.

Feature **1** is due to a single transition at 5.2 meV. But already feature **2** arises from two very different transitions with similar intensities, which happen to be close in energy, at 14.9 meV and at 15.2 meV. The picture grows more complicated with higher energies. Thus, feature **3** comprises four T–R transitions that have significant intensities, with energies from 19.9 to 21.3 meV; for some final states, the quantum numbers could not be assigned, except for *λ*. This is characteristic also for the other features. It is clear that the ability of theory to compute accurate INS intensities, in addition to their energies, is crucial for correct interpretation and assignment of the higher temperature INS spectra.

### (d) Temperature dependence of the inelastic neutron scattering spectra: NE gain

Figure 3 shows the INS spectra of HD@C_{60}, mostly in NE gain, calculated at 20, 80, 120 and 240 K for the incident neutron wavelength *λ*_{n}=3.0 Å. This wavelength has been used to record the INS spectra of H_{2} and HD in C_{60} in NE gain [15]. Some of the peaks in the NE gain arise from the reversals of the transitions present in NE loss at low temperature, but many more are associated with previously unobserved transitions originating in higher lying T–R states.

The peak at −10.0 meV arises from the rotational transition of HD, ; it is the NE gain partner of the peak at 10.0 meV in NE loss at 0 K (figure 1). At 20 K, the −10.0 meV peak is the only one visible in NE gain. It is barely observable because of the very low population of the excited (0,1,1,0) state, but its amplitude grows rapidly with increasing temperature. This peak is present in the INS spectra of HD@C_{60} measured at 120 K [15]. The slight asymmetry of this peak apparent at 240 K is due to the transition at −9.6 meV, which attains significant intensity at this temperature.

The prominent peak at −19 meV is in part the NE gain partner of the peak in NE loss at 18.7 meV corresponding to . However, in NE gain, besides the reverse transition at −18.7 meV, another transition, at −18.2 meV, contributes to this peak. The ratio of the intensities of these two transitions varies strongly with temperature, from 3 : 1 at 80 K it drops to 1 : 1 at 240 K, reflecting the changes in the relative populations of the initial states of the two transitions.

The remainder of this broad band, which at 240 K extends to −46 meV, is likely to be poorly resolved in the experimental INS spectra. The transitions responsible for some discernible features at 240 K are indicated in figure 3. For example, a shoulder around −24 meV corresponds to two transitions, the dominant at −23.48 meV and the weaker transition at −23.52 meV. Two weak features in the tail of this band arise from the transitions at −33.5 meV and at −39.3 meV, respectively; both represent reversals of the transitions visible in NE loss at 0 K (figure 1).

Between −10 and 0 meV, two smaller peaks are visible at 240 K, and the transitions which give rise to them are shown in figure 3. The peak at −6 meV comprises mainly two transitions, at −6.5 meV and at −5.2 meV. Their relative intensities change significantly with temperature, from 3 : 1 at 80 K to 1.3 : 1 at 240 K, as the populations of the respective initial states shift.

Finally, the NE loss region between 0 and 10 meV shows four features at 240 K whose assignments are given in the caption of figure 3.

It should be noted that for *ΔE* between 0 and −50 meV there are nearly 300 distinct INS transitions, of which at any temperature only a small fraction have intensities that are potentially observable in the INS spectra. Consequently, in the absence of reliable theoretical information about the INS intensities, assigning the features present in the experimental INS spectra is bound to be plagued with uncertainties.

## 4. Conclusions

We have reported rigorous quantum simulations of the INS spectra of HD@C_{60} in both NE loss and NE gain, and how they vary with temperature. They were performed using the sophisticated methodology for calculating the INS spectra of a hydrogen molecule inside a nanoscale cavity, which has been introduced by us recently [17,18]. The extension of the methodology to heteronuclear guest molecules such as HD, which was required for the computations in this work, has also been presented 26. The procedure is computationally very demanding, because calculating the intensity of each of the several thousand distinct INS transitions considered involves integration over the neutron momentum transfer (*κ*) values, as well as over all directions of the wavevector ** κ**, in order to achieve the most realistic simulation possible of the INS spectra taken from powdered samples [18]. The T–R energy levels and wave functions, which serve as input for the INS spectral simulations, were calculated on the pairwise additive five-dimensional PES whose parameters were optimized by fitting to the IR spectroscopic data [7] for H

_{2}(

*v*=0) in C

_{60}.

The low-temperature INS spectrum of HD@C_{60} in NE loss, calculated for 0 K and the incident neutron wavelength *λ*_{n}=1.1 Å, covers the transitions from the ground T–R state of HD to excited T–R eigenstates with energies up to 70 meV. All transitions can be assigned with the quantum numbers (*n*,*j*,*λ*,*l*) of their final states. The exception are the transitions above 60 meV, for which *n*,*j* and *l* cannot be determined with confidence owing to strong mixing of the rotational basis functions, and *λ* remains the only good quantum number.

The low-temperature INS spectrum in NE loss shows a clear imprint of the T–R coupling scheme developed by us [2,3]. Most of the peaks arise from the transitions originating in the ground T–R state of HD to the components of the multiplets with (*n*=1,*j*=1), (*n*=2,*j*=1) and (*n*=1,*j*=2), whose degeneracy is split as a result of the T–R coupling into patterns which are well understood [2,3]. However, the INS calculations have resulted in a new and unexpected finding. For each of the multiplets above, a transition to one of the sublevels is missing from the computed spectrum, because its intensity is exactly zero. It turns out that these particular transitions are strictly forbidden according to the novel and elegant selection rule formulated in this paper; its derivation is complicated and lengthy, and will be presented in a future publication.

INS spectra have also been calculated at 120 and 240 K in NE loss, and at 20, 80, 120 and 240 K in NE gain. They reveal strong temperature dependence of the spectra in both energy domains. With increasing temperature, numerous new features associated with transitions from a growing number of excited states which become populated appear in the INS spectra and gain in intensity, whereas amplitudes of the peaks originating in the ground T–R state diminish. Our calculations show that most peaks in the higher temperature INS spectra are associated with two, three or more transitions which are close in energy and often have similar intensities. Moreover, these transitions are typically interspersed with others whose intensities are small or negligible. In these circumstances, reliable interpretation and assignment of the experimental spectra critically depends on the ability of theory to provide accurate energies and intensities of the INS transitions.

We expect that the accurate INS spectra computed in this study for a wide range of temperatures will provide valuable guidance to experimental investigations of HD@C_{60} in the future.

The novel quantum methodology for computing the INS spectra [17,18]; 26, whose power has been demonstrated in this paper, is applicable to a broad range of systems of both fundamental and practical significance, where a hydrogen molecule is confined inside the nanocavities of the materials. Besides H_{2}@C_{60}, the theoretical study of which will be reported shortly, this includes molecular hydrogen in clathrate hydrates [19,29,30], metal–organic frameworks [31,32,33,34,35] and zeolites [36,37], among others. Finally, methodological extensions are in progress in our group to the INS spectra of nanoconfined polyatomic molecules such as methane [38] and water [39].

## Funding statement

The computational resources at NYU used in this work were supported in part by the NSF MRI (grant no. CHE-0420810).

## Acknowledgements

Z.B. and N.J.T. thank the NSF for generous support of this research through grants nos. CHE-1112292 and CHE-1111398, respectively. Z.B. has also been supported by the China 111 Project B1202.

## Footnotes

One contribution of 13 to a Theo Murphy Meeting Issue ‘Nanolaboratories: physics and chemistry of small-molecule endofullerenes’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.