In this perspective, I review the current status of the theoretical investigations of the quantum translation-rotation (TR) dynamics and spectroscopy of light molecules encapsulated inside fullerenes, mostly C_{60} and C_{70}. The methodologies developed in the past decade allow accurate quantum calculations of the TR eigenstates of one and two nanoconfined molecules and have led to deep insights into the nature of the underlying dynamics. Combining these bound-state methodologies with the formalism of inelastic neutron scattering (INS) has resulted in the novel and powerful approach for the quantum calculation of the INS spectra of a diatomic molecule in a nanocavity with an arbitrary geometry. These simulations have not only become indispensable for the interpretation and assignment of the experimental spectra but are also behind the surprising discovery of the INS selection rule for diatomics in near-spherical nanocavities. Promising directions for future research are discussed.

## I. INTRODUCTION

Ever since the serendipitous discovery of C_{60} by Kroto *et al.*,^{1} the exceptionally high symmetry and elegance of this fullerene have captured the attention and imagination of scientists. While having an unrivaled iconic status, C_{60} is only one member of the family of fullerenes, molecules consisting of an even number of carbon atoms (C_{70}, C_{72}, C_{76}, C_{80}, C_{82}, etc.) arranged to form closed cages that have pentagonal and hexagonal faces only. Regardless of the size of the carbon cage, the number of pentagons is always 12, while the number of hexagons is variable; thus, C_{60} has 20 hexagons, and C_{70} has 25. Moreover, fullerenes are subject to the so-called isolated-pentagon rule,^{2} stating that two pentagonal faces are never adjacent to each other; C_{60} is the smallest fullerene that satisfies this rule. The capability of the interior cavity of C_{60} and other fullerenes to encapsulate atoms and small molecules was noted immediately.^{1} Such supramolecular complexes are referred to as endohedral fullerenes (endofullerenes) and denoted as X@C_{n},^{3} where X is the entrapped atom or molecule and *n* is the number of carbon atoms in the fullerene cage. Endofullerenes offer a unique possibility for exploring the behavior of molecular species under the conditions unlike any other studied previously—confinement in nanoscale cages of high and variable symmetry (depending on the fullerene), that are nonpolar, homogeneous, and which to a high degree, but not entirely, isolate the guest species from the environment beyond the host cage.

The first endofullerene to be isolated actually did not involve C_{60}, but the fullerene C_{82} encapsulating an atom of lanthanum, i.e., La@C_{82}.^{3} It was followed by the preparation of a wide variety of other endohedral metallofullerenes (EMFs) where one or more metal atoms, and often nonmetals (N, P, C, etc.), are trapped inside fullerenes C_{2n}, with *n* ranging from 68 to 100.^{4,5} The synthesis of EMFs entails the vaporization of graphite rods containing metal oxides by arc discharge in a helium atmosphere. The metal and nonmetal atoms, or their clusters, are entrapped during the complex process of the fullerene cage formation. Using the same approach, or by heating fullerenes with noble gases at high pressure, it has been possible to generate another class of endofullerenes, one that has a noble-gas atom encapsulated inside the fullerene cage, e.g., He@C_{60} and Ne@C_{60}.^{6–8} The insertion of two He or two Ne atoms inside the larger fullerene C_{70} was achieved in the same way as well, resulting in the endofullerenes He_{2}@C_{70} and Ne_{2}@C_{70}.^{9,10} The yield of endofullerenes produced by these procedures is low, generally well below 1%.

There was much interest in going beyond EMFs and noble-gas endofullerenes and achieving the encapsulation of small light molecules inside the fullerene cages. But the physical processes available for production, involving extreme temperatures and high pressures, were not suited for this task. The badly needed, crucial alternative was put forth by Rubin,^{11–13} who proposed the “molecular surgery” procedure in which (i) an orifice is opened in the fullerene cage through a sequence of chemical reactions, (ii) the guest molecule is then inserted through the orifice into the cavity, under elevated temperature and pressure, and (iii) the pristine fullerene cage is reformed, leaving the guest molecule permanently trapped in its interior. This technique proved remarkably successful, leading within the past decade to the synthesis of the light-molecule endofullerenes (LMEFs) H_{2}@C_{60},^{14,15} H_{2}O@C_{60},^{16} and most recently, HF@C_{60}.^{17} The structure of H_{2}O@C_{60} is depicted in Fig. 1. Going one step further, with molecular surgery it has been possible to encapsulate two light molecules inside the larger cage of C_{70}, resulting in the extraordinary endofullerenes (H_{2})_{2}@C_{70},^{18} (H_{2}O)_{2}@C_{70},^{19} and even the mixed endofullerene (H_{2}O·HF)@C_{70}.^{20} In addition, larger light molecules such as NH_{3}^{21,22} and CH_{4}^{22,23} have been incorporated into various open-cage fullerenes derived from C_{60} by creating a larger orifice in its shell, but their closure to reform the C_{60} cage (and obtain NH_{3}@C_{60} and CH_{4}@C_{60}) has not been accomplished yet. Owing to these feats of organic synthesis, light-molecule endofullerenes, hereafter denoted as LMEFs, are now available in macroscopic quantities and with high purity, enabling a wide range of spectroscopic studies of their fundamental properties and their potential for practical applications.^{24}

A major distinction between EMFs and LMEFs pertains to the sharp difference in the dynamical behavior of the guest species in these two classes of endofullerenes. The guest species in the EMFs, individual metal atoms or their small clusters, are very heavy. As a result, their properties are discussed largely in terms of preferential position and orientation within the cage, as well as their classical motions in the cage interior.^{4,5} Dynamical quantum effects are viewed as playing a minor role at best and are seldom considered.

By contrast, the dynamics of light molecules, such as H_{2}, HF, and H_{2}O (and their respective isotopologues), confined in the interior of C_{60} and other fullerenes is dominated by strong quantum effects.^{25} This is true in particular for the low temperatures (typically ranging from 1.5 K to about 30 K) at which the spectroscopic measurements on LMEFs are usually carried out. In fact, the degree and the clarity with which the LMEFs manifest many of the fundamental, textbook principles of quantum mechanics are virtually unrivaled among realistic molecular systems. The quantum effects arise from several sources, one of which is the quantization of the translational center-of-mass (c.m.) degrees of freedom of the light molecules due to their confinement inside the fullerene (particle-in-a-box effect). Since the confinement is tight and light molecules have low molecular masses, the energy differences between the resulting translational eigenstates are large relative to *kT* (where *k* is the Boltzmann constant). The same is true for the quantized rotational states of the light molecules, owing to their large rotational constants. Finally, the confining potential of the fullerene cage couples the quantized translational and rotational motions of the guest light molecule, giving rise to a sparse translation-rotation (TR) (or “rattling”) energy level structure that, as discussed later, displays the fingerprints of the surprisingly rich and intricate quantum dynamics.

New and prominent quantum features are introduced in the TR dynamics by the phenomenon of nuclear spin isomerism, exhibited by guest molecules with (at least) two symmetrically equivalent nuclei. For H_{2} and H_{2}O, the two identical nuclei are ^{1}H; they are fermions since their nuclear spin is 1/2. In this case, the Pauli principle demands that the total wave function of the caged molecule, its spatial and spin components, must be antisymmetric with respect to the exchange of the two fermions. This requirement results in the particular entanglement of the spin and spatial quantum states, giving rise to nuclear spin isomers, *para* and *ortho*, of both H_{2} and H_{2}O, having total nuclear spins *I* = 0 and 1, respectively. For *para*-H_{2}, only even-*j* rotational levels are allowed (*j* = 0, 2, …), while *ortho*-H_{2} occupies only odd-*j* rotational levels (*j* = 1, 3, …), with *j* being the rotational quantum number of diatomic molecules. The rotational states of H_{2}O are conventionally labeled with the asymmetric top quantum numbers $jkakc$; for *para*-H_{2}O, *k*_{a} + *k*_{c} has even parity, while for *ortho*-H_{2}O, *k*_{a} + *k*_{c} has odd parity.^{26} The restrictions placed by the Pauli principle on the rotational quantum numbers of the spin isomers of H_{2} and H_{2}O make their already sparse TR level structure inside the fullerenes even sparser and push the TR dynamics deeper in the quantum regime.

