We present a methodology that, for the first time, allows rigorous quantum calculation of the inelastic neutron scattering (INS) spectra of a triatomic molecule in a nanoscale cavity, in this case, H_{2}O inside the fullerene C_{60}. Both moieties are taken to be rigid. Our treatment incorporates the quantum six-dimensional translation–rotation (TR) wave functions of the encapsulated H_{2}O, which serve as the spatial parts of the initial and final states of the INS transitions. As a result, the simulated INS spectra reflect the coupled TR dynamics of the nanoconfined guest molecule. They also exhibit the features arising from symmetry breaking observed for solid H_{2}O@C_{60} at low temperatures. Utilizing this methodology, we compute the INS spectra of H_{2}O@C_{60} for two incident neutron wavelengths and compare them with the corresponding experimental spectra. Good overall agreement is found, and the calculated spectra provide valuable additional insights.

## I. INTRODUCTION

The fullerene C_{60} has fascinated scientists with its exceptionally high symmetry, and the potential for encapsulating atoms and small molecules in its cavity, from the moment it was serendipitously discovered by Kroto *et al.*^{1} This tantalizing potential became a reality with the introduction of the “molecular surgery” procedure.^{2–4} The approach proved to be remarkably successful and versatile, and within a decade, it led to the synthesis of the light-molecule endofullerenes (LMEFs) H_{2}@C_{60},^{5,6} H_{2}O@C_{60},^{7} HF@C_{60},^{8} and, most recently, CH_{4}@C_{60}.^{9} These endofullerenes are now available with high purity and in macroscopic quantities, making possible a wide range of spectroscopic studies of their fundamental properties and potential practical applications.^{10}

A distinctive feature of H_{2}O@C_{60} and other LMEFs that has prompted a large number of experimental^{10} and theoretical studies^{11,12} are strong nuclear quantum effects (NQEs), which dominate the dynamics and spectroscopy of the guest molecules, particularly for low temperatures (typically ranging from 1.5 K to about 30 K), at which the spectroscopic measurements are typically performed. These species are the embodiments of many of the fundamental principles of quantum mechanics to the degree and with a clarity that are unrivaled in molecular systems amenable to experimental investigations. NQEs have multiple sources. One of them is the quantization of the translational center-of-mass (c.m.) degrees of freedom (DOFs) of the encapsulated molecules arising from their confinement inside the fullerene cavity (the textbook particle-in-a-box effect). Because of the tight confinement and low molecular masses of the guest molecules, the energy differences between their translational eigenstates are large relative to *kT* (*k* is the Boltzmann constant). The same holds for the quantized rotational states of these light molecules having one or more hydrogen atoms, due to their large rotational constants. Furthermore, the quantized translational and rotational DOFs of the guest molecule are coupled by the confining potential of the fullerene interior, resulting in a sparse and intricate translation–rotation (TR) energy level structure.^{11,12}

When a molecule, such as H_{2}O and H_{2}, has two symmetrically equivalent H atoms, it exhibits the phenomenon of nuclear spin isomerism that introduces additional prominent NQEs in its TR dynamics inside the fullerene cages. The two identical ^{1}H nuclei are fermions since their nuclear spin is 1/2. Satisfying the Pauli principle requires that the total molecular wave function of the caged molecule, its spatial and spin components, must be antisymmetric with respect to the exchange of the two fermions. A consequence of this is the particular entanglement of the spin and spatial quantum states, which results in nuclear spin isomers, *para* and *ortho*, of H_{2}O (and H_{2}), having total nuclear spins *I* = 0 and 1, respectively. 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.^{13} The constraints that the Pauli principle places on the rotational quantum numbers of the spin isomers of H_{2}O make its already sparse TR level structure inside the fullerenes even sparser, pushing the TR dynamics even deeper in the quantum regime.

Theoretical studies have thoroughly characterized the TR level structure of H_{2}O@C_{60}.^{11,12} The TR eigenstates of *para*- and *ortho*-H_{2}O in an (isolated) C_{60} cage with *I*_{h} symmetry have been obtained by means of fully coupled quantum calculations in 6D, treating both H_{2}O and C_{60} as rigid,^{14} and also in 9D, for flexible H_{2}O and its excited inter- and intramolecular vibrational states inside C_{60}^{15} (the 9D calculations in Ref. 16 were performed only for H_{2}O in the ground intramolecular vibrational state). These calculations^{14,15} demonstrated that the TR level structure of H_{2}O@C_{60} exhibits all of the key qualitative features previously identified for H_{2}@C_{60}.^{11,12,17,18} The purely translational eigenstates can be assigned with the quantum numbers of the 3D isotropic harmonic oscillator (HO), 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. For the H_{2}O rotational excitations, the asymmetric top quantum numbers $jkakc$ are employed. When the TR eigenstates are excited both translationally and rotationally, the orbital angular momentum **l** of the c.m. of H_{2}O and the rotational angular momentum **j** of H_{2}O couple to give the total angular momentum ** λ** =

**l**+

**j**, with

*λ*=

*l*+

*j*,

*l*+

*j*− 1, …, |

*l*−

*j*|. The values of

*l*are those allowed for the quantum number

*n*of the 3D isotropic HO. The TR states having the same quantum numbers

*n*and $jkakc$ are split into as many distinct levels as there are different values of

*λ*, each with the degeneracy 2

*λ*+ 1.

^{14}A remarkable and comprehensive infrared (IR) spectroscopic study probing these and other features of intra- and intermolecular vibrations of H

_{2}O in C

_{60}was published recently.

^{19}

Much of what is known experimentally about the quantum TR dynamics of molecules having one or more hydrogen atoms, such as H_{2}, HD, HF, and H_{2}O, encapsulated inside C_{60} has come from the inelastic neutron scattering (INS) spectroscopy of these endohedral complexes.^{8,20–25} The exceptional power and utility of the INS stem from two unique characteristics. One of them is the unusually large cross section for the incoherent neutron scattering from the hydrogen (^{1}H) nucleus; it is ∼15 times greater than for the nucleus of any other elements, including the isotope deuterium (D). As a result, INS is an exquisitely selective probe of the dynamics of the confined H-containing molecules. The peaks present in the measured INS spectra arise from the transitions between the TR eigenstates of the caged molecules induced by the incident neutrons. The second important advantage of the INS is that neutrons, unlike photons, can change the total nuclear spin of the guest molecule, allowing the observation of rotational Δ*j* = 1 transitions that are forbidden in optical, IR and Raman, spectroscopy, since they involve the *ortho*–*para* conversion.

The rich information content of the INS spectra of endohedral complexes has been underutilized for a long time. Its quantitative extraction is possible only with the help of high-level theory, given the intricacy of the TR dynamics. Complete analysis and assignment of the experimental INS spectra of the entrapped light molecules require theory that is capable of computing rigorously not only the TR excitation energies but also the intensities of the INS transitions. Together, they constitute a unique spectroscopic fingerprint of the system considered, enabling direct and more conclusive comparison with the experimental INS data. However, rigorous methodology for such calculations simply did not exist until rather recently.

In the past decade, we have developed an approach that enables rigorous and efficient quantum calculation of the INS spectra, i.e., the energies and the INS intensities of the TR transitions, of a diatomic molecule, initially homonuclear, such as H_{2}, confined in a nanocavity of an arbitrary shape, assuming that the guest molecule and the cavity are rigid.^{26,27} What sets this approach apart from all others in the literature is that it uses the quantum 5D TR wave functions of the entrapped diatomic molecule, which fully couple its translational and rotational DOFs, as the spatial parts of the initial and final states of the INS transitions. As a result, the simulated INS spectra exhibit an unprecedented degree of realism and reflect the complexity of the TR dynamics of the guest molecule in nanoconfinement. The first demonstration of this novel methodology involved the calculation of the INS spectra of *para*- and *ortho*-H_{2} in the small cage of the sII clathrate hydrate,^{26} in remarkably good agreement with the experimental data.^{28} Subsequent application involved the computation of the INS spectra of H_{2}@C_{60},^{29} for a range of temperatures. In the next step, this approach was extended to a heteronuclear diatomic molecule entrapped inside a nanoscale cavity.^{30} It was implemented to compute the INS spectra of HD in the small cage of the sII clathrate hydrate, which agreed well with the measured INS spectra of this system,^{28,31} and also in C_{60}.^{32} The most recent application was to the quantum simulation of the INS spectra of HF@C_{60},^{33} which was instrumental in the assignment of the measured spectrum.^{8} An alternative approach for calculating the INS spectrum of a homonuclear diatomic molecule inside a spherical nanocavity was reported by Mamone *et al.*^{34} based on the expansion of the confining potential into multipoles of the coupled rotational and translational angular variables. It was applied to an impressive simulation of the INS spectrum of H_{2}@C_{60}.

In addition to being indispensable for the analysis of the experimental INS spectra, the new methodology has led to some genuine conceptual advances. The most important among them is the surprising discovery of the INS selection rule, in the field where it has long been taken for granted, as an article of faith, that the incoherent INS spectroscopy of the vibrations of discrete molecular systems is not subject to any selection rules.^{35–38} However, rigorous quantum calculations revealed that certain transitions in the INS spectra of H_{2}@C_{60}^{29} and HD@C_{60}^{32} have zero-intensity (to 4–5 significant figures), strongly suggesting that they are forbidden according to some as yet unknown, and unsuspected, INS selection rule. This motivated the analytical derivation of the first ever selection rule in the INS vibrational spectroscopy, applicable to transitions originating from the ground TR state of *para*-H_{2}^{29} or HD^{32} confined in a near-spherical nanocavity, such as that of C_{60}. Within a year of its publication, this precedent-setting selection rule was validated by a joint experimental and theoretical INS investigation of H_{2}@C_{60}.^{39} Shortly afterward, the generalized selection rule was derived, which defines all possible forbidden INS transitions, originating from a variety of TR states, ground and excited, of *ortho*-H_{2}, *para*-H_{2}, and HD, inside a near-spherical nanocavity.^{40} Almost simultaneously, the same generalized selection rule was derived for H_{2}@C_{60}, but not HD@C_{60}, using group-theory arguments.^{41} A recent comprehensive review of the methodology for accurate quantum calculations of the INS spectra of diatomic molecules inside nanocavities and its applications is available.^{42}

Once the quantum methodology for computing the INS spectra of diatomic molecules in nanoconfinement was developed, it was natural to take the next step and extend the formalism to the INS spectra of triatomic molecules inside nanocavities, for which such a treatment is not available. Additional strong motivation for embarking on this endeavor came from the INS spectroscopic measurements of H_{2}O@C_{60} by Horsewill and co-workers.^{24,25} In particular, Ref. 25 presented the INS spectra of H_{2}O@C_{60} taken for several incident neutron wavelengths and a range of temperatures, as well as their analysis. These INS measurements 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, clear evidence that the symmetry of the environment felt by the caged H_{2}O must be lower than *I*_{h}. Similar symmetry breaking was observed for solid H_{2}@C_{60}^{43} and HF@C_{60}.^{8} Quantitative first-principles explanation of this phenomenon in terms of an electrostatic interaction model was developed by us in Refs. 11, 44, and 45.

In this paper, the methodology is introduced, which makes possible rigorous quantum computation of the INS spectra of a triatomic molecule in a nanocavity, specifically H_{2}O inside the C_{60} cage. The present approach builds on our previous treatment of the INS spectra of nanoconfined diatomic molecules.^{26,27,30,33} However, going from a diatomic molecule to a (nonlinear) triatomic molecule inside a cavity increases the complexity of the treatment significantly, which also accounts for the effects of symmetry breaking. With this methodology, the INS spectra of H_{2}O@C_{60} are calculated for two incident neutron wavelengths and compared to their experimental counterparts in Ref. 25. There is good overall agreement between the theory and experiment, and the former sheds additional light on the measured INS spectra.

## II. THEORY

### A. Basic equations

As in our previous treatments of the INS spectra of nanoconfined diatomic molecules (H_{2}, HD, and HF),^{26,27,30,33} the starting point of the methodology for calculating the INS spectrum of H_{2}O@C_{60} in this work is the standard expression for the neutron scattering double differential cross section in the first Born approximation,^{46,47}