LMEFs have been the subject of intense research activity, by both experimentalists and theorists. The fullerene cavities have come to be viewed as nanolaboratories in which the intriguing quantum behavior of light endohedral molecules can be probed under the conditions of high, although not complete, isolation from the environment outside the host cage. To date, the investigations have been driven primarily by the desire to characterize experimentally the fascinating quantum TR dynamics and spectroscopy of these highly unusual physical systems and to achieve quantitative theoretical understanding of their distinctive features. A review of these topics, with the emphasis on spectroscopy, is available,^{24} while the other articles published in the same issue of Phil. Trans. R. Soc. A devoted to LMEFs convey the breadth and vibrancy of the research pursued in this field across the globe.^{27} Another review, focused on the quantum TR dynamics of light molecules in C_{60} and its rigorous treatment (but pertinent to other fullerenes as well), was published recently.^{25} There has also been growing interest in extending these studies to consider the coupling between the endohedral species, in particular those with permanent electric dipoles (H_{2}O, HF), residing in neighboring fullerene cages.^{28–31} At low temperatures, the many-body correlations between the caged dipolar molecules can in principle result in the emergence of dipole-ordered phases and ferroelectricity,^{32} opening the possibilities for innovative nanoelectronics applications.

The hallmark of the research involving the LMEFs has been the exceptionally close and stimulating interplay, indeed a rarely seen synergy, between the high-level theory and spectroscopic experiments. In this perspective, I will survey the evolution and the steadily growing scope of the accurate quantum dynamical treatments of LMEFs over the past decade in three main areas:

Fully coupled multidimensional quantum calculations of the TR eigenstates of a broad range of paradigmatic LMEFs aimed at gaining a detailed understanding of their quantum TR dynamics at a fundamental level and elucidating its salient features. This line of research has its origin in the quantum five-dimensional (5D) calculations of the TR eigenstates of H

_{2}and its isotopologues in C_{60}^{33}and C_{70}.^{34}Rigorous quantum calculations of the inelastic neutron scattering (INS) spectra of H

_{2}and HD in C_{60}, including their temperature dependence, and the key role that the theory has played in the interpretation and assignment of the experimental INS spectra. These simulations rely on the recent developments in the general methodology for accurate quantum computations of the INS spectra of a diatomic molecule confined inside a nanoscale cavity.^{35–37}Explaining, and quantitatively accounting for, the symmetry breaking observed experimentally in the solid LMEFs M@C

_{60}(M = H_{2}, H_{2}O, HF) at low temperatures and its manifestations in the infrared (IR) and INS spectra of these systems.^{38,39}

I will also highlight how closely interconnected the venues (1)–(3) are, in ways that were largely not anticipated when this research program was initiated with (1). Finally, the rare pleasure of being totally surprised will be illustrated by the unexpected discovery of the selection rule in the INS spectroscopy of H_{2}^{40} and HD^{41} confined inside near-spherical nanocavities such as that of C_{60}, made possible by (1) and (2). This selection rule, the first ever to be established for the INS spectroscopy of discrete molecular compounds (and which according to the INS literature should not exist^{42–45}), was soon confirmed experimentally.^{46}

## II. QUANTUM TRANSLATION-ROTATION DYNAMICS

In the past decade, theory has been able to provide deep insight into the quantum TR dynamics of LMEFs at a fundamental level and describe it quantitatively from first principles. This has proved indispensable for the interpretation and analysis of the many experiments probing the dynamical features of LMEFs. Theory owes much of this success to the fact that the intermolecular TR motions of the guest molecule couple only weakly to both its intramolecular vibration(s) and those of the host fullerene, a consequence of several fortuitous properties of LMEFs: an order-of-magnitude difference between the frequencies of the TR and intramolecular mode(s) of the guest molecule, generally weak interactions between the entrapped molecule and the fullerene cage, and large disparity between the masses of the guest molecule and the host fullerene. As a result, both the caged molecule(s) and the host cage can be treated as internally rigid without introducing significant errors, leading to a large reduction in the dimensionality of the quantum dynamics problem that needs to be solved, moving it from the virtually intractable category to highly challenging but solvable. The problem is simplified further by the very weak interactions of the guest molecule inside a fullerene cage with the environment beyond it, the guest molecules occupying different fullerene cages and the other cages themselves. Hence, at least initially, an LMEF can be treated as an isolated system to a very good approximation, providing an excellent starting point for gradual inclusion of the weak interactions with the extended environment, the consequences of which are experimentally observable at higher resolution.

### A. Single encapsulated molecule

#### 1. Quantum calculations of fully coupled TR eigenstates

Taking the guest molecule, denoted as A, and the host fullerene to be rigid, the quantum TR dynamics of A can be described using two sets of coordinates, **R** and *ω*. **R** is the vector pointing from the center of the cage, that is also the origin of the space-fixed (SF) Cartesian axis system attached to the cage, to the center of mass (c.m.) of A. When A is a linear molecule, *ω* ≡ (*θ*, *ϕ*) are the polar and azimuthal angles, respectively, that **r**, the internuclear vector of A, makes with the SF system. For nonlinear molecules, *ω* ≡ (*ϕ*, *θ*, *χ*) are the Euler angles^{47} that specify the orientation of a body-fixed (BF) axis system affixed to A relative to the SF frame. The position vector **R** can be expressed in terms of the SF Cartesian coordinates {*x*, *y*, *z*}; these are well-suited for A inside cages of arbitrary geometries, spherical and nonspherical.^{25,33,34,48,49} An alternative representation of **R** is in the spherical polar coordinates {*R*, Ω}, where *R* ≡ |**R**|, Ω ≡ (Θ, Φ) are the polar and azimuthal angles of **R** relative to the SF axes. They are appropriate when the confining cage is nearly spherical, such as C_{60}, where they work well computationally and offer advantages in the analysis of the results.^{25,38,39,50–54}

In the rigid-monomer approximation, and assuming that the host fullerene cage is infinitely heavy and nonrotating (since it is embedded in the solid fullerene), the TR Hamiltonian of the nanoconfined A is given by

where $T^$ is the kinetic-energy operator associated with the TR motions of rigid A and *V*_{A−fullerene} is the intermolecular potential energy surface (PES) for A inside the fullerene cage (5D when A is a linear molecule and 6D for nonlinear A), discussed later in the article. When A is linear

where ∇^{2} is the Laplacian associated with **R** (in Cartesian or spherical polar coordinates), *Ĵ*^{2} is the operator corresponding to the square of the rotational angular momentum of A, *M* is the mass of A, and *I* is its moment of inertia. For nonlinear A,

Here ∇^{2} and *M* have the same meaning as in Eq. (2), $\u0134k2$ corresponds to the square of the rotational angular momentum of A about its *k*th principal axis, and *I*_{k} is the moment of inertia of A about its *k*th principal axis.

The TR energy levels and eigenfunctions of *Ĥ* in Eq. (1) are obtained by constructing its matrix representation in a suitable basis and diagonalizing it. In the early work on this topic, dealing with the quantum TR dynamics of H_{2} and isotopologues inside C_{60},^{33,34,48} and the anisotropic cages of C_{70}^{34} and an open-cage derivative of C_{60},^{49} the position vector **R** was expressed in Cartesian coordinates {*x*, *y*, *z*}, and the 3D direct-product sinc-discrete variable representation (DVR) basis^{55} was used for them. The DVR, introduced by Light and co-workers,^{56} offers several distinct advantages: first, it eliminates the need for multidimensional numerical integration of the potential matrix elements and second, it can be pruned and tailored to the features of the PES by retaining only those DVR basis functions which lie in the relevant, energetically accessible regions of the PES. In addition, the DVR enables the implementation of the sequential diagonalization and truncation procedure^{57–60} that has proved very effective in reducing drastically the size of the Hamiltonian matrix without the loss of accuracy. The spherical harmonics $|j,mj=Yjmj(\theta ,\varphi )$ have been used as the angular basis. Diagonalization of the TR Hamiltonian matrix formed in this 5D basis yields the eigenstates of A.

More recent work focused primarily on the light molecules H_{2},^{38,39,51,53,61} H_{2}O,^{38,39,52} and HF^{38,39,54} in C_{60} has utilized the eigenfunctions |*n*, *l*, *m*_{l}⟩ of the 3D isotropic Harmonic oscillator (HO) for the **R** degrees of freedom and the rigid-rotor eigenstates to cover the *ω* degrees of freedom. For A = H_{2}, HF, the latter are the usual $|j,mj=Yjmj(\omega )$, and for A = H_{2}O, they are the symmetric-top eigenfunctions of the form $|j,mj,k=2j+1/8\pi 2[Dmj,k(j)(\omega )]*$. Together, they form a direct-product 3D-HO/rigid rotor basis in 5D (linear A) or 6D (nonlinear A). An efficient method for computing the TR eigenstates of *Ĥ* in such a basis is the Chebyshev version^{62} of filter diagonalization^{63} that has been used in a number of recent studies of these systems.^{25,28,29,39,52,53}

#### 2. Intermolecular potential energy surfaces

Accurate first-principles calculation of the intermolecular rigid-monomer PES for LMEFs [*V*_{A−fullerene} in Eq. (1)] in 5D or 6D is highly challenging and has been largely an elusive goal so far. The difficulties mainly stem from two sources. One is the large size of these systems with 60, 70, or more carbon atoms, that only the density functional theory (DFT) can cope with, in principle. Moreover, the interactions between the guest molecule and the fullerene interior are dominated by dispersion, in which the standard variants of DFT fail to describe.^{64} While many approaches have been developed to treat dispersion within DFT,^{64–66} the intermolecular (5D) PES for only one LMEF, HF@C_{60}, has been calculated using DFT to date.^{54}

For the other LMEFs considered, H_{2}@C_{60},^{33,34,38–40,48,53} H_{2}@C_{70},^{34,67} and H_{2}O@C_{60},^{38,39,52} the intermolecular PESs *V*_{A−fullerene} in Eq. (1) (5D for H_{2} and 6D for H_{2}O) have been constructed as sums of pairwise-additive Lennard-Jones (LJ) potentials between the sites on the guest molecule A and each C atom of the confining fullerene

Here, *i* runs over the *m* sites on A (*m* = 3 for H_{2} and H_{2}O), *k* runs over the nuclear positions of the *n* C atoms of the fullerene, *r*_{ik} is the distance between site *i* and site *k*, while *ϵ*_{i} and *σ*_{i} are the LJ parameters for the site *i* on A. In the case of H_{2}O@C_{60},^{52} the LJ parameters were taken from the literature.^{68} Also appearing in Eq. (4) are the weights $wi$ that for H_{2}O@C_{60} are all equal to 1.^{52} However, for H_{2}@C_{60}, achieving a near-quantitative match between the calculated TR energy levels and those observed in the IR spectra of H_{2}@C_{60}^{51,61} required the introduction of an additional site 1 at the midpoint of the H–H bond^{34} (besides the sites 2 and 3 at the H nuclei, with $w2$ = $w3$ = 1). The weight $w1$ associated with site 1 and the LJ parameters *ϵ* and *σ* were optimized so as to reproduce very closely the low-lying TR levels measured in the IR spectra of H_{2}@C_{60}.^{51,61} The same functional form of the intermolecular PES was used for H_{2}@open-cage fullerene (ATOCF^{69}), although additional terms were required because of the presence of noncarbon (N, O, S, H) atoms in ATOCF.^{49}

#### 3. Near-spherical confinement: M@C_{60} (M = H_{2}, H_{2}O, HF)

The paradigmatic endofullerene H_{2}@C_{60} was the first LMEF to be synthesized.^{14,15} The quantum TR dynamics of H_{2} and its isotopologues HD and D_{2} caged in C_{60} has been explored by a large number of theoretical^{25,33,34,40,41,48,51,53} and experimental studies, primarily inelastic neutron scattering (INS)^{70–73} and infrared (IR).^{51,61,74,75} Theoretical studies^{33,34,48} were the first to reveal that the apparent simplicity of H_{2}@C_{60} is deceptive and actually conceals surprisingly rich and intricate quantum dynamics whose salient features arise from the TR coupling induced by the confining C_{60} cage with icosahedral (*I*_{h}) symmetry. Figure 2 shows the lower lying TR energy levels of *para*-H_{2} and *ortho*-H_{2} molecules inside C_{60} from the quantum 5D calculations.^{46} A conspicuous feature of this TR energy level structure is numerous closely spaced multiplets, each exhibiting a rather peculiar degeneracy;^{33,34,48} the splitting between the members of a multiplet is much smaller than the energy separation between different multiplets. As elucidated by theory,^{33,34,48} the TR level structure of H_{2}@C_{60} can be understood in terms of the model with the following ingredients: (*i*) the purely translational eigenstates of the entrapped H_{2} can be assigned with the quantum numbers of the 3D isotropic HO, suggested by the high symmetry of C_{60}; these are the principal quantum number *n* = 0, 1, 2, …, and the orbital angular momentum quantum number *l* = *n*, *n* − 2, …, 1 or 0, for odd and even *n*, respectively. (*ii*) The rotational energy levels of the caged H_{2} can be assigned in terms of the quantum numbers *j* = 0, 1, 2, …, of a linear rigid rotor.

However, (*i*) and (*ii*) alone cannot account for the fact that in Fig. 2 the TR levels with the same (nonzero) quantum numbers *n* and *j* are clumped into distinct multiplets. This prominent feature requires adding another key element to the model:^{33,34,48} In the case of the TR eigenstates which are excited both translationally and rotationally, the orbital angular momentum ** l** associated with the translational c.m. motion of H

_{2}and the rotational angular momentum

**couple vectorially to give the total angular momentum**

*j***=**

*λ***l**+

**j**having the values

*λ*=

*l*+

*j*,

*l*+

*j*− 1, …, |

*l*−

*j*|, with the degeneracy of 2

*λ*+ 1. The values of

*l*are those allowed for the quantum number

*n*. The clumps in Fig. 2 correspond to the eigenstates with the same (nonzero) quantum numbers

*n*and

*j*that, due to the TR coupling, are split into closely spaced levels, each labeled with one of the possible values of

*λ*, and exhibiting the degeneracy of 2

*λ*+ 1.

^{33,34,48}We will refer to these splittings as

*λ*splittings; they are on the order of several wavenumbers. The predicted splitting of the (nine)

*n*=

*l*= 1,

*j*= 1 TR eigenstates, associated with the translational fundamental of

*ortho*-H

_{2}, into three levels corresponding to

*λ*= 1, 2, and 0, respectively,

^{33,34,48}was subsequently observed in both the IR

^{51}and INS spectra

^{71}of H

_{2}@C

_{60}, thus validating this model of the TR coupling.

It is worth pointing out that although the translational (*j* = 0) excitations of H_{2} (and HD, D_{2}) inside C_{60} can be assigned using the quantum numbers *n* and *l* of the 3D isotropic HO, they are in fact *not* harmonic^{33,34,48} since (*i*) their energies depend not only on *n* (as they do for the 3D isotropic HO) but also on *l* and (*ii*) the energy separation between successive translational (*j* = 0) levels with *n* = 1, 2, … *increases* with *n*, evidence of their *negative* anharmonicity.