where

and

In Eqs. (1)–(3), |*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}, $\kappa \u20d7=k\u20d7\u2212k\u20d7\u2032$, with $k\u20d7$ and $k\u20d7\u2032$ being the wave vectors of the incident and the scattered neutrons, respectively, *ℏω* = *E* − *E*′ = *ℏ*^{2}(*k*^{2} − *k*′^{2})/(2*m*), with *m* being the neutron mass, $B\u0302n$ is the scattering length operator, and $r\u20d7n$ is the position of nucleus *n*. For a system of protons,

where $\sigma \u20d72$ represents the neutron spin and $S\u20d7n$ is the spin of proton *n*. The coefficients *b*_{1} and *b*_{2} can be found in Ref. 47.

### B. Treatment of H_{2}O@C_{60}

The coordinates used in the treatment of the INS spectrum of H_{2}O@C_{60} described below are depicted in Fig. 1. They are as follows:

In Eq. (5), $R\u20d7$ fixes the position of the H_{2}O c.m. relative to a space-fixed Cartesian system with the origin at the center of the C_{60} cage, while $\rho \u20d7RO$, $\rho \u20d7RH1$, and $\rho \u20d7RH2$ connect the c.m. of H_{2}O to O, H_{1}, and H_{2} atoms, respectively. In addition, the vector pointing from H_{1} to H_{2} is denoted as $\rho \u20d7H2$, and the vector from the c.m. of H_{2}O to the midpoint of $\rho \u20d7H2$ is denoted as $\rho \u20d7RHH$. It is easy to see that

It should be noted that $\rho \u20d7RO$, $\rho \u20d7RH1$, $\rho \u20d7RH2$, $\rho \u20d7H2$, and $\rho \u20d7RHH$ all lie in the plane defined by the three atoms of H_{2}O.

With the above, $Mfi$ in Eq. (3) can be written as

where

The above equations already reflect that for ^{16}O atoms, $b2O=0$.^{47} Moreover, the ^{17}O and ^{18}O isotopes have negligible abundance.

For H_{2}O confined inside the C_{60} cavity, the total wave functions of the initial (|*i*⟩) and final states (|*f*⟩) in Eqs. (1)–(3) can be written in the following product form:

In Eq. (12), *ϕ*, *θ*, and *χ* are the Euler angles that specify the orientation of the body-fixed (BF) axes of H_{2}O relative to the axes fixed to the C_{60} frame.^{14} $|I\u0302\tau \u3009$ and $|\Psi \tau 6D(R\u20d7,\varphi ,\theta ,\chi )\u3009(\tau =i,f)$ are the nuclear-spin and 6D spatial wave functions, respectively, of the caged H_{2}O molecule, respectively. $|\Psi \tau 6D(R\u20d7,\varphi ,\theta ,\chi )\u3009$ are the eigenstates of the 6D TR Hamiltonian for rigid H_{2}O molecules in (rigid) C_{60} obtained as described below. Hereafter, $|\Psi \tau 6D(R\u20d7,\varphi ,\theta ,\chi )\u3009$ are denoted as $|\Psi \tau 6D\u3009(\tau =i,f)$.

Taking advantage of the product form of the initial and final states in Eq. (12), Eq. (7) can be written as

The above equation for $Mfi$ contains two distinct types of matrix elements: (*i*) those involving the nuclear-spin wave functions |*I*_{τ}⟩ (*τ* = *i*, *f*) and (*ii*) those involving the 6D spatial (TR) wave functions $|\Psi \tau 6D\u3009(\tau =i,f)$. In the following, they are referred to as the spin and spatial matrix elements, respectively, and are evaluated separately below.

### C. Evaluation of the spin matrix elements

For H_{2}O, the nuclear-spin state $|I\u0302\u3009H2O$, either *para*- (*I* = 0) or *ortho*- (*I* = 1), can be written as

The nuclear-spin state |0⟩_{O} of the oxygen atom is a singlet, while the nuclear-spin wave functions $|I\u0302\u3009H2$ are the same as those of the H_{2} molecule, singlet *para*- (*I* = 0) or triplet *ortho*- (*I* = 1) states. Consequently, for simplicity, $|I\u0302\u3009H2O$ is hereafter referred to as $|I\u3009H2$.

In Eq. (13), it is necessary to compute two spin matrix elements that describe the coupling between the neutron spin and the nuclear spins of the H_{1} and H_{2} protons,

In the above equation, $S\u20d7n$ is replaced by $i\u20d7n(n=1,2)$, and the total spin $I\u20d7\u2261S\u20d7H1+S\u20d7H2$ is defined. In addition, the nuclear-spin wave function |*I*_{τ}⟩ is written explicitly as $|I\tau MI\tau \u3009(\tau =i,f)$.

Following closely the treatment of the spin matrix elements for H_{2} in a nanocavity,^{27} one obtains the equations for the nuclear-spin transition *I*_{i} → *I*_{f} of the caged H_{2}O molecule,

where

Equation (16) represents the final expressions for the neutron scattering cross sections for the transitions between different nuclear-spin states of H_{2}O. What remains is the evaluation of the spatial matrix elements in Eq. (17).

In Eq. (16), the coherent and incoherent scattering cross sections for the proton, $\sigma cohH$ and $\sigma incH$, respectively, and the coherent scattering cross section for the O atom, $\sigma cohO$, are^{48}

The spatial matrix elements in Eq. (17) are evaluated in Sec. II E. However, in order to set the stage for this, in Sec. II D, we first describe the computational approach used to calculate the 6D TR eigenstates $|\Psi \tau 6D\u3009(\tau =i,f)$ of H_{2}O inside the C_{60} cage. They are the spatial components of the states |*i*⟩ and |*f*⟩ in Eq. (12), which are the initial and final states, respectively, of the INS transitions.

### D. Calculation of the 6D translation–rotation eigenstates of H_{2}O@C_{60} in the presence of symmetry breaking

The 6D TR eigenstates $|\Psi \tau 6D\u3009(\tau =i,f)$ of H_{2}O@C_{60} employed in our calculations of the INS spectra of this system incorporate the effects of symmetry breaking. The electrostatic interaction model, which we have developed to accomplish this and account quantitatively for the magnitude of the 1:2 splitting of the 1_{01} ground state of *ortho*-H_{2}O observed in the INS spectra of solid H_{2}O@C_{60} at low temperatures,^{24,25} is described in detail in Refs. 11, 44, and 45, and only its salient features are summarized here.

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.^{49} In the dominant one, referred to as the P orientation, the electron-rich double bonds shared by two hexagons (denoted 6:6) of one C_{60} unit face directly the electron-poor pentagons of the neighboring cages.^{49,50} The other, H orientation has the electron-rich 6:6 bonds of one C_{60} immediately adjacent to the electron-poor hexagonal faces of the neighboring units.^{49} Below 90 K, and at ambient pressure, the relative proportion of molecules in the P and H orientations is 5:1.^{49} 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},^{43,51} capable of splitting the three-fold degeneracy of the 1_{01} ground state of *ortho*-H_{2}O.

Our approach^{11,44,45} focuses on the smallest fragment of solid C_{60} that can give rise to *S*_{6} symmetry in its center, comprised of 12 nearest-neighbor (NN) cages around the central cage, in either P or H crystal orientation. The cages are treated as rigid and assumed to have *I*_{h} symmetry. Only the central cage is occupied by the H_{2}O molecule, while the 12 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*) (M = H_{2}O, H_{2}, HF) 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 density functional theory (DFT) calculation, and the charge density on M is related to its body-fixed quadrupole moments reported in the literature. Thus, the M–NN electrostatic interaction, denoted as *V*_{ES}, contains no adjustable parameters.

For this electrostatic interaction model, the 6D TR Hamiltonian of (rigid) H_{2}O inside the central cage of the (rigid) 13-cage fragment of H_{2}O@C_{60}(*s*) can be written as

where $T\u0302$ is the operator associated with the TR kinetic energy of H_{2}O, $VH2O\u2212C60$ is the part of the 6D intermolecular PES arising from the H_{2}O–central cage interaction, and *V*_{ES} is the H_{2}O–NN electrostatic interaction above. $VH2O\u2212C60$ is defined in our previous study of H_{2}O@C_{60};^{14} it is constructed as a sum over the pairwise interactions, modeled with the Lennard-Jones 12-6 potentials, of each atom of H_{2}O with each atom of C_{60}. In principle, one could try to represent $VH2O\u2212C60$ by a polynomial/multipole expansion, as was done for He@C_{60}.^{52} However, it remains to be seen how well would this approach work for a 6D PES. Only the leading term in the multipole expansion of the H_{2}O–NN electrostatic interaction is retained in *V*_{ES}.^{44,45} Owing to the *S*_{6} point-group symmetry of the 13-cage (central + NN) crystal fragment,^{53,54} that term is the quadrupolar one,

where $Qm(2)$ constitute the electric-quadrupole spherical tensor of H_{2}O@C_{60} and $Im(2)$ constitute the electric-field-gradient tensor arising from the charges on the NN cages. It has been shown that retaining just the leading quadrupole terms in *V*_{quad} suffices to quantitatively account for the manifestations of symmetry breaking observed in all the three M@C_{60}(*s*) (M = H_{2}O, H_{2}, and HF) species considered thus far.^{44,45} We demonstrated that the P crystal orientation, besides being dominant at low temperatures, causes much larger TR level splittings than the H crystal orientation^{44,45} and is therefore the one employed in the present calculations. The geometry of the 13-cage fragment of H_{2}O@C_{60}(*s*) is fully defined in the supplementary information of Ref. 44.

The only, but computationally significant, modification in the above procedure made in this work is that $R\u20d7$, the position vector of the c.m. of H_{2}O, is expressed in terms of Cartesian coordinates {*x*, *y*, *z*}, instead of the spherical polar coordinates {*R*, *β*, *α*} (where $R\u2261|R\u20d7|$) used in our past work on H_{2}O@C_{60}.^{14,15} This allows us to use the 3D direct-product discrete variable representation (DVR)^{55} basis {|*X*_{α}⟩|*Y*_{β}⟩|*Z*_{γ}⟩} in the *x*, *y*, *z* coordinates (the sinc-DVR^{56}) as the basis for the translational c.m. degrees of freedom (DOFs) of H_{2}O; it is labeled by the grid points {*X*_{α}}, {*Y*_{β}}, and {*Z*_{γ}}, at which the respective DVR basis functions are localized. We showed previously^{27} that the use of this 3D DVR is crucial for the efficient computation of the spatial matrix elements, as it greatly simplifies the otherwise rather forbidding equations and drastically reduces the computational effort involved in their evaluation. The symmetric-top eigenfunctions, i.e., the normalized Wigner rotation functions ${|j,m,k\u3009\u2261(2j+1)/8\pi 212Dmkj(\varphi ,\theta ,\chi )*}$, serve as the basis for the angular DOFs of H_{2}O. The Euler angles *ϕ*, *θ*, and *χ* specify the orientation of the body-fixed (BF) axes of H_{2}O relative to the axes fixed to the C_{60} frame. The BF *z* axis of H_{2}O is taken to be the symmetry axis of the molecule, the *x* axis is taken to be the other principal axis of the molecule lying in the molecular plane, and the *y* axis is perpendicular to the molecular plane and is such that the axis system is right-handed.

Together, the above translational and rotational bases constitute the 6D direct-product basis set {|*X*_{α}⟩|*Y*_{β}⟩|*Z*_{γ}⟩|*jmk*⟩}. The eigenstates $|\Psi i6D\u3009$ obtained by diagonalizing in this basis the 6D TR Hamiltonian for H_{2}O@C_{60} in Eq. (19) are

### E. Evaluation of the spatial matrix elements

Let the direction of $\kappa \u20d7$ in Eq. (13) be defined by the angles $(\theta \kappa \u20d7,\varphi \kappa \u20d7)$. On the other hand, the directions of the three vectors $\rho \u20d7\eta $ (*η* = RO, RH_{1}, and RH_{2}) in Fig. 1 are defined by the angles (*θ*_{η}, *ϕ*_{η}), respectively, in the cage-fixed Cartesian axis system. Therefore, we make use of the well known expression as follows:^{57}