The quantum 5D calculations^{33,34,48} have also revealed that the 2*λ* + 1 degeneracy of the TR levels of H_{2}@C_{60} for *λ* > 2 is split by the “crystal field” of the icosahedral (*I*_{h}) environment of C_{60}. Here, we denote these splittings as *m*_{λ} splittings. Thus, in accordance with a group-theoretical prediction,^{76} *λ* = 3 and 4 levels are split into closely spaced pairs of levels with the degeneracies 3 and 4 for *λ* = 3, and 4 and 5 for *λ* = 4. The small magnitude of the *m*_{λ} splittings, typically less than 1 cm^{−1}, reflects the weakness of the *I*_{h} corrugation, i.e., the small deviation from spherical symmetry of the interaction potential for H_{2} in C_{60}.

It is of considerable interest to find out whether the magnitudes of the *λ* splittings and *m*_{λ} splittings, as well as the energy ordering of the *λ* components in given (*n*, *l*, *j*) multiplet, can be connected with any particular part(s) of the TR Hamiltonian. For this purpose, borrowing from the work of Olthof *et al.*^{50} and Mamone *et al.*,^{51} the 5D intermolecular PES *V*(*R*, Ω, *ω*) of H_{2}@C_{60} in Eq. (1), where *R* ≡ |**R**|, Ω ≡ (Θ, Φ), and *ω* ≡ (*θ*, *ϕ*), is decomposed into a sum over bipolar spherical tensors,^{25,53}

where the $GL,JK,Q(R)$ are scalar functions of *R* and the $TL,J(K,Q)(\Omega ,\omega )$ are the bipolar spherical tensors. Variational quantum 5D calculations of the TR eigenstates of H_{2}@C_{60} utilizing this expansion of the intermolecular PES^{53} have established that the *λ* splittings are essentially determined by just the single term, $T2,2(0,0)$, in Eq. (5), consistent with the work Mamone *et al.*^{51} This has prompted going one step further and using the $T2,2(0,0)$ term as the perturbation in the first-order perturbation theory (PT) treatment;^{53} in general, this PT approach accounts quantitatively for both the magnitudes and the patterns of the *λ* splittings in H_{2}@C_{60}. Moreover, it turns out that the analysis of the $T2,2(0,0)$ term provides a remarkably simple physical explanation for the energy ordering of the *λ* components within a given (*n*, *l*, *j*) multiplet that has been a long-standing puzzle.^{53}

Are the salient features of the TR level structure of H_{2}@C_{60} unique to this species, or are they common to other light-molecule endohedral complexes involving C_{60}? Recent computational results point to the latter. The TR energy level structure of H_{2}O@C_{60}, from fully coupled quantum 6D calculations,^{52} exhibits all the above qualitative features of H_{2}@C_{60}. These include the use of 3D HO quantum numbers *n* and *l* to assign the translational excitations of H_{2}O, while for the H_{2}O rotational excitations the asymmetric top quantum numbers $jkakc$ are employed. The orbital angular momentum of the c.m. of H_{2}O couples to the rotational angular momentum of the molecule, as in H_{2}@C_{60}, giving rise to the characteristic *λ* multiplets. In addition, the TR levels of H_{2}O@C_{60} with *λ* ≥ 3 manifest the same fine *I*_{h} “crystal field” splitting.^{52} The TR eigenstates of HF@C_{60} from the recent quantum 5D calculations^{39,54} are also assignable in terms of the quantum numbers (*n*, *l*, *j*, *λ*) utilized for H_{2}@C_{60}^{33,48} and HD@C_{60}.^{48} One can conclude that the quantum TR dynamics of the three light-molecule endofullerenes M@C_{60} (M = H_{2}, H_{2}O, HF) investigated so far can be described within a unified theoretical framework. It should be emphasized here that these calculations have assumed that M is in the confining environment with *I*_{h} symmetry, that of an isolated C_{60} cage. But, as discussed in Sec. IV, this assumption is not strictly valid for solid LMEFs at low temperatures, which introduces additional important experimentally observable features.

Besides the physically motivated quantum numbers (*n*, *l*, *j*, *λ*), the TR eigenstates of H_{2}@C_{60} can each be given a symmetry assignment in terms of one of the irreducible representations (irreps) of the symmetry group $Ih(12)=Ih\u2297S2$.^{77} The same holds for the TR energy levels of H_{2}O@C_{60}.^{52} The group-theory-based assignments are rigorous and account fully for the degeneracies of the computed TR levels. On the other hand, (*n*, *l*, *j*, *λ*) assignments provide the essential physical insight and information that the symmetry-based labels lack entirely: they specify the nature and the degree of the excitations for each TR level and account for the clumps, or multiplets, that are ubiquitous in the TR energy level structure. Therefore, the two types of assignments are truly complementary.

Finally, how good are the quantum numbers *n*, *l*, *j* ($jkakc$ for H_{2}O), and *λ*? Recent calculations^{39} have demonstrated that the total angular momentum quantum number *λ* is close to an integer, and therefore a good quantum number, for the lower lying TR eigenstates of M@C_{60} (M = H_{2}, H_{2}O, HF). The situation is somewhat less uniform with respect to the quantum numbers *n*, *l*, *j* for H_{2} and HF, and $n,l,jkakc$ for H_{2}O. In the case of H_{2} and H_{2}O, for the lower energy TR eigenstates considered, the contribution of the dominant direct-product 3D-HO/rigid rotor basis function |*n*, *l*, *m*_{l}, *j*, *m*_{j}⟩ (for H_{2}) or $|n,l,ml,jkakc,mj$ (for H_{2}O), described earlier in the paper, is close to 1. The implication is that in this range of excitation energies, the quantum numbers *n*, *l*, *j* are good for H_{2} and the quantum numbers $n,l,jkakc$ are good for H_{2}O.^{39} By contrast, for HF, apart from the few lowest eigenstates, the contribution of the dominant basis function |*n*, *l*, *m*_{l}, *j*, *m*_{j}⟩ is typically in the range 0.5–0.6, implying that *n*, *l*, *j* are not nearly as good quantum numbers for HF as they are for H_{2}. This has been attributed to the combined effects of the large mass asymmetry of HF and the fact that its translational fundamental is close in energy to the lowest-lying rotational excitation of the species.^{39}

#### 4. Reduced-symmetry confinement: H_{2} inside C_{70} and an open-cage fullerene

The symmetry of the confining nanocage has a strong influence on the quantum TR dynamics of the guest molecule. It ultimately determines the patterns of the TR energy levels and the translational quantum numbers appropriate for their assignment, the splittings of the rotational excitations, and the overall TR coupling scheme. Thus, in the case of H_{2}@C_{60}, the use of the quantum numbers of the 3D isotropic HO to assign the translational excitations and the particular TR coupling responsible for the distinctive *λ* multiplets are both consequences of the high, near-spherical (*I*_{h}) symmetry of C_{60}. What happens when the symmetry of the confining fullerene cage is reduced? This question was the primary motivation for the theoretical studies of the TR level structure of H_{2} in C_{70} (*D*_{5h})^{34} and an open-cage fullerene (ATOCF^{69}) that has no symmetry,^{49} discussed in the following.

In the case of H_{2} in C_{70}, the symmetry of the cage (*D*_{5h}) is appreciably lower than that of C_{60} (*I*_{h}). C_{70} has a long axis (*z*) coinciding with its *C*_{5} axis of rotation that is distinct from the two symmetrically equivalent short axes (*x* and *y*, respectively) perpendicular to the *C*_{5} axis. The profile of the intermolecular PES along the long axis, highly anharmonic and almost flat over an extended central region of C_{70}, differs qualitatively, in shape and spatial extent, from the profiles along the two short axes. Hence, the 5D intermolecular PES of H_{2}@C_{70} is highly anisotropic with respect to the direction, *z* vs. *x* and *y*, in which the c.m. of H_{2} moves away from the center of the cage, in contrast to the PES of H_{2}@C_{60} whose directional, or radial, anisotropy is very weak. Quantum 5D calculations^{34} have revealed very weak coupling of the quasi-1D translational mode along the *z* axis to the excitations in the plane of the short *x* and *y* axes. The *z*-mode excitations are readily assignable in terms of the Cartesian quantum number $vz$, and translational excitations in the *xy* plane can be assigned using the quantum numbers $v$ and *l* of the 2D isotropic HO, where $v$ denotes the number of quanta and *l* is the vibrational angular momentum along the *z* axis; for a given $v$, *l* can take $v$ + 1 values, −$v$, −$v$ + 2, …, $v$ − 4, $v$ − 2, $v$.^{78}