where $\rho \eta \u2261|\rho \eta \u20d7|$ and *j*_{l}(*κρ*_{η}) are the spherical Bessel functions.^{58}

We now consider explicitly the term in Eq. (13) associated with the O atom. Then, following Eq. (22), we can write

Taking the initial angles of $\rho \u20d7RO$ to be $(\theta \rho RO0,\varphi \rho RO0)=(\pi ,\pi )$, under the rotation Ω ≡ (*ϕ*, *θ*, *χ*), we have

where $Dmm\u2033l(R)*$ is the complex conjugate of an element of the Wigner D-matrix. With this, Eq. (23) can be written as

The change to upper-case indices *L*, *M*, and *M*″ in the above equation is made in order to avoid possible confusion with the indices used in the further derivation.

In the same vein, expressions analogous to that in Eq. (25) for $\rho \u20d7RO$ can be obtained for $\rho \u20d7RH1$ and $\rho \u20d7RH2$. From the information in Table I, one can deduce that their initial angles are $(\theta \rho \u20d7RH10,\varphi \rho \u20d7RH10)=(\theta \rho \u20d7RH10,0)$ and $(\theta \rho \u20d7RH20,\varphi \rho \u20d7RH20)=(\theta \rho \u20d7RH10,\varphi \rho \u20d7RH10+\pi )=(\theta \rho \u20d7RH10,\pi )$, respectively. With this, one can obtain the following expressions:

and

Based on Eqs. (26) and (27), and noticing that $\rho RH1=\rho RH2$, the following factors in Eq. (13) can be expressed as

Label . | X . | Y . | Z . |
---|---|---|---|

O | 0.0000 | 0.0000 | −0.065 638 07 |

H | 0.7575 | 0.0000 | 0.520 861 93 |

H | −0.7575 | 0.0000 | 0.520 861 93 |

Label . | X . | Y . | Z . |
---|---|---|---|

O | 0.0000 | 0.0000 | −0.065 638 07 |

H | 0.7575 | 0.0000 | 0.520 861 93 |

H | −0.7575 | 0.0000 | 0.520 861 93 |

In the following, we make use of the equation

which follows from the well-known expression for the integral over the product of three rotation matrices.^{57,59} This allows us to obtain the matrix elements involving the expressions in Eqs. (25) and (28) for the initial and final angular states |*jmk*⟩ and |*j*′*m*′*k*′⟩, respectively,

In order to make the notation more compact, we now define the following quantities:

where the matrix elements on the right-hand side are given in Eq. (30).

From the properties of the spherical harmonics,

Taking into account Eqs. (21) and (34), and the DVR points $X\alpha ,Y\beta ,Z\gamma $, the following final expressions can be derived for the entities on the left-hand side of Eq. (17):

Equation (35) represents the key result of this derivation, as it allows the evaluation of the matrix elements in Eq. (17). These in turn enter Eq. (16), which gives the final equations for the neutron scattering cross sections for the transitions between different nuclear-spin states of the encapsulated H_{2}O.

### F. Calculation of the INS spectra of H_{2}O@C_{60}

The experimental INS spectra of H_{2}O@C_{60}, with which in Sec. III we compare our computed INS spectra, were taken from a polycrystalline sample,^{25} in which the fullerene cages are randomly oriented with respect to the incoming neutron beam. Therefore, to achieve a more realistic comparison with the measured spectra, the INS spectra computed by the approach described above are averaged over all possible orientations of C_{60} by the procedure described previously^{27} and implemented in our quantum simulations of the INS spectra of H_{2}/HD@C_{60}^{29,32} and HF@C_{60}.^{33} The numerical integration involved in the powder averaging of $S(\kappa \u20d7,\omega )$ in Eq. (2) over the directions of $\kappa \u20d7$, defined by the angles $(\theta \kappa \u20d7,\varphi \kappa \u20d7)$, to yield *S*(*κ*, *ω*), is carried out using the ten-point Gauss–Legendre quadrature in *θ*_{κ} and the 20-point Gauss–Chebyshev quadrature in *ϕ*_{κ}. Moreover, these spectra were measured on the IN4C spectrometer, which allows the spectra to be recorded for a broad range of $\kappa \u2261|\kappa \u20d7|$ values. This means that *S*(*κ*, *ω*) needs to be integrated also over *κ* in order to obtain the desired *S*(*ω*), as described in Ref. 32.

The above integration of the computed INS spectrum over both the angles $(\theta \kappa \u20d7,\varphi \kappa \u20d7)$ and *κ* is performed for each INS transition in the range of excitation energies considered. The procedure is computationally time consuming but necessary and results in a realistic simulation of the experimental INS spectra.

In the calculations of the 6D TR eigenstates of H_{2}O@C_{60} described in Sec. II D, the mass of the H_{2}O molecule equal to 18.010 56 amu is used, together with the moments of inertia of H_{2}O about the BF axes, *I*_{x} = 0.604 72, *I*_{y} = 1.8156, and *I*_{z} = 1.161 64 amu Å^{2} (taken to be those associated with the measured rotational constants of H_{2}O in the ground vibrational state:^{60} *A* = 27.877, *C* = 9.285, and *B* = 14.512 cm^{−1}, respectively). The dimension of the sinc-DVR basis is 14 for each of the three Cartesian coordinates *x*, *y*, and *z*, spanning the range from −1 to +1 Å in each coordinate. In addition, all angular |*j*, *m*, *k*⟩ basis functions for 0 ≤ *j* ≤ 6 are included in the 6D basis.

## III. RESULTS AND DISCUSSION

The energies of the transitions in the INS spectra discussed here do not exceed 10 meV, while the first translationally excited state of *para*-H_{2}O on the PES used in this work is at 21.4 meV. Therefore, all of the TR states involved are purely rotational, in the ground translational state. For them, the translational quantum numbers *n* and *l* of the 3D isotropic HO^{14,15} are equal to zero and are not included in the assignments. Hence, hereafter, the TR states are assigned as $(jkakc,\lambda ,|m\lambda |)$, where |*m*_{λ}| is the absolute value of the component of *λ* along the *C*_{3} axis of the *S*_{6} point group.^{45} Both *λ* and |*m*_{λ}| remain very good quantum numbers in the presence of the quadrupolar H_{2}O–NN interactions, and the latter labels the sublevels arising from the symmetry breaking.^{45}