The frequency of the *z*-mode fundamental, 54.15 cm^{−1}, is significantly lower than that of the 2D *xy*-mode fundamental, 138.74 cm^{−1},^{34} reflecting the greater stiffness of the intermolecular PES in the *x* and *y* directions than along the *z* axis. This is in good agreement with the corresponding experimental results for the 1D and 2D translational modes, 56 and 151 cm^{−1}, respectively, from the IR spectra of H_{2}@C_{70}.^{75} Moreover, all translational modes are characterized by negative anharmonicity. Concerning the rotational excitations of H_{2} in C_{70}, the calculations predict the 1:2 splitting (by 7.4 cm^{−1}) of the triply degenerate *j* = 1 ground state of *ortho*-H_{2} into a non-degenerate lower energy state and a doubly degenerate higher energy state.^{34} The calculated magnitude of the splitting agrees remarkably well with the experimental value of about 7 cm^{−1} from IR^{75} and nuclear magnetic resonance (NMR) spectroscopy^{79} of H_{2}@C_{70}. However, both NMR and IR results support the 2:1 splitting pattern of the ground state of *ortho*-H_{2}, opposite from the calculated one. It is very likely that a more refined description of the H_{2}-nanocarbon interaction is needed to reproduce this subtle difference.

In the endohedral complex of H_{2} with the aza-thia-open-cage fullerene (ATOCF),^{69} denoted as H_{2}@ATOCF, the (open) cage has no symmetry elements. Therefore, quantum 5D calculations of the TR eigenstates of H_{2}@ATOCF^{49} represent the next logical step in investigating how systematically decreasing the cage symmetry affects the salient features of the quantum TR dynamics of the encapsulated H_{2}. Due to the lack of symmetry in ATOCF, the profiles of the intermolecular PES along the three principal axes *x*, *y*, and *z* (that points in the direction of the cage orifice) differ substantially, with the cut along *z* being the most anharmonic. The plots of the 3D reduced-probability densities of the computed eigenstates reveal three well-defined quasi-1D translational modes, one along each of the principal axes that can be assigned in terms of the Cartesian quantum numbers ($vx$, $vy$, $vz$). The strong radial anisotropy of the intermolecular PES completely splits the translational fundamental into three excitations with very different frequencies; that of the *z* mode, toward the cage opening, is the lowest. In the same vein, the three-fold degeneracy of the *j* = 1 ground state of *ortho*-H_{2} in ATOCF is lifted entirely by the angular anisotropy of the PES.^{49} The agreement between the theoretical results and the INS measurements^{80} is semiquantitative.

### B. Two encapsulated molecules: (H_{2})_{2}@C_{70}

The molecular surgery approach^{11–13} has enabled the encapsulation of two light molecules inside the cage of C_{70}, yielding the endofullerenes (H_{2})_{2}@open-cage C_{70},^{81} (H_{2})_{2}@C_{70},^{18} (H_{2}O)_{2}@C_{70},^{19} and (H_{2}O·HF)@C_{70}.^{20} Such endohedral complexes present a rare opportunity for investigating various aspects of the quantum TR dynamics of two molecules confined in the nanocavity of rather high symmetry and well defined geometry. However, a rigorous theoretical treatment of this problem presents major challenges for theory, stemming primarily from the high dimensionality of the calculations involved. In the rigid-monomer approximation, each linear encapsulated molecule adds five (coupled) degrees of freedom that must be included and six if the molecule is nonlinear. Thus, (H_{2})_{2}@C_{70} presents a 10-dimensional quantum problem that requires a fully coupled approach; for (H_{2}O)_{2}@C_{70} the quantum calculations would be 12-dimensional. Both cases are close to the limits of what is currently feasible. To date, only two such studies have been reported,^{67,82} both for (H_{2})_{2}@C_{70}, and are reviewed below.

#### 1. Ground-state properties: Quantum effects in the energetics, maximum H_{2} occupancy, and spatial distribution

Diffusion Monte Carlo (DMC) calculations of the quantum TR dynamics have been performed for up to three *para*-H_{2} molecules encapsulated in C_{70} and, for comparison purposes, for one and two *para*-H_{2} molecules inside C_{60}.^{67} Their objective was to elucidate the role of quantum effects in the ground-state properties of these endohedral complexes: their energetics and the number of H_{2} molecules that can be stabilized in the C_{70} and C_{60} cages, the role that the TR zero-point energies (ZPEs) play in both, and for C_{70}, the spatial distributions of one and two confined *para*-H_{2} molecules. The energy of the global minimum on the intermolecular PES is negative for one and two H_{2} molecules in C_{70} but has a high positive value when the third H_{2} is added, implying that at most two H_{2} molecules can be stabilized inside C_{70}. By the same criterion, in the case of C_{60}, only the endohedral complex with one H_{2} molecule is energetically stable. Both results are in accord with the fact that (H_{2})_{n}@C_{70} (*n* = 1, 2)^{18} and H_{2}@C_{60}^{14,15} have been synthesized, but not (H_{2})_{3}@C_{70} or (H_{2})_{2}@C_{60}.

Inclusion of quantum effects via the DMC calculations leaves unchanged the conclusions regarding the at most double H_{2} occupancy of C_{70} and the single H_{2} occupancy of C_{60}. However, the ZPE of the coupled TR motions does play a major role in the energetics since it grows rapidly with the number of encapsulated *para*-H_{2} molecules and amounts to a substantial fraction of the well depth of the intermolecular PES. In the case of (*para*-H_{2})_{n}@C_{70}, the ZPE represents 11% of the global minimum for *n* = 1 and 52% for *n* = 2. For H_{2}@C_{60}, the ZPE is 16% of the potential well depth. Moreover, the inclusion of the ZPE increases by almost a factor of two the energetic penalty for adding the third *para*-H_{2} molecule in C_{70}, and by 66% for the second *para*-H_{2} in C_{60}. Consequently, accurate calculation of the TR ZPE is essential for reliable theoretical predictions regarding the stability of the endohedral fullerene complexes with hydrogen molecules and their maximum H_{2} content.^{67}

Spatial distributions of one and two H_{2} molecules in C_{70} are characterized by two probability distribution functions (PDFs) from the DMC calculations: the cage center to *para*-H_{2} c.m. distance (*R*) PDF *P*(*R*) and the *para*-H_{2}–*para*-H_{2} c.m. distance (*r*) PDF *P*(*r*). Both PDFs are shown in Fig. 3. Particularly illuminating is *P*(*R*). In going from *n* = 1 to *n* = 2, the peak of *P*(*R*) shifts from 0.9 to 2.3 a.u. and its width decreases markedly. This narrowing of *P*(*R*) shows that two H_{2} molecules are much more restricted in their motions inside the cage than one. Consistent with this, the fact that *P*(*R*) is close to zero for *R* < 1.5 a.u. reveals that the two guest molecules are virtually excluded from the central region of the cage, due to mutual repulsion, and are confined to a narrow range of *R* values.^{67}

#### 2. Excited translation-rotation eigenstates