In the analysis of the INS spectra, below, we will often refer to Table 3 of Ref. 45, which gives the energies and other properties of the calculated TR eigenstates of H_{2}O@C_{60} in the presence of symmetry breaking, computed as described in Sec. II D.

The INS spectra of H_{2}O@C_{60}(*s*) computed for the incident neutron wavelengths *λ*_{n} = 3.0 Å and *λ*_{n} = 2.3 Å, at 0 K, are shown in Figs. 2 and 3, respectively. They correspond to the experimental INS spectra of Goh and co-workers^{25} recorded for the same neutron wavelengths and at 1.6 K, which are displayed in Figs. 2 and 5, respectively, of Ref. 25. The benefits of using different incident neutron wavelengths are apparent from the comparison of the INS spectra, calculated or measured, for these two incident wavelengths. Decreasing the wavelength, i.e., increasing the energy of the incoming neutrons, enables excitation of higher-energy TR states of the guest H_{2}O. However, this also lowers the energy resolution of the INS spectrum, which, as discussed below, leaves some important features of the TR energy level structure unresolved—which are resolved for longer incident wavelengths.

The stick spectra appearing in Figs. 2 and 3 are calculated for the 3:1 ratio of *ortho*-H_{2}O to *para*-H_{2}O. In order to facilitate visual comparison with the measured INS spectra, the stick spectra in both figures are convolved with the instrumental resolution function, resulting in smooth spectra depicted by the full red line. Also shown in Figs. 2 and 3 are smooth INS spectra depicted by the dashed blue line; they correspond to the stick spectra simulated for the 1:1 ratio of *ortho*-H_{2}O to *para*-H_{2}O (not shown for clarity) convolved with the instrumental resolution function. The purpose of computing the INS spectra for 3:1 and 1:1 ratios of *ortho*-H_{2}O to *para*-H_{2}O is to roughly mimic the corresponding experimental INS spectra recorded in the first hour and eighth hour after cooling, respectively.^{25} Spontaneous interconversion from *ortho*-H_{2}O to *para*-H_{2}O over several hours observed in H_{2}O@C_{60} decreases the initial 3:1 *ortho*- to *para*-H_{2}O ratio.^{25}

The quantity appearing on the abscissa of Figs. 2 and 3 is the neutron energy transfer Δ*E* = *E* − *E*′, where *E* and *E*′ are the energies of the incident and scattered neutron, respectively. Δ*E* is positive in neutron-energy (NE) loss, when the incoming neutron transfers a part of its energy to the sample and excites the caged H_{2}O from a lower to a higher TR energy level. In NE gain, when Δ*E* is negative, the scattered neutron gains energy when H_{2}O makes a downward transition from a higher to a lower TR energy level.

The peaks in Figs. 2 and 3 have labels 1–6 and 1–7, respectively. The corresponding transitions 1–6 are given in Table II.

Label . | ΔE (meV)
. | Intensity (b) . | Initial $(jkakc,\lambda ,|m\lambda |)$ . | Final $(jkakc,\lambda ,|m\lambda |)$ . | Final g
. |
---|---|---|---|---|---|

1 | 0.52 | 12.11 | (1_{01}, 1, 0) | (1_{01}, 1, 1) | 2 |

2 | 1.91 | 7.67 | (1_{01}, 1, 0) | (1_{11}, 1, 1) | 2 |

1.96 | 4.67 | (1_{01}, 1, 0) | (1_{11}, 1, 0) | 1 | |

3 | 2.42 | 38.94 | (1_{01}, 1, 0) | (1_{10}, 1, 1) | 2 |

2.94 | 0.000 | (1_{01}, 1, 0) | (1_{10}, 1, 0) | 1 | |

4 | 2.61 | 97.10 | (0_{00}, 0, 0) | (1_{01}, 1, 0) | 1 |

3.13 | 171.07 | (0_{00}, 0, 0) | (1_{01}, 1, 1) | 2 | |

5 | 5.81 | 19.37 | (1_{01}, 1, 0) | (2_{02}, 2, 0) | 1 |

5.93 | 29.52 | (1_{01}, 1, 0) | (2_{02}, 2, 1) | 2 | |

6.35 | 0.04 | (1_{01}, 1, 0) | (2_{02}, 2, 2) | 2 | |

6 | 6.92 | 11.37 | (1_{01}, 1, 0) | (2_{12}, 2, 0) | 1 |

7.05 | 15.51 | (1_{01}, 1, 0) | (2_{12}, 2, 1) | 2 | |

7.40 | 0.76 | (1_{01}, 1, 0) | (2_{12}, 2, 2) | 2 |

Label . | ΔE (meV)
. | Intensity (b) . | Initial $(jkakc,\lambda ,|m\lambda |)$ . | Final $(jkakc,\lambda ,|m\lambda |)$ . | Final g
. |
---|---|---|---|---|---|

1 | 0.52 | 12.11 | (1_{01}, 1, 0) | (1_{01}, 1, 1) | 2 |

2 | 1.91 | 7.67 | (1_{01}, 1, 0) | (1_{11}, 1, 1) | 2 |

1.96 | 4.67 | (1_{01}, 1, 0) | (1_{11}, 1, 0) | 1 | |

3 | 2.42 | 38.94 | (1_{01}, 1, 0) | (1_{10}, 1, 1) | 2 |

2.94 | 0.000 | (1_{01}, 1, 0) | (1_{10}, 1, 0) | 1 | |

4 | 2.61 | 97.10 | (0_{00}, 0, 0) | (1_{01}, 1, 0) | 1 |

3.13 | 171.07 | (0_{00}, 0, 0) | (1_{01}, 1, 1) | 2 | |

5 | 5.81 | 19.37 | (1_{01}, 1, 0) | (2_{02}, 2, 0) | 1 |

5.93 | 29.52 | (1_{01}, 1, 0) | (2_{02}, 2, 1) | 2 | |

6.35 | 0.04 | (1_{01}, 1, 0) | (2_{02}, 2, 2) | 2 | |

6 | 6.92 | 11.37 | (1_{01}, 1, 0) | (2_{12}, 2, 0) | 1 |

7.05 | 15.51 | (1_{01}, 1, 0) | (2_{12}, 2, 1) | 2 | |

7.40 | 0.76 | (1_{01}, 1, 0) | (2_{12}, 2, 2) | 2 |

The only states populated at 0 K in the calculations (and 1.6 K in the experiments^{25}) are the (0_{00}, 0, 0) ground state of *para*-H_{2}O and the (1_{01}, 1, 0) ground state of *ortho*-H_{2}O. The latter is the lower-energy component of the three-fold degenerate 1_{01} ground state of *ortho*-H_{2}O split in the 1:2 pattern by the symmetry breaking.^{25,44} Therefore, all INS transition present in the calculated spectra in Figs. 2 and 3, as well as in the measured INS spectra shown in Figs. 2 and 5 of Ref. 25, must originate from one of these two states.

We now discuss the INS spectrum in Fig. 2 calculated for the incident neutron wavelength *λ*_{n} = 3.0 Å. Only one peak is present in the NE gain part of the spectrum, at −2.61 meV; its origin is in the transition from the (1_{01}, 1, 0) ground state of *ortho*-H_{2}O to the (0_{00}, 0, 0) ground state of *para*-H_{2}O. This provides clear evidence for the existence of two nuclear spin isomers of water, *ortho*-H_{2}O, and *para*-H_{2}O. The same peak is present in the measured INS spectrum^{25} at −2.56 meV. Its intensity decays with time, evident in Fig. 2 of Ref. 25, as a result of slow conversion from *ortho*- to *para*-H_{2}O. This agrees with the fact that the height of the peak at −2.61 meV in Fig. 2 computed for the 3:1 ratio of *ortho*-H_{2}O to *para*-H_{2}O is greater than that obtained for the 1:1 ratio of these two spin isomers.

In the NE loss, the lowest-energy calculated peak at 0.52 meV, labeled 1 in Fig. 2, arises from the transition (**1**_{01}, 1, 0) → (**1**_{01}, 1, 1) between the lower- and higher-energy components of the 1:2 split ground state of *ortho*-H_{2}O. This peak is not identified in the measured INS spectrum since it appears as barely a shoulder of the massive peak centered at 0 meV.^{25}

Label 2 applies to the peak, which according to Table II corresponds to the pair of *ortho*–*para* transitions, (**1**_{01}, 1, 0) → (**1**_{11}, 1, 1) and (**1**_{01}, 1, 0) → (**1**_{11}, 1, 0), at 1.91 and 1.96 meV, respectively.

The peak at 2.42 meV labeled 3 is associated with the *ortho*–*ortho* transition (**1**_{01}, 1, 0) → (**1**_{10}, 1, 1) (Table II). However, it is evident from Table 3 of Ref. 45 that the final state 1_{10} is actually split by the symmetry breaking into *two* components, (1_{10}, 1, 1) at 5.04 meV and (1_{10}, 1, 0) at 5.55 meV. Therefore, another peak should be present in the spectrum at 2.94 meV due to the transition (**1**_{01}, 1, 0) → (**1**_{10}, 1, 0), but it is *absent* since its calculated intensity is *zero* (0.000000). This is reminiscent of the situation for H_{2}@C_{60}, where our calculations revealed that certain INS transitions have zero intensity as well. This led us to postulate and derive the first ever, and entirely unexpected, selection rule in the INS spectroscopy of discrete molecular compounds.^{29} Soon thereafter, the predicted precedent-setting selection rule was validated by a joint experimental and theoretical INS investigation of H_{2}@C_{60}.^{39} It remains to be seen whether a similar INS selection rule exists for H_{2}O@C_{60}.

In the corresponding measured INS spectrum in Fig. 2 of Ref. 25, this *ortho*–*ortho* transition, marked as 1_{10}, is assigned to a shoulder just below 2 meV. However, as discussed above, according to our calculations, this is where the two *ortho*–*para* transitions (**1**_{01}, 1, 0) → (**1**_{11}, 1, 1) and (**1**_{01}, 1, 0) → (**1**_{11}, 1, 0) appear (label 2).

Label 4 refers to the two peaks at 2.61 and 3.13 meV, respectively, associated with the transitions out of the (**0**_{00}, 0, 0) ground state of *para*-H_{2}O to the two components of the ground state of *ortho*-H_{2}O, (**1**_{01}, 1, 0) and (**1**_{01}, 1, 1), resulting from the symmetry breaking. In the experimental INS spectrum, these two transitions are observed at slightly different energies, 2.50 and 3.02 meV, respectively.^{25} However, the measured and calculated splittings of the ground state of *ortho*-H_{2}O are the same, 0.52 meV.^{44,45}

Label 5 applies to three *ortho*–*para* transitions at 5.81, 5.93, and 6.35 meV, respectively: (**1**_{01}, 1, 0) → (**2**_{02}, 2, 0), (**1**_{01}, 1, 0) → (**2**_{02}, 2, 1), and (**1**_{01}, 1, 0) → (**2**_{02}, 2, 2), respectively (Table II). They arise from the symmetry breaking of the final rotational state 2_{02} of *para*-H_{2}O into three components (Table 3 of Ref. 45). Only two of these INS transitions, at 5.81 and 5.93 meV, are visible in the calculated spectrum in Fig. 2 since the intensity of the third at 6.35 meV is extremely small (Table II). In the INS spectrum measured for *λ*_{n} = 2.3 Å, this *ortho*–*para* transition is assigned to the peak at 5.68 meV (Table 1 in Ref. 25).