Excited TR eigenstates of (H_{2})_{2}@C_{70} have been calculated by a nuclear-orbital/configuration-interaction method^{82} that is similar in spirit to the self-consistent field approach to molecular vibrations.^{83} In this method, one first computes iteratively a set of single-H_{2} “nuclear orbitals” which approximate the TR states of that moiety interacting with C_{70} and the mean field of the second H_{2} molecule. Next, configurations that approximate the TR states of (H_{2})_{2}@C_{70} are constructed as bilinear products of these orbitals populated by the two H_{2} moieties, respectively. These configurations serve as the basis in which the matrix of the full TR Hamiltonian of the system in the rigid-monomer approximation is computed and diagonalized, yielding the TR eigenstates of the species. This approach has a number of attractive features, including the observation that single configuration basis states generally constitute very good zeroth-order approximations to the final TR eigenstates. This reduces the size of the configuration basis required to compute well-converged TR eigenstates and is very helpful in their interpretation. Analysis of the results shows that the low-energy TR eigenstates of (H_{2})_{2}@C_{70} can be described in terms of a small number of basic H_{2} excitation types.^{82}

## III. RIGOROUS QUANTUM TREATMENT OF INELASTIC NEUTRON SCATTERING SPECTRA

INS spectroscopy is an extremely valuable and versatile tool for probing the quantum TR dynamics of molecules with one or more hydrogen atoms, e.g., H_{2}, HD, H_{2}O, and CH_{4}, encapsulated in the nanoscale cavities of a wide variety of host materials, such as fullerenes, clathrate hydrates, metal-organic frameworks, and zeolites, among others.^{42,43} Its power and broad applicability stem from two unique features. One of them is the exceptionally large cross section for the incoherent neutron scattering from the hydrogen (^{1}H) nucleus,^{84} ∼15 times greater than for any other nucleus, including the isotope deuterium (^{2}H). This makes INS an exquisitely selective probe of the dynamics of entrapped H_{2}, HD, and other H-containing molecules. The second distinctive feature of the INS is the ability of neutrons to change the total nuclear spin *I* of the guest molecule, something that photons cannot do. As a result, rotational Δ*j* = 1 transitions can be observed that are forbidden in the optical, IR, and Raman spectroscopy, such as *j* = 0 → 1 of *para*-H_{2} (*I* = 0) and *j* = 1 → 0 and *j* = 1 → 2 of *ortho*-H_{2} (*I* = 1) that interconvert *para*- and *ortho*-H_{2}.

In the context of LMEFs, INS spectroscopy has been used to investigate the TR dynamics of H_{2}@C_{60},^{70–72} H_{2}@ATOCF,^{80} H_{2}O@C_{60},^{85,86} and HF@C_{60}.^{17} The measured INS spectra encode a wealth of information about the dynamical behavior of the guest molecule and its anisotropic interactions with the host cage. However, given the complexity of the TR dynamics, complete and quantitative extraction of this rich information content is possible only with the help of high-level theory. The experimental spectra can be partly analyzed by aligning the TR excitation energies from quantum 5D (or 6D) calculations with the measured INS transitions. But at higher energies, as the spectrum becomes denser, spectral assignment based solely on the comparison of calculated and observed TR excitations becomes increasingly unreliable and ambiguous. The analysis and assignment of the experimental INS spectra of LMEFs would be much more complete and reliable if in addition to accurate TR excitation energies, the theory was also capable of yielding the intensities of the INS transitions. This would generate the complete spectroscopic fingerprint of the system considered for direct and more conclusive comparison with the experimental INS data. However, the methodology for such calculations was not available.

Rather recently, this gap began to close with the development of the approach that enables accurate and efficient quantum calculation of the INS spectra, i.e., the energies and the INS intensities of the TR transitions, of an H_{2} molecule confined in a nanocavity of an arbitrary shape, assuming the rigid-monomer approximation.^{35,36} The starting point of our quantum treatment is the standard basic equations of the INS theory.^{87,88} What sets it apart from all other treatments in the literature is the incorporation of the quantum 5D TR wave functions of the entrapped hydrogen molecule, calculated by the methods described in Sec. II A 1, as the spatial parts of the initial and final states of the INS transitions. The INS spectra simulated using this approach exhibit an unprecedented degree of realism and reflect in full the complexity of the TR dynamics of the guest molecule in nanoconfinement. The power of the new methodology was demonstrated first by calculating the INS spectra of *para*- and *ortho*-H_{2} in the small cage of the sII clathrate hydrate,^{35} in remarkably good agreement with the experimental data.^{89} Subsequently, this approach was extended to general heteronuclear diatomic molecules entrapped inside a nanoscale cavity.^{37} It was used to compute the INS spectra of HD molecule in the small cage of the structure II clathrate hydrate that agrees very well with the measured INS spectra of this system.^{89,90}

Another, a conceptually different approach to the calculation of the INS spectra of homonuclear diatomic molecules confined within spherical nanocavities has been formulated and implemented on H_{2}@C_{60}.^{73} In this methodology, the intermolecular PES is not calculated, or constructed, *a priori*. Instead, it is defined in terms of a few physical parameters, such as the harmonic term, anharmonic corrections, and TR coupling, that are determined by matching the simulations to the experimental spectrum.

In the following, the quantum methodology for computing the INS spectra of diatomics inside nanocavities is outlined, together with the application to H_{2}@C_{60}.^{40} Also covered is the totally unexpected discovery of the very first selection rule in the INS spectroscopy of discrete molecular compounds,^{40} that was validated soon thereafter,^{46} and generalized.^{77,91}

### A. Methodology and application to H_{2}@C_{60}

For the case of an entrapped homonuclear molecule such as H_{2}, the details of the formalism are given in Ref. 36; its extension to the caged heteronuclear diatomics, e.g., HD, is presented in Ref. 37. Here, we only define the key quantities, give the most essential equations, and elaborate on some important aspects of our quantum simulation of the INS spectra.

Let *E* and *E*′ be the energies of the incident and scattered neutron, respectively, and $k\u2192$ and $k\u2192\u2032$ be their respective wave vectors. Then the neutron energy transfer can be written as *ℏω* = Δ*E* = *E* − *E*′, and the neutron momentum transfer can be written as $\u210f\kappa \u2192=\u210fk\u2192\u2212\u210fk\u2192\u2032$. The double differential cross section for neutron scattering in the first Born approximation is^{87,88}

where

and

In Eqs. (6)–(8), $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}, $b^n$ is the scattering length operator (Ref. 36), and $r\u2192n$ is the position of nucleus *n*. In our case, the spatial components of |*i*⟩ and |*f*⟩ are different TR states of H_{2} inside C_{60} and are represented by the 5D TR eigenfunctions of the Hamiltonian in Eq. (1). The evaluation of $S(\kappa \u2192,\omega )$ in Eq. (7) using the quantum 5D TR eigenstates is described in Ref. 36.

It should be noted that the expressions in Eqs. (6)–(8) involve summation over the initial TR eigenstates of the system weighted by their relative populations according to the Boltzmann distribution for the given temperature. Therefore, one can explore computationally the temperature dependence of the INS spectra, which was done for H_{2}@C_{60}^{40} and HD@C_{60}.^{41}

The experimental INS spectra of H_{2}@C_{60} and other LMEFs are taken from powdered samples, where the fullerene cages are randomly oriented with respect to the incoming neutron beam. Consequently, for a more realistic comparison with the measured spectra, the computed INS spectra are averaged over all possible orientations of C_{60}.^{36} Moreover, if the INS spectrometer used in the measurements, e.g., IN4C,^{17,71} allows the spectra to be recorded for a broad range of $\kappa \u2261|\kappa \u2192|$ values, then *S*(*κ*, *ω*), obtained by averaging $S(\kappa \u2192,\omega )$ in Eq. (7) over all directions of the wave vector $\kappa \u2192$,^{36} needs to be integrated also over *κ* to obtain the desired *S*(*ω*).^{41}