Finally, label 6 denotes the three *ortho*–*ortho* transitions at 6.92, 7.05, and 7.40 meV, respectively: (**1**_{01}, 1, 0) → (**2**_{12}, 2, 0), (**1**_{01}, 1, 0) → (**2**_{12}, 2, 1), and (**1**_{01}, 1, 0) → (**2**_{12}, 2, 2), respectively (Table II). They are due to the symmetry breaking of the final rotational state 2_{12} of *ortho*-H_{2}O into three components (Table 3 of Ref. 45). Of these three INS transitions, only two, at 6.92 and 7.05 meV, are visible in the calculated spectrum in Fig. 2 since the intensity of the third at 7.40 meV is very small (Table II). In the INS spectrum measured for *λ*_{n} = 2.3 Å, this *ortho*–*ortho* transition is assigned to the peak at 6.51 meV (Table 1 in Ref. 25).

In the following, we discuss briefly the INS spectrum in Fig. 3 calculated for the incident neutron wavelength *λ*_{n} = 2.3 Å. In its main features, including the transitions with labels 1–6, it bears close resemblance to the spectrum in Fig. 2 computed for *λ*_{n} = 3.0 Å. An additional band appears, labeled 7, which is centered at 9.6 meV. It is associated with the peaks at 9.53, 9.66, and 10.01 meV, respectively, corresponding to the three transitions out of the (**0**_{00}, 0, 0) ground state of *para*-H_{2}O to the final states (**2**_{12}, 2, |*m*_{λ}|), |*m*_{λ}| = 0, 1, and 2 of *ortho*-H_{2}O. In the INS spectrum measured for *λ*_{n} = 2.3 Å, this *para*–*ortho* transition is assigned to the peak at 9.20 meV (Table 1 in Ref. 25).

As mentioned earlier, the ability to access higher-energy TR states comes at price of a lower resolution of the INS spectrum. In particular, the feature labeled 4 arising from the transitions out of the (**0**_{00}, 0, 0) ground state of *para*-H_{2}O to the (**1**_{01}, 1, 0) and (**1**_{01}, 1, 1) components of the ground state of *ortho*-H_{2}O, which appears as two well-resolved peaks in the *λ*_{n} = 3.0 Å spectrum in Fig. 2, is unresolved in the *λ*_{n} = 2.3 Å spectrum. The same is true for the corresponding experimental INS spectra shown in Figs. 2 and 5 of Ref. 25.

Overall, there is good correspondence between the INS spectra calculated in this work and the measured INS spectra reported in Ref. 25. The energies of the calculated transitions are typically higher by a fraction of a wavenumber than the corresponding ones in the experimental INS spectra. This undoubtedly reflects the shortcomings of the PES employed in the calculations, not surprising given its rather crude pairwise-additive character. On the other hand, the calculated INS spectra provide valuable insights that are obscured in the measured spectra. Thus, they reveal that most bands in the measured spectra in fact comprise 2–3 unresolved individual transitions, largely caused by symmetry breaking, whose intensity can vary greatly. They also reveal transitions, such as (**1**_{01}, 1, 0) → (**1**_{01}, 1, 1) at 0.52 meV, that are hidden as shoulders of much more intense transitions and are therefore difficult to identify with confidence in the measured INS spectra. Moreover, one calculated transition, (**1**_{01}, 1, 0) → (**1**_{10}, 1, 0), is shown to have zero intensity, hinting at the possible existence of an INS selection rule in this system, akin to the one we uncovered for H_{2}@C_{60}.^{29,39}

## IV. CONCLUSIONS

We have presented a methodology that allows rigorous quantum calculation of the INS spectra of H_{2}O inside C_{60}, with both moieties taken to be rigid. With relatively straightforward modifications, the method can be applied to simulating the INS spectra of other triatomic molecules confined inside arbitrarily shaped nanocavities. This methodology represents a natural extension of the quantum treatment pioneered by us previously to compute the INS spectra of diatomic molecules in nanocavities,^{26,27,30,33} which was applied successfully to diatomics nanoconfined inside the cages of C_{60} and clathrate hydrates.^{11,42} The treatment in this paper incorporates the quantum 6D TR wave functions of the entrapped water molecule, which couple its translational and rotational DOFs, as the spatial parts of the initial and final states of the INS transitions. This allows the simulated INS spectra to reflect in full the complexity of the TR dynamics of the guest molecule in nanoconfinement. Moreover, the 6D TR eigenstates are calculated in a way that accounts quantitatively for the effects of the symmetry breaking.^{44,45} As a result, the computed INS spectra exhibit characteristic features present in the INS spectra measured for solid H_{2}O@C_{60} at low temperatures,^{25} which arise from symmetry breaking.

In this paper, the newly developed methodology is utilized to compute the INS spectra of H_{2}O@C_{60} for two incident neutron wavelengths, *λ*_{n} = 3.0 Å and *λ*_{n} = 2.3 Å, at 0 K. Their comparison with the experimental INS spectra^{25} recorded for the same incident neutron wavelengths at 1.6 K reveals very good overall match between the theory and experiment. However, the calculated INS spectra lead to deeper understanding and interpretation of their experimental counterparts. What emerges from the former is that underneath most of the bands present in the measured spectra, there are 2–3 individual transitions, mostly stemming from symmetry breaking, which can differ greatly in intensity, that the experiments do not resolve. Moreover, the simulated spectra are able to reveal transitions that in the measured spectra are hidden in the shoulders of more intense transitions, and whose existence, and assignment, therefore cannot be established unambiguously. An intriguing result of our calculations is that one calculated *ortho*–*ortho* transition, (**1**_{01}, 1, 0) → (**1**_{10}, 1, 0), has zero intensity. This suggests the possibility that an INS selection rule exists in this system that remains to be identified, analogous to the one we uncovered for H_{2}@C_{60}.^{29,39} Future research will explore this fundamental issue in detail.

The accuracy of the spectra computed by the methodology developed in this work is limited largely by the quality of the intermolecular PES employed. In the case of H_{2}O@C_{60}, clearly, there is room for improvement, and we hope that the present results will add to the motivation for computing a high-quality endohedral PES of this system. Moreover, the recent synthesis of CH_{4}@C_{60}^{9} suggests further methodological development that would enable quantum simulations of the INS spectra of this endohedral complex as well.

## ACKNOWLEDGMENTS

Z.B. and P.M.F. acknowledge the National Science Foundation for its partial support of this research through Grant Nos. CHE-2054616 and CHE-2054604, respectively. P.M.F. is grateful to Professor Daniel Neuhauser for his support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.