The performance and the utility of this methodology are demonstrated in Fig. 4. Shown there is the INS spectrum of H_{2}@C_{60} measured at 2.5 K^{46} that greatly extends the range of energy transfer recorded previously^{71} from 60 meV to beyond 250 meV. Superimposed on the measured INS spectrum is the spectrum of H_{2}@C_{60} calculated for 2.5 K^{46} and the statistical 3:1 mixture of *ortho*- and *para*-H_{2}, expected for the room temperature (or higher) at which H_{2}@C_{60} is prepared.^{92} Subsequent cooling of the sample in the absence of an external catalyst does not alter this equilibrium ratio. At 2.5 K, all INS transitions, observed or calculated, must originate in the ground TR states of either *ortho*- or *para*-H_{2} since they are the only ones populated at this low temperature. Already a glance at Fig. 4 allows one to make two observations: (*i*) the measured and calculated INS spectra are in remarkable agreement over the entire energy range shown, a testimony to the accuracy of the computational methodology and the high quality of the PES employed. This permits the interpretation and assignment of the entire recorded INS spectrum, as was done previously^{40} for the INS transitions up to ∼60 meV. (*ii*) The majority of observed bands corresponds to a number of unresolved individual TR transitions. This implies that in the absence of the accurately calculated INS spectrum, reliable assignment of the experimental spectrum over a substantial energy range would be virtually impossible.

While theory is crucial for the interpretation of even the low-temperature INS spectra at higher excitation energies, the studies of H_{2}@C_{60}^{40} and HD@C_{60}^{41} have shown that it becomes even more indispensable for the analysis of the INS spectra taken at higher temperatures. Then, a growing number of excited TR states becomes appreciably populated, greatly increasing the number of transitions capable of giving rise to new peaks in the INS spectrum that are often blended into broad bands, thus making it much more congested.

### B. Selection rule that should not exist

From the calculations of the INS spectra of H_{2}@C_{60},^{40} it emerged that transitions out of the ground state of *para*-H_{2} to certain excited TR states have intensities that are *zero* to at least 4-5 significant figures. Zero-intensity transitions are present in the calculated INS spectra of HD@C_{60} as well.^{41} This strongly suggested that such transitions are in fact forbidden according to some as yet unknown INS selection rule. But this conjecture ran counter to the widely held view, expressed in numerous reviews^{42–44} and textbooks^{45} on the subject, that there are *no* selection rules in the INS spectroscopy of discrete molecular compounds; all vibrations are supposed to be active in INS and, in principle, measurable.^{45} Nevertheless, the computational evidence for the existence of a selection rule seemed compelling and we undertook to derive one pertaining to *para*-H_{2} and HD inside a near-spherical nanocavity.^{40}

The derivation starts from the spatial component of the INS transition matrix element,^{35,36}

$|\Psi i5D$ and $|\Psi f5D$ in Eq. (9) are the 5D eigenstates of the TR Hamiltonian in Eq. (1), representing the spatial parts of the initial (|*i*⟩) and final states (|*f*⟩) of the INS transition, respectively. $R\u2192cm$ is the position vector of the c.m. of H_{2}, $\rho \u2192$ is the vector connecting the two H atoms, and $\kappa \u2192$ is the wave vector of the neutron momentum transfer.^{35,36} One needs to establish when $Ifi$ in Eq. (9) is equal to zero.

The TR eigenstates $|\Psi \tau 5D\u2009(\tau =i,f)$ can be written as

where |*n*, *l*⟩ is the radial (*R*) part of the TR eigenstate. In a key step of the derivation, since the nanocavity is assumed to be near-spherical, such as that of C_{60}, |*n*, *l*⟩ is represented by $Rnl(R)$, the radial part of the eigenfunction of the 3D isotropic HO. Then, after a rather lengthy derivation, the following INS selection rule is obtained:^{40}

For a *para*-H_{2} or HD molecule confined in a near-spherical nanocavity, the INS transitions between its ground TR state and excited TR states assigned as (*n*, *l*, *j*, *λ*) are *forbidden* for *λ* = *j* + *l* − 1.

This is the first ever selection rule established in the INS spectroscopy of aperiodic, discrete molecular compounds. All transitions that in the calculated INS spectra of H_{2}@C_{60}^{40} and HD@C_{60}^{41} have zero intensities are forbidden according to this selection rule. Within a year of its publication, the validity of this precedent-setting selection rule was confirmed by a joint experimental and theoretical INS investigation of H_{2}@C_{60};^{46} several transitions predicted to be forbidden by the selection rule are absent from the measured INS spectrum.

This selection rule covers only the INS transitions in which one of the two states involved is the ground TR state of *para*-H_{2} or HD. It says nothing about the transitions from the TR states of *ortho*-H_{2} or the excited TR states of *para*-H_{2} and HD. Shortly afterwards, we managed to remove this limitation by deriving the generalized selection rule, that defines *all* possible forbidden INS transitions, originating in a variety of TR states, ground *and* excited, of *ortho*-H_{2}, *para*-H_{2}, and HD, inside a near-spherical nanocavity.^{91} The INS selection rule derived previously^{40} is a special case of this generalized version.

Relying solely on group-theoretical considerations, Poirier has almost simultaneously derived the same generalized selection rule and showed that it is rigorous, or strict in the group-theory sense,^{26} for the case of H_{2} in a spherical cavity, but not in an *I*_{h} environment.^{77} This insight is valuable. On the other hand, derivations that do not rely on group-theory arguments,^{40,91} but of course invoke the near-spherical character of the nanocavity, result in the selection rules whose somewhat greater scope is apparent from the outset, such as the applicability to the INS transitions of the caged HD molecule.

## IV. SYMMETRY BREAKING IN SOLID LIGHT-MOLECULE ENDOFULLERENES

In the discussion so far, it has been generally assumed that the confining environment has *I*_{h} symmetry of the isolated *C*_{60} molecule. In such an environment, the *j* = 1 ground state of *ortho*-H_{2} would maintain the three-fold degeneracy that it has in the gas phase. However, beginning in 2009, three lines of investigation of solid H_{2}@C_{60} at low temperatures, involving (*i*) the specific heat anomaly observed around 0.6 K for H_{2}@C_{60}(*s*),^{93} (*ii*) the temperature dependence of the line corresponding to the INS transition from the ground state of *ortho*-H_{2} to the ground state of *para*-H_{2}, in the range 0.06–35 K,^{94} and (*iii*) low-temperature (6 K) near-IR (NIR) spectroscopy of H_{2}@C_{60}(*s*),^{51,75} revealed that the *j* = 1 triplet of *ortho*-H_{2} is split by about 1 cm^{−1} into a lower energy non-degenerate level and a higher energy doubly degenerate level. These findings imply that in the low-temperature H_{2}@C_{60}(*s*), the symmetry of the environment felt by the caged H_{2} must be lower than *I*_{h}.

Symmetry breaking has also been observed in the dipolar endofullerenes H_{2}O@C_{60}(*s*) and HF@C_{60}(*s*) at low temperatures. The INS measurements of H_{2}O@C_{60}(*s*)^{85,86,95} at 1.6 K showed that the three-fold degenerate 1_{01} ground state of *ortho*-H_{2}O is split by 4.19 cm^{−1} into a non-degenerate lower energy state and a doubly degenerate higher energy state. For HF@C_{60}(*s*), in the far- and mid-IR spectra measured at 5 K, the symmetry breaking manifested as the splitting of the *j* = 0 → *j* = 1 transition by 3.9 cm^{−1}.^{17}

The fact that the three LMEF solids investigated so far, M@C_{60}(*s*) (M = H_{2}, H_{2}O, HF), all display manifestations of symmetry breaking, strongly suggests that this feature is not accidental but intrinsic to such systems. Naturally, this raises the question about the origin, or mechanism responsible for the symmetry breaking. Having a permanent electric dipole moment is clearly not a prerequisite since symmetry breaking is present regardless of whether the endohedral molecule is dipolar (H_{2}O, HF) or not (H_{2}). One of the possibilities proposed was the distortion of the host C_{60} cage away from the *I*_{h} symmetry, due to the interaction with the guest molecule as it moves in its interior.^{17,85,86,94} The problem is that no quantitative predictions have been made based on this and other postulated mechanisms, making it impossible to assess their viability by comparison with the measured splittings.

Faced with this unsatisfactory state of affairs, and desiring to gain quantitative understanding of the symmetry breaking in the LMEFs M@C_{60}(*s*) (M = H_{2}, H_{2}O, HF) at low temperatures, we have recently developed the approach^{38} outlined below that proved highly successful, based on the following considerations. At temperatures below 90 K, the fullerene molecules in the solid C_{60} are locked in, and coexist as, two orientationally ordered configurations of neighboring units.^{96} In the dominant one, referred to as the P orientation, associated with the global minimum of the inter-cage potential, the electron-rich double bonds shared by two hexagons (denoted as 6:6) of one C_{60} unit face directly the electron-poor pentagons of the neighboring cages.^{96,97} Each C_{60} unit has six electron-rich 6:6 bonds and six electron-deficient pentagons facing its twelve nearest neighbors. The slightly higher local minimum on the inter-cage potential corresponds to what is denoted as the H orientation, where electron-rich 6:6 bonds of one C_{60} are immediately adjacent to the electron-poor hexagonal faces of the neighboring units.^{96} Below 90 K, and at ambient pressure, the relative proportion of molecules in the P and H orientations is 5:1.^{96}

The P and H orientations are relevant for symmetry breaking since for both the point-group symmetry of the environment at the center of a C_{60} cage is *S*_{6},^{93,94} capable of splitting the three-fold degeneracy of the ground states of *ortho*-H_{2}, *ortho*-H_{2}O, and the *j* = 1 state of HF into a non-degenerate and a doubly degenerate state, as observed experimentally. But for the consequences of this symmetry lowering to be observable, the guest molecule inside one C_{60} must have appreciable interaction with the neighboring C_{60} units, not just its own cage. The nature of the symmetry-breaking interaction has been unresolved.

In the recent paper,^{38} the origin of this crucial interaction in the endofullerenes M@C_{60}(*s*) (M = H_{2}, HF, H_{2}O) was identified and verified by comparing its predictions with the experimental results. This was accomplished by focusing on the smallest fragment of solid C_{60} that can give rise to *S*_{6} symmetry in its center, comprised of twelve nearest-neighbor (NN) cages around the central cage, with all the cages being in either the P or the H orientation (see Fig. 5). The cages are treated as rigid and assumed to have *I*_{h} symmetry. Only the central cage is occupied by the guest molecule M, while the twelve NN cages are left empty. At the heart of this approach is the proposition that the predominant source of the symmetry breaking in M@C_{60}(*s*) at low temperatures is the electrostatic interaction between the charge densities on the NN C_{60} cages and that on M in the central cage. The electron density on the NN cages is obtained from a first-principles DFT calculation, and the charge density on M is related to its body-fixed (BF) quadrupole moments reported in the literature. Thus, the M-NN electrostatic interaction contains no adjustable parameters.

For this electrostatic-interaction model, the TR Hamiltonian of M confined inside the central cage of the thirteen-cage fragment of M@C_{60}(*s*) can be written as

The first two terms in Eq. (11) correspond to the TR Hamiltonian of M inside the central C_{60} cage, given in Eq. (1) The third term, *V*_{ES}, is the electrostatic interaction between the charge density distribution on M and that on the twelve nearest-neighbor (NN) C_{60} cages that surround the central cage. Only the leading term in the multipole expansion of the M-NN interaction is retained in *V*_{ES}.^{38} Owing to the *S*_{6} point-group symmetry of the 13-cage (central + NN) crystal fragment,^{98,99} that term is the quadrupolar one

where the $Qm(2)$ constitute the electric-quadrupole spherical tensor of M@C_{60} and the $Im(2)$ constitute the electric-field-gradient tensor arising from the charges on the NN cages. The key proposition of the model is that retaining just the leading quadrupole terms in *V*_{quad} suffices to quantitatively account for the manifestations of symmetry breaking observed in all three M@C_{60}(*s*) (M = H_{2}, HF, H_{2}O) species considered.^{38}

The splittings of the ground states of *ortho*-H_{2}, *ortho*-H_{2}O and the *j* = 1 level of HF inside the central C_{60} and under the quadrupolar interaction with the twelve NN cages have been computed for the P and H orientations by means of variational calculations, as well as by the first-order perturbation theory.^{38} For the P orientation, the calculated splittings agree extremely well, to within the experimental error bars, with those measured for the corresponding M@C_{60}(*s*) systems. In addition, the calculations yield the 1:2 splitting patterns in agreement with the experiments. This quantitative success of the model provides a strong argument in favor of its central tenet that the electrostatic M-NN interaction, dominated by the leading quadrupole terms, is the main mechanism for the low-temperature symmetry breaking in the M@C_{60}(*s*) species considered. Surprisingly, the symmetry-induced splittings calculated for the H orientation are about a factor of 30 smaller than the corresponding ones in the P orientation, for the three endofullerene systems considered.^{38}

In an extension of this study, we have investigated the effects of the symmetry breaking on the excited TR eigenstates of the endofullerenes M@C_{60}(*s*) (M = H_{2}, HF, H_{2}O), as well as its manifestations in the simulated near-IR (NIR) and far-IR (FIR) spectra of these species.^{39} Utilizing the computed TR eigenstates, the low-temperature NIR (for M = H_{2}) and FIR (for M = HF, H_{2}O) spectra were simulated for both the P and H orientations and also combined as their weighted sum, 0.15H + 0.85P (corresponding to the measured abundances of the two crystal types). The weighted-sum spectra computed for M = H_{2} and HF match quantitatively the corresponding measured spectra, while for M = H_{2}O the weighted-sum FIR spectrum predicts features that can potentially be observed experimentally.^{39}

## V. CONCLUSIONS AND FUTURE DIRECTIONS

The theoretical advances made over the past ten years have resulted in deep understanding of the fundamental aspects of the TR dynamics of light molecules encapsulated inside fullerenes and other types of nanocavities, as well as the ability to describe the TR level structure accurately, from first principles. This includes the recent quantitative account of the symmetry breaking observed in the solid light-molecule endofullerenes at low temperatures that eluded explanation for a decade. Combining accurately calculated TR eigenstates of nanoconfined molecules with the INS formalism has resulted in the novel approach for the quantum simulation of the INS spectra of diatomic molecules inside an arbitrarily shaped nanocavity, that exhibit a high degree of realism and accuracy. Besides becoming an indispensable partner in the analysis of the experimental INS spectra, the new methodology has led to some genuine conceptual advances, including the surprising discovery of the INS selection rule, in the field that none was supposed to exist.

Going forward, several research directions appear to be promising. One pertains to extending the quantum TR dynamics investigations to more complex molecular dimers inside larger fullerenes, possibly starting with (HF)_{2}@C_{70}, and moving on to (H_{2}O)_{2}@C_{70} and (H_{2}O·HF)@C_{70} that have been synthesized already. Another aims to investigate the quantum many-body correlations between caged dipolar molecules such as H_{2}O and HF in H_{2}O@C_{60} and HF@C_{60} crystals, respectively, and explore the possibility of the emergence of dipole-ordered phases and ferroelectricity.^{32} These explorations need not be restricted to M@C_{60} (M = HF, H_{2}O) crystals but should include lower dimensional, 1D and 2D, arrays of these moieties, arranged in different geometries. Initial steps in this most intriguing direction have been taken,^{28,29} but virtually all the work still lies ahead. Finally, in trying to anticipate the main foci of future research, it is useful to keep in mind the following quote usually attributed to Niels Bohr: “Prediction is very difficult, especially about the future.”

## ACKNOWLEDGMENTS

Z.B. is grateful to the National Science Foundation for its partial support of this research through the Grant No. CHE-1566085. Joseph Cendagorta is gratefully acknowledged for preparing Fig. 1.

## REFERENCES

_{2}, HF, and H

_{2}O inside the fullerene C

_{60}

A. J. Horsewill has informed us that the observed 101 splitting in H_{2}O@C_{60} is such that the lower energy level is nondegenerate and the upper one is doubly degenerate. Due to a misprint, the opposite ordering was erroneously reported in Goh *et al*.^{86}