Hydrogen confined in single-wall carbon nanotubes: Anisotropy effects on ro-vibrational quantum levels

The energy levels of a hydrogen molecule embedded in the cavity of single-walled carbon nanotubes with different morphologies are studied using quantum dynamics simulations. All degrees of freedom of the confined molecule are explicitly included in our model, revealing that the vibrational motion is notably affected by the presence of a confining potential. The most relevant effects are nevertheless found in the rotational motion of the molecule and the appearance of a quantized translational motion. We further analyze the dependence of the confinement effects on the interaction potential, considering different parameters for the carbon-hydrogen interaction.


I. INTRODUCTION
During the past ten years, the effects of nanoconfinement have generated an increasing interest both at experimental and theoretical levels. [1][2][3][4][5] Research on carbon nanotubes has developed rapidly since their discovery by Iijima in 1991, 6 due to their variety of unique mechanical, physical, and chemical properties. 7 Their capacity for encapsulating molecules and confining them in nearly uni-dimensional structures has motivated recent literature on the adsorption of gases in structures composed of a bunch of single-walled carbon nanotubes (SWCNTs). [8][9][10][11][12][13] Particular interest has been devoted to hydrogen physisorption, [14][15][16][17][18] in view of the potentiality of such nanostructures as storage devices. 19,20 The present work aims at exploring the changes induced in the spectrum of quantum ro-vibrational states of hydrogen molecules when confined in nanostructures of SWCNTs. To this aim, we model the quantum dynamics of a single hydrogen molecule embedded in a rigid carbon nanotube. This system is referred to as H 2 @SWCNT throughout the paper. The accurate simulation of the whole H 2 @SWCNT system, including explicitly all the carbon atoms of the nanotube in a fully quantum formalism, is numerically unfeasible. On the other hand, classical formalisms cannot reproduce quantumnature mechanisms that govern the processes on such smallsized systems. An alternative with a reasonable computational cost is to restrict the quantum formalism to the description of H 2 , keeping the carbon atoms of the nanotubes frozen. The positions of the carbon atoms are, however, explicitly considered at the level of the interaction.
The strategy we follow here is similar to the one used by Lu et al. 17 and Yildirim and Harris, 21 the last authors using a generic cylindrical confining potential. In our contribution, however, we propose an even more realistic model. The main novelty we introduce is the assumption of a fullquantum approach of H 2 trough the use of a five-dimensional a) Electronic mail: jaime.suarez@uam.es. b) Electronic mail: fermin.huarte@ub.edu.
Cartesian coordinate system. This allows us to simultaneously consider (1) hydrogen's internal vibration, (2) rotation of the confined molecule, and (3) its center of mass translation inside the nanotube. Up to date, theoretical studies of such systems have eluded the first point by considering a rigid hydrogen molecule, based on the assumption that the vibrational energy level spacing is too large with respect to the perturbation induced by the nanotube walls. Following this reasoning, the confinement would result in a significant coupling between rotational and translational motions without affecting the vibrational pattern. However, a serious limitation appears in the interpretation of results when comparing with experimental data. Well established optical techniques such as low-temperature infrared (IR) spectroscopy perform the analysis at the ν = 1 vibrationally excited level of hydrogen. In a recent publication concerning the encapsulation of hydrogen inside fullerenes, 22 Xu et al. proposed a basic extrapolation procedure to optimize a realistic potential energy surface for the H 2 @C 60 system. They inferred a generic vibration-less energy spectrum from vibrationally excited experimental IR data, 23 taking the effective rotational constant of the unconfined H 2 molecule. Although to a much smaller extent than for rotation, we believe that the nanoconfinement must somehow affect also the energy range and internal structure of each vibrational band.
Standard quantum molecular dynamics methods usually approach the nuclear problem by using those curvilinear coordinates which better adapt to the processes that will be reproduced, making the methodology problem specific. Regarding their computational implementation, the extrapolation to bigger systems or their adaptation to other type of processes is all but routine. At this point, the development of a Cartesian coordinates-based code which can act as a black-box has clear advantages, not the least of which are the simplicity of the algorithms (no cross-derivatives in the kinetic term) and the absence of singularities in the Hamiltonian. On the other hand, Cartesian coordinates schemes usually involve three more dimensions than curvilinear coordinates schemes, where the global rotation can be explicitly written in the Hamiltonian operator, and therefore separated off in the absence of external fields. This is not the case, though, for a hydrogen molecule embedded in a rigid and static nanotube. A reasonable model can then be achieved by considering a fivedimensional Cartesian coordinate system, with the additional advantage of an exact treatment of trans-vibro-rotational couplings. The heavy computational cost of such an approach can be overcome by the parallelization of the work over large clusters of computers. We employ the parallelized Grid time dependent Schrödinger equation (TDSE) computational tool, 24 co-authored by one of us, to solve the corresponding TDSE that rules the dynamics of the H 2 @CNT system. The structure of the paper is as follows. In Sec. II we detail the formalism that lies under the GridTDSE code. We also describe the coordinate system adopted, the selection of an appropriate initial wave function and the generation of the potential energy surface (PES) that rules the nuclear quantum dynamics. In Sec. III we give the results of the propagations. In order to facilitate the interpretation of the spectral lines of the full-dimensional 5D Hamiltonian, we first introduce the results of two-dimensional and three-dimensional models. Our main conclusions are summarized in the final section.

A. GridTDSE code
GridTDSE is a parallel Fortran 95 code that performs quantum molecular dynamics calculations by means of wave packet propagations in Cartesian coordinates. The code is based on Grid methods, which discretize the wave function ( x, t) over the whole coordinate space, with the vector x denoting the N d mass-weighted Cartesian coordinates. The TDSE that rules the dynamics is then written aŝ whereĤ and ( x, t) are, respectively, the Hamiltonian and the total wave function of the system. The computational implementation of Eq. (1) renders a non-diagonal Hamiltonian matrix, due to the non-local kinetic term in the spatial representation. The sparsity, however, can be efficiently increased by using the variable order finite difference (VOFD) method, 25,26 which considers a reduced number n s of neighboring grid points (n s is called the "stencil") to calculate the second derivatives that appear in Cartesian coordinates in the Laplacian operator. The combination of a Cartesian coordinates scheme with VOFD methods makes the parallelization of the matrix very efficient in computer memory and time. Regarding the potential energy operator, the sparsity can be additionally increased by considering an energy cut-off for the process, V c .
Considering a time-independent potential energy operator, the solution of Eq. (1) can be written in terms of an exponential evolution operator U (t,Ĥ ), which applies to the initial The time evolution operator can be then approximated by using a second order difference method (SOD), 27 which preserves norm and energy, obtaining a final two-term recursive formula. In the SOD propagation scheme, the most cumbersome operation from a computational point of view is a matrix-vector multiplication. GridTDSE uses the portable extensible toolkit for scientific computation library, 28 very efficient at carrying out parallel matrix-vector operations with sparse matrices. Further details on the code can be found in Ref. 24.
Once the integration of the TDSE yields the value of the wave function at any time, t, all molecular properties, such as power spectra and eigenfunctions, can be readily calculated from a Fourier transformation of the time-correlation function C(t), obtaining, For infinite integration time, the magnitude I(E) consists of a manifold of peaks that mimic a set of delta functions. The spectral resolution depends inversely on the total time of integration, while the spectral window decreases with the timestep. GridTDSE is also provided with an implementation of the iterative diagonalization Lanczos method. 29 In such a technique, the extraction of specific eigenvalues and eigenvectors is performed more efficiently than in standard Fourier transform methods. However, the accuracy of the method is restricted to a small-width window in the energy spectrum.
In the H 2 @CNT system, the combination of a high-energetic vibrational motion with lower energetic motions, such as rotation and translation of the center of mass, renders a broad energy range which cannot be accurately covered by this technique. Consequently, Lanczos calculations have been restricted to extract specific energies at the lowest part of the fundamental vibrational band, while the Fourier transform of C(t) has been used to obtain the general spectrum.

B. Coordinate system and initial wave function
In order to use GridTDSE to study the dynamics of H 2 embedded in a SWCNT, we have implemented a fivedimensional Cartesian coordinate Hamiltonian system. Three Cartesian coordinates (ρ x ,ρ y ,ρ z ) are devoted to characterize the internal motion of H 2 (rotation and vibration), described by the internuclear vector ρ. Two additional Cartesian coordinates (X cm ,Y cm ) refer to the confined translation of the H 2 center of mass on the plane perpendicular to the nanotube axis, associated to the vector R. Here we have made use of the assumption that the CNT is an infinite-length structure. The three types of motion previously cited (vibration, rotation, and translation) have very different characteristic energy ranges. As we will see later, the first vibrationally excited state of the confined H 2 molecule is 4145 cm −1 above the ground state, approximately one order of magnitude larger than the energies associated to rotation and translation inside the nanotube. Consequently, the density of confined states in the energy spectrum will be enormous, and it will extend over a very long energy range (thousands of wavenumbers). In such a pattern, the identification of the energy levels by graphical inspection of every single eigenstate is unfeasible.
At this point, the selection of a initial wave function with specific symmetry conditions may act as a filtering technique for clarifying the molecular spectrum. All the spectral information of the bound states gathered from Eq. (4) is related to the initial wave function 0 ( x) introduced in Eq. (2). In fact, the height of the lines correlate to the squared absolute value of the coefficients in the spectral decomposition of 0 ( x), while their positions correspond to the respective eigenenergies. The standard procedure for building the initial wave packet consists in introducing a direct product of gaussian functions, one for each quasi-uncoupled degree of freedom, exciting the specific modes by displacing the corresponding gaussian function from its equilibrium position. Alternatively, one can impose specific symmetry configurations for the global initial wave function, so that only those modes with the same symmetry will appear in the spectrum.
In an equivalent way to the procedure adopted by Lu et al. 17 for the identification of the eigenstates, we introduce a translation/(rotation-vibration) separable initial wave function, ( x) = ( R) ( ρ). The implicit assumption here is to initially decouple translation and treat the vibrationalrotational coupling as in the case of an unconfined H 2 molecule. In general terms, such a product will not be the solution of the five-dimensional equation, and therefore we should not expect a single-line spectrum, but it is a valuable graphical tool to find out the main characteristic motion in each energy level. The separable initial wave function can be employed as an even more powerful filtering tool, using the solutions of a two-dimensional ( R) or a three-dimensional ( ρ) model. In order to do so we will first study the case of a structure-less H 2 molecule confined in the interior of an infinitely long SWCNT and build the solutions of the corresponding 2D Hamiltonian. Next, the three-dimensional case of the unconfined H 2 molecule will be solved. The solutions of both models will serve as initial wave functions in the fulldimensional calculations, trying to extract selected information from the dynamics simulation. Such will be the general strategy in the presentation of results that appears in Sec. III.

C. Potential energy surface
The first step, however, concerns the generation of the PES for the model system, a diatomic homo-nuclear molecule inside an infinite-length rigid SWCNT. In view of the small fluctuations found for the potential energy when the center of mass (c.o.m) of H 2 moves along the CNT axis (here Z axis) compared to the energy variations long a perpendicular direction, we have circumscribed the motion of the c.o.m. to the Z cm = 0 plane. The explicit position of the C atoms in the nanotube is instead considered when treating the internal motion of the H 2 molecule. In other words, the potential is Z-invariant considering the center of mass ( R) but not considering the ρ coordinate.
There are typically two main possibilities when a hydrogen molecule encounters a carbon based nanostructure: It can either dissociate into hydrogen atoms and form strong C-H bonds or it can keep its molecular form and bind only weakly with the carbon structure. In the non-dissociative process of hydrogen physisorption by SWCNTs, the main responsible for the interactions are instantaneous electron correlations, which can be modeled by Van der Waals forces. At present, non-reactive process of such weakly bound systems can readily be formulated from semi-empirical atom-atom interactions. In this way, we can detach the internal potential for H 2 from the interaction of each hydrogen atom with the nanotube, resulting in a two-term potential function, The associated five-dimensional Hamiltonian describes simultaneously the rotation, vibration, and translation of H 2 inside the SWCNT. The first term of the potential in Eq. (5) accounts for the length of the H-H covalent bond, and is modeled employing a Morse potential 30 with the standard parameters 18, 31 for H 2 . The second term describes the atomatom interaction of each hydrogen with the C atoms of the rigid nanotube, and consists of a sum of C-H Lennard-Jones potential energy functions, where r C i −H j is the C i -H j distance determined by R and ρ. Most part of the calculations presented in this paper have been carried out with a set of values σ = 2.82 Å and = 0.0605 kcal/mol previously employed by Huarte-Larrañaga and Albertí to simulate hydrogen physisorption in CNTs. 18 Our molecule-nanotube interaction model has evident limitations in terms of accuracy and clearly a more accurate potential energy surface, maybe developed from the recent family of dispersion corrected density functional theory (DFT) methods, 15, 32 would mean a significant breakthrough for the accuracy of quantum dynamics simulations. However, such an advance can only be supported by broad experimental confirmation, while only a few neutron scattering experiments have been successfully carried out in the past. 1,33,34 Despite the general qualitative approach in the present work, we cannot ignore the considerable discussion on the election of the semi-empirical potential parameters in the literature. 16,17,22,35 In this respect, Lu et al. 17 showed the influence of σ and parameters in the distribution of the confined energy levels, revealing the key role that the equilibrium distance characteristic of the C-H interaction plays in quantum sieving. They proposed the use of a Franck-Brenner potential with the set of parameters σ = 3.08 Å and = 0.055 kcal/mol (labelled FB). On the other hand, a threecenter semi-empirical C-H potential was modelled by Bacić

A. Two-dimensional model: {X cm , Y cm }
Concerning the translation of the center of mass of H 2 inside the nanotube, an illustrative two-dimensional potential energy surface U 2D ( R) can be obtained as a function of the {X cm , Y cm } coordinates by minimizing the potential over all possible orientations of the inter-nuclear vector ρ. The confined two-dimensional translation of the molecule inside the carbon nanotube originates discretized energy levels E (n) , with the associated wave functions (n) ( R). In previous contributions 17, 21 it was shown that, as a general feature, the variation of the CNT radius resulted in two different types of confining potentials U 2D ( R). Narrowtype CNTs exhibit parabolic potentials centered at the nanotube axis, while for wider nanotubes the potential minimum is displaced towards the nanotube walls, exhibiting a ring structure as a result of the cylindrical symmetry. The CNTs considered in this work, with chiral indices (n = 8,m = 0) and (n = 10,m = 0), have been specifically chosen as representative examples of parabolic and ring structure type potentials, respectively. The "narrow" (8,0) nanotube is 6.26 Å in diameter, while the "broad" (10,0) one is 7.84 Å. They have both been generated concatenating 20 unit cells along the nanotube axis, assuring that the ends of the nanotube do not play any role.
The energies E n and eigenstates (n) of the 2D model were obtained from the TDSE equation using the Lanczos diagonalization implementation in the GridTDSE code. In the case of the parabolic potential created by the (8,0) SWCNT, the potential energy surface presents a minimum of −1628 cm −1 located at the axis of the nanotube ({X cm = 0, Y cm = 0}). The calculated spectrum resembles the structure of a quasi-harmonic potential, with the energy levels approximately equally separated. We superimpose in Fig. 1 (left) the energy levels over a Y cm = 0 cut of the potential. Figure 2 shows the probability density plots of the edge functions corresponding to the first four levels. The eigenenergies of the first two levels are E (n = 0) = 214.13 cm −1 and E (n = 1) = 453.44 cm −1 above the minimum of the PES, the second being doubly degenerated due to the symmetry between X and Y axis. Unlike the 2D harmonic oscillator case, the third state is only doubly degenerated E (n = 2) = 714.64 cm −1 , being the fourth (non-degenerated) value slightly more energetic (739.54 cm −1 ). As Yildirim and Harris explained in their paper concerning cylindrical-type potentials, 21 this shift is the result of the anharmonicity in the confining potential produced by the walls of the CNT.
The structure of the confining potential changes significantly in the case of the (10,0) SWCNT. As expected for a broader nanotube, a ring structure is now observed in the PES, with the potential energy minimum located 0.97 Å away from the axis of the nanotube. Only very small oscillations of the potential (±0.21 cm −1 ) along the azimuthal direction can be detected, yielding a quasi-freeĴ 2D z orbiting about the nanotube axis. In the right panel of Fig. 1, we display a one-dimensional cut of the PES, including the first energy levels found for this system. The mentioned free orbiting is at the root of the B 2D l(l + 1) progression found for the lowest energies of the two-dimensional model, l being the quantum number associated to theĴ 2D z operator. In fact, the rotational constant associated to the classical orbiting of a 2m H mass particle along the R min = 0.97 Å circumference is B 2D 0 =¯2/(2m H R 2 min ) = 8.92 a.u., coherent with the energy structure obtained for the first four states (E n = 0 = 80.12 cm −1 , E n = 1 = 99.88 cm −1 , E n = 2 = 139.9 cm −1 , and E n = 3 = 195.01 cm −1 ). As expected, the probability density of the corresponding eigenstates, shown in Fig. 3, is mainly localized over the ring structure. At the same time, the orbiting pattern combines with a parabolic anharmonic effect induced by the CNT wall when the molecule approaches to it. The coupling between both motions yields a more complicated picture for the whole energy range, as it is shown in Fig. 1.

B. Three-dimensional model: Unconfined H 2 in Cartesian coordinates
A spectral analysis of the vibrational and rotational levels of unconfined H 2 has been performed using GridTDSE in a three-dimensional Cartesian system (ρ x ,ρ y ,ρ z ). Despite the apparent simplicity of a diatomic molecular system, usually described by a Morse potential, the absence of exact analytical solutions for rotating molecules (J = 0) is still leading to new work on the subject. 36,37 In standard curvilinear coordinatesbased methods, where the vibrational-rotational coupling is explicit in the Hamiltonian, the usual procedure to solve the TDSE is to employ effective J-dependent potentials in either variational or second-order perturbation schemes. Both paths present limitations for the calculations of large-angular momentum states, since the perturbative term is then larger. On the contrary, a Cartesian coordinates formalism includes the vibrational-rotational coupling implicitly in the Laplacian of the Hamiltonian operator. In order to extract the eigenstates associated to a specific angular momentum values, projection operators are introduced at the level of the initial wave function, 38 resulting a new wave function J ( x) = P J ( x).
In contrast to general curvilinear methods, the projection procedure is exact and does not involve significant changes in the computational implementation for higher J values. Explicit analytical formulae of such projections for the diatomic case can be found in the bibliography. 39,40 The procedure for introducing the angular momentum J in the calculations consists of propagating an arbitrary wave packet-which is not eigenstate of the total angular momentum operator-under the global Cartesian Hamiltonian. Throughout the propagation, one evaluates the overlapping with the J-projection of the initial wave function J (t = 0), obtaining an autocorrelation function C J (t) = J (0)| (t) which contains all the relevant spectral information associated to the specific J value. We can, therefore, perform the spectral analysis for several J values from a single wave packet propagation. This technique was successfully employed in the past to obtain accurate energies for highly rotationally excited diatomic (H 2 ) and triatomic (H 2 O) molecules. 24 The zero point energy (ZPE) of the unconfined H 2 molecule is located at 2166.87 cm −1 , being the first rotationally excited levels 119.54 cm −1 , 357.45 cm −1 , 711.57 cm −1 , and 1178.58 cm −1 above this value. The first vibrational overtone is 4145.7 cm −1 , just above the seventh-excited rotational level. As we see in the ro-vibrational spectrum depicted in Fig. 4 (upper panel), the agreement of the J = 0 levels with the analytical solutions of a Morse potential (vertical dotted lines) is perfect. The large energy gap between vibrational levels, compared to the rotational and translational (seen in Sec. II A) ones, queries the effect that the CNT may have on the vibrational of the molecule.

C. Five-dimensional model: Full QMD
In the five-dimensional model that fully represents the quantum molecular dynamics of H 2 embedded in SWCNTs, the nanotube is taken as a rigid structure whose axis lays always along the Z direction. Under this assumption, the total angular momentum, considering the whole system nanotube + molecule, is not well defined. However, useful information can be gathered by introducing the rotational angular momentum quantum number j of the H 2 molecule, although the associated rotational operator does not exactly commute with the total molecular Hamiltonian. But it can be considered as a good quantum number in a first approximation. From now on, we will substitute J by j, referring to the internal rotation of H 2 .

Parabolic confinement CNTs: (8,0)
We first examine the effects of nanoconfinement produced by the carbon nanotube of chiral index (8,0). The in- put parameters of the GridTDSE code employed for the 5D calculations of this work can be found in Table I. The convergence of the calculations with the number of points in the grid, the number of stencil points, and the potential cut-off has been verified. The time dependent Schrödinger equation has been integrated up to 48.4 ps, with a resolution in the energy spectrum of ∼0.6 cm −1 . When referring to the energy values of specific molecular eigenstates, additional accuracy is provided by performing Lanczos calculations with a similar implementation for the Hamiltonian matrix. The power spectrum obtained from the full-dimensional simulation of the confined molecule is shown in the bottom panel of Fig. 4. The localized gaussian function that has been used as initial wave function represents the H 2 molecule with its c.o.m located 1 a.u. away from the bottom of the potential well and with the inter-nuclear vector leaning 60 • with respect to the nanotube axis. The energy range in the lower panel has been displaced by 372.6 cm −1 with respect to the unconfined system (upper panel) in order to match the first line of the spectrum of both panels. The overall structure of the spectrum in the 5D confined situation is much more complicated and  II. Energies (E in cm −1 ) of the lowest line for each vibrational band, corresponding to states with fundamental rotation and translation for both unconfined H 2 (3D calculation) and H 2 @(8,0) (5D calculation). Energy gaps between the first ortho-and para-levels ( E ortho-para ) are included in cm −1 in an adjacent column. Energy differences between consecutive vibrational overtones are displayed in parenthesis.

Unconfined 3D
Confined 5D dense due to the combination of translational and rotational lines, but we can still clearly distinguish the corresponding vibrational bands. From graphical inspection we note that as the vibrational number increases, the separation in energy between adjacent vibrational bands increases with respect to the 3D unconfined system. The effect of confinement on the vibrational structure can be better quantified in Table II, where the energies of the lowest eigenstate for each vibrational band are specified for both the unconfined (tri-dimensional) and confined (five-dimensional) models. The difference in energy between adjacent vibrational levels, displayed in parenthesis in both calculations, is typically ∼25 cm −1 larger in the fivedimensional calculations than in the three-dimensional ones. This is an overall pattern that can be explained by considering that as the H 2 vibrational excitation increases, the total wave function reaches outer regions of the phase space, where the repulsive effect of the confining nanotube walls is enhanced. One noticeable effect of the confinement induced by the nanostructure is to raise the ZPE by 372.6 cm −1 . This increase can be partially attributed to the residual energy added by the translation of the center of mass (Sec. III A). However, our calculations for a (reduced) two-dimensional model predicted a value of 213 cm −1 for the lowest energy value. The remaining 160 cm −1 gap is due to a significant coupling among translational and rotational motions, which is hindered in the two-dimensional case (we recall here that the 2D potential energy surface is obtained minimizing the energy along ρ). A first filtering of the dense spectrum previously shown in Fig. 4 can be efficiently achieved by introducing the eigenstates (n) (X cm , Y cm ) of the 2D model in the translational part of the initial wave function. Concerning the rotational part, we can further elucidate the structure by means of the H 2 angular momentum j, assuming specific projections P j (ρ x , ρ y , ρ z ) in the "internal" part of the initial wave function. Figure 5 represents the spectral lines obtained for several j values and using the lowest eigenfunction (0) (X cm , Y cm ) in the translational part of the initial wave function. The figure has been divided in three panels, corresponding to a 1600 cm −1 extract of the first three vibrational bands, with similar energy scales in order to better compare the effects of the confinement. Should j be a constant of motion, each spectral line would be uniquely correlated to the diatomic rotational quantum number j, or in other words, to the corresponding j-projection of the internal part of the initial wave function (labelled as j ). In our model, however, the rotation of the nanotube is strictly forbidden, and consequently the operator j cannot commute with the total Hamiltonian of the system. This causes multiple contributions for a specific spectrum line, the contributions stemming from different j-projected initial wave functions. Due to the homo-nuclear symmetry of the H 2 molecule, the mixing occurs only between lines with the same j-parity 17 (e.g., the curves associated to j = 1 and j = 3, or those associated to j = 0 and j = 2). For most of the lines, however, there is a dominant contribution, and we can identify five different rotational bands in each panel.
Concerning the composition of each rotational band, we see how the (2j + 1)-degeneracy of the unconfined H 2 system (coarse blue dotted lines), characteristic of an isotropic potential, breaks under the cylindrically symmetric external confining potential into j + 1 different spectral lines. The identification of the axis of the nanotube (Z in our Cartesian coordinate system) with the axis of the azimuthal momentum operator, represented by the azimuthal quantum number m, is key for understanding this behavior. In fact, we can refer only to its absolute value |m|, since the symmetry of the nanotube indicates that +m and −m values correspond to states with the same energy (hence the appearance of j + 1 different energy values). Additional valuable information can be gathered from varying the orientation of the inter-nuclear vector in the initial wave function. Figure 6 shows the spectral lines obtained using three initial gaussian wave packets (g) differently oriented with respect to the nanotube axis, namely, 0 o (along Z axis), 90 o (along X axis), and 60 • (the one used in the previous two figures). Contrary to the oblique orientation (60 • ), ρ x -and ρ z -oriented initial wavefunctions keep in part the symmetry of the potential and originate filtered spectra. The reason for this can be viewed adopting the classical counterpart. In this formalism, the angular momentum vector aligns in a plane perpendicular to the initial orientation of the molecule. A molecule initially oriented along the Z axis can only contribute to m = 0 states. That explains the single line obtained at the lowest part of each rotational band when a Z-aligned initial wave function is employed. In the case of X-alignment, the molecule could rotate in any plane which contains the X axis. The gaussian wave function introduced is in fact a combination of functions with different |m| = 0 values associated. Such a combination includes the |m| = j value and those states with the same symmetry, which is determined by the parity of j − |m|. As a result of it, the only spectral lines that are not filtered out now are associated to |m| = j, j − 2, . . . values. We can therefore conclude that high |m| values are associated with wave functions that spread in the XY plane, while lower ones reflect a localization along the Z axis. With the center of mass of the molecule mainly localized at the nanotube axis as corresponds to the fundamental translational motion ( (0) ), the proximity between hydrogen atoms and the carbon walls enhances the confining energy 2600 2800 3000 3200 3400 3600 3800  for rotational states associated to high-|m| values, increasing their energy with respect to the unconfined case. In the low-|m| case, the hydrogen atoms tend to be localized along the Z axis, where the confining potential is minimal and negative (see Fig. 1), and the energies are even lower than in the isolated H 2 molecule. The energy range covered by each rotational band is then mainly due to the orientation of the total angular momentum with respect to the nanotube axis, and thus does not change significantly from one band to another. 2600 2800 3000 3200 3400 3600 3800 FIG. 6. Same as Fig. 5 including three different orientations of H 2 internuclear vector: ρ = (0.866, 1.5, 1) a.u., ρ = (0, 0, 2) a.u. and ρ = (2, 0, 0) a.u.
As one moves to higher vibrational levels (lower panels), each j-manifold becomes broader. Two competitive effects are responsible for this: On the one hand, the vibrationalrotational coupling equally red-shifts all the eigenstates of a certain j-band. It is the consequence of a higher effective momentum of inertia, and thus a lower rotational constant. On the other hand, in a similar way as we explained before, the combination of high-|m| values and a longer vibrational amplitude implies a shorter distance of the hydrogen atoms to the nanotube walls. This consequence blue-shifts those spectral lines associated with highest energies inside each rotational band. For low-|m| values, however, the vibration takes place mainly along the nanotube axis, and the magnitude of the confining energy is mainly unaffected throughout the vibration. To sum up, the vibrational-rotational coupling pulls down the energy of all the levels in the manifold, while the confining potential lifts only the most energetic ones.
This pattern is also at the root of the decreasing gap between the first two spectral lines in each vibrational band as ν increases. The two levels correspond in fact to the lowest state of para-(characterized by even-j values) and ortho-hydrogen (odd-j values), respectively. Treating the confinement as a perturbation U , the internal part of the ground (para-) state correlates in the U − → 0 limit to the spherical harmonic Y 0 0 , which is a hollow sphere centered at ρ 0 = 1.4 a.u. The presence of the nanotube constrains elliptically the unperturbed orbital, originating the mixing between {even-j,|m = 0|} levels. In such a quasi-spherical structure, more quanta on the vibrational spectrum significantly affects the thickness of the orbital. The wave function spreads now on a part of the phase space which lays closer to the nanotube wall, and thus the confining energy increases. On the other hand, the first-excited (ortho-) rotational state correlates in the U − → 0 limit to a p z -type orbital. Upon vibrational excitation, such an orbital elongates along the Z axis, along which the confining energy U ( ρ, R) stays constant at its minimum (negative) value. As a result, the energies of the first two states then approach each other, and the structure of both orbitals differ in the presence of a node in Z = 0 for the antisymmetric first-excited state. This claim is evidenced in the E ortho-para column of Table II. Following the mentioned column, the rotational energy gap in the confined case drops from 60 cm −1 for ν = 0 to 12 cm −1 for ν = 4, therefore, at a faster rate (∼12 cm −1 between consecutive vibrational levels) than in the unconfined case (∼6 cm −1 ). Due to the effect of confinement inside the nanotube, lowest ortho-and para-hydrogen states reach almost similar energy values in the ν = 5 vibrational band. This leads us to expect an interchange in the order of ortho-and para-levels for higher vibrational bands, but this could not be assured since the computation of such highly excited states demands a much more dense and extended grid than the one used in our fivedimensional calculations.
A similar effect appears in the case of excitations on the translation of the H 2 center of mass. A very efficient way to analyze this effect consists of introducing in the translational part of the initial wave function the first ( (n = 1) ), second ( (n = 2) ), and third ( (n = 3) ) eigenstates of the 2D Hamiltonian. The excitation introduced in the translational part of the wave function increasingly shifts the lowest ro-vibrational (internal) state towards higher energies. Table III lists the energy gap between the first three excitations in the translational motion and the n = 0 ground state, for given vibrational quantum numbers. The first row of the column contains the corresponding energy differences in the bi-dimensional model, as a reference. The gap is notably lower in the 2D model than the one we obtained in the ν = 0 vibrational band of the 5D case. Even for the lowest vibrational band, the difference between two-and five-(full-)dimensional models is ∼45 cm −1 for the eigenstate associated to n = 1, increasing to ∼90 cm −1 for the next two levels, n = 2 and n = 3, which would TABLE III. Energy gap (cm −1 ) between states with n-excited and (n = 0)fundamental translational motion and associated to rotational labels (j, m) = (0, 0) are displayed left to right for the H 2 @(8,0) system. Each row corresponds to a different ν-vibrational band. The energy differences found for the reduced 2D model (Sec. III A) are included in the first row as a reference. have similar energies for harmonic confining potentials. We recall here that the potential energy in the two-dimensional model was minimized with respect to all possible orientations of the inter-nuclear vector ρ. In the five-dimensional model, however, the translation does not take place under this minimized PES, exhibiting higher energy values than in the twodimensional model. Obviously, the coupling is more effective for n = 2 and n = 3 than for n = 1 since the wave function spreads in a broader region of the subspace {X cm , Y cm }. From Table III it is possible to analyze also how the vibrational pattern is affected by successive excitations in the translational part. The energy gap E ν − E ν − 1 uniformly increases with ν along a specific n-column, but in merely a few wavenumbers. Instead, for the ν = 4 vibrational band, the shift suddenly becomes much more significant, and ranges from ∼40 cm −1 for the n = 1 case to ∼70 cm −1 for the spectrum obtained with n = 2 and n = 3 initial states. At this stage, it seems that the vibrational amplitude of the diatom is elongated enough (∼2.5 a.u.) as to produce, combined with the displacement of the center of mass of H 2 , a significant increase of the repulsive energy due to the short distance of the H atoms to to the nanotube walls.
We focus now on the effect that excitations on the translational part induce in the rotational pattern. The three panels of Fig. 7 show the ν = 0 vibrational band of the power spectra obtained from initial states of the type j (n) , with n ranging from n = 1 to n = 3. Here again, P j projections are taken for the internal part of the wave function j . As in the case of n = 0, shown in Fig. 5, we can still distinguish nonoverlapping rotational bands. However, the structure of each one is more complex now, with new lines stemming from the highest m-levels of each band as a consequence of the new symmetry in the translational part of the wave function. We note also that the lowest state, corresponding to m = 0, is basically unperturbed in all cases. We can better understand the new rotational pattern by inspection of the corresponding 2D eigenstates shown in Fig. 2. The ground state (n = 0) is a quasi-gaussian function centered at the axis of the nanotube. Because of the cylindrical symmetry of this function there are two equivalent orientations of the angular momentum (j x and j y ) and this yields a spectrum with only (j + 1) lines per band, being those states with j = 0 doubly degenerated. For the case n = 1, the two-dimensional eigenstate has the shape of a p xtype orbital (or p y , due to the degeneracy), with one single node along the X-(Y-)direction. The asymmetry between Xand Y-rotations in this case unveils the j + 1-degeneracy and originates new lines in the spectrum. A simple classical interpretation in this case is not enough to justify the unfolding of degeneracy and a fully quantum view is required.  isotropic potential proposed by Yildirim and Harris the lines are exactly equidistant). The number of spectral lines that appear in the j = 2 and j = 3 bands increases to six and eight, respectively. In a similar "parabolic" potential obtained from a (3,6)-CNT, Lu et al. could identify seven different lines for the j = 2 band. However, two of these states where so close in energy that it is very feasible that they would conform one single line under the specific potential parameters chosen in our work. The results obtained in the second panel when the translational part of the initial wave function is replaced by the second-excited eigenstate (n = 2) show again an increase in the number of states associated to each j-band. This is not surprising provided the four-lobbed structure of the translational part of the wave function (third panel of Fig. 2). In particular, the rotational bands j = 1, 2, 3 consist now mainly of five, six, and seven spectral lines. As it has been pointed out earlier in Sec. III A, the energies of n = 2 and n = 3 translational states in the 2D model are very similar, and the difference between them is only due to the anharmonic character of the confining potential. Instead, the spectra they generate are remarkably different. A notable simplification in the 5D spectral curves is observed for the case n = 3, where the wave function regains the cylindrical symmetry in the translational part of the wave function, as depicted in the fourth panel of Fig. 2. The rotational bands j = 1, 2, 3 consist in this case of three, four, and six peaks, respectively. We must mention here that, although both n = 0 and n = 3 2D states are cylindrically symmetric, the second presents also a node in the radial direction of the XY plane which makes its spectral pattern rather different to the n = 0 case. We finally note that the energy gap between the lowest ortho-and pararotational levels in the power spectrum (labelled j = 0, |m| = 0, black line, and j = 1, |m| = 0, first red line) is notably reduced in the highest translational levels. The reasoning that supports this behavior was already exposed when commenting Fig. 6. Similar to vibrational excitations, more quanta on the translational state implies shorter distances to the nanotube walls, and therefore stronger interactions, pulling j = 0 and j = 1, |m| = 0 together.
Previous works showed that the separation between the lowest levels assigned to j = 1 and j = 0 (ortho-and paralevels) can be related to the intensity of quantum sieving effects. 16,17 We have used GridTDSE for a better visualization of the dependency of the spectrum on the choice of parameters , σ used for the modelization of the molecule-wall interaction. The results are displayed in Fig. 8 for two different potential parameterizations. The upper panel represents spectral curves obtained for the potential "FB" of Ref. 17. This choice of parameters, with a longer C-H equilibrium distance, originates extreme 1D confinement, resulting in strong quantum sieving. In our calculations, this reveals in a dramatic spread of the levels included in each j-manifold, with the consequence of a quasi-degeneracy between the first two levels (correlated to j = 0, |m| = 0 and j = 1, |m = 0| states of the unperturbed system). We also note stronger mixing between different j contributions for each line, and consequently j cannot be considered as a good quantum number in such a system. The reference lines from the original paper of Lu et al. 17 have been blue-shifted 78.3 cm −1 in order to match our fundamental state. This might be the reflection of the absence of H 2 vibrational motion in their model. Still, we find small 3000 3200 3400 3600 3800 4000 4200 4400  22 in the study of H 2 embedded in C 60 . The nanoconfinement does not distort significantly the molecular energy structure of the diatom, as it did for the FB potential and, to a lesser extent, for the potential parameters used in the previous calculations. Although the breaking of j-degeneracy still appears in clear and distinguishable rotational bands, there is no significant mixing among different j-contributions for each molecular eigenstate. In this case j can be considered as a good quantum number, and the interaction of the molecule with the wall can be treated as a first-order perturbation of the unconfined situation. The notable discrepancy with the other two potentials might be a consequence of optimizing C-H interactions from a spherical symmetry system (fullerene) instead of a cylindrical structure (nanotube), using a three-center potential for each pair of molecule-carbon interaction.

Ring confinement CNTs: (10,0)
The next example of nanoconfinement concerns the embedding of hydrogen in a wider carbon nanotube of chiral indices (10,0). The GridTDSE code was modified to implement the corresponding potential energy function, which in this case presents a ring structure for the minimum of the (minimized) confining potential U 2D ( R) as it can be seen in Fig. 1 (right), characteristic of wider nanotubes. 21 Compared to the case of H 2 @(8,0), the energy spectrum is expected to be considerably more complicated due to the quasi-free orbiting associated to this ring structure. This type of angular translational motion on the {X cm , Y cm } plane can be represented by a L z operator, which will efficiently couple those rotational states of H 2 with non-zero azimuthal number (m = 0). Most of the structures in these spectra can be assigned FIG. 9. Power spectra for H 2 @(10,0) along the first three vibrational bands. The translational part of the initial wave function is a gaussian centered at R cm = (1, 0) a.u. For the "internal" part of the wave function we have taken j-projections of a gaussian initially oriented on ρ = (0.866, 1.5, 1) a.u. Vertical dotted lines correspond to ro-vibrational levels of the unconfined H 2 molecule labelled as (ν, j) and translated to the ground state energy of the H 2 -CNT system.
to one-dimensional confinement along the circle of minimum energy. We start by displaying in Fig. 9 the power spectra for the first three vibrational bands, obtained from projections in the angular momentum j of H 2 and placing the initial translational wave packet away from the minimum energy zone. In this manner, most of the excitations referred to the translational part will be reproduced. As a general feature, we confirm that the spectral density is increased with respect to the (8,0) case. At the same time, the energy shift of the ground state with respect to the unconfined system is 95.7 cm −1 , approximately 275 cm −1 smaller than for the (8,0) SWCNT system due to the wider confining potential well. The structure of the spectra is conserved in general terms as we move up to excited vibrational bands. Nevertheless, there is a very subtle increase in the energy gap between the lowest line and the first-and second-excited lines as we move up in the vibrational level, i.e., going from the top panel to the bottom panel. This change is hardly noticeable in the figure, with energy gaps of +2 cm −1 and +6 cm −1 for ν = 1 and +11 cm −1 and +18 cm −1 for ν = 2, both taken with respect to the ν = 0 band. The reason for this increase is that the region of minimum energy is located close to the walls of the nanotube and a slight increase in the inter-nuclear distance already results in higher confining energies. At the same time, the rotational constant B associated to j = 1 levels decreases compared to j = 0 levels as a result of the vibrational-rotational coupling. Since the third and fourth states of each vibrational band are characterized by j = 0 and j = 1 numbers, respectively, the energy gap that we find between them decreases as we move up on ν.
In the wide (10,0) SWCNT, the nanoconfining effect can be treated as a perturbation in very good approximation. As a consequence, the mixing between different j-contributions found in the spectrum (being j the rotational quantum number defined in the isolated H 2 molecule) is almost inexistent now, and we can safely consider j as a good quantum number. Combining the j-labelling with the identification of the translational part according to the 2D model introduced in Sec. III A, one can easily assign the first energy levels as displayed in Table IV. The energies associated to the translation of the center of mass are now typically smaller than those associated to the rotational motion. In Fig. 10, we focus on the ν = 0 vibrational band of the spectra obtained from three different initial wave functions, whose translational part is built from the 2D eigenstates (n = 0, 1, 2) . The first two excited states correspond to the excitations in the orbital translation of the H 2 c.o.m. around the axis of the nanotube, which explains the l(l + 1) TABLE IV. Energies (cm −1 ) associated to the first three translational levels n = 0, 1, 2, aligned according to the rotational label j, are displayed for the fundamental vibrational band of H 2 @(10,0) nanotubes. pattern found for the energies of such eigenstates, and already pointed out in Sec. III A. The spreading inside each rotational band is much smaller than it was for the H 2 @(8,0) system. In fact, since the perturbation is now smaller, the degeneracy between rotational energy levels of same j value is increased. As a matter of fact, the j = 1, 2, 3 and 4 rotational bands associated to the ground translational state (n = 0) consist only of two spectral lines, hardly noticeable in Fig. 10 (consider the 0.5 cm −1 precision of the spectra obtained), in contrast to the j + 1 lines that we obtained for the (8,0) nanotube. Although this might appear contradictory to the dense spectra that we obtained for H 2 @(10,0), the degeneracy is counter balanced by the more complicated translational structures (orbiting + vibration inside the nanotube). As soon as the translational part includes phonon excitations (n = 0) the number of spectral lines is amplified (except for the j = 0 case, which consists always of one single line).

IV. CONCLUSIONS
We have simulated the quantum dynamics of a hydrogen molecule confined in the interior of two CNTs with different diameters. The narrower nanotube involves a parabolic-type potential, while the potential energy surface of the wider one is characterized by a ring-structure. Unlike previous studies on the subject, we have carried out a full-dimensional treatment, as far as the confined molecule is concerned, based on a Cartesian coordinate system. Such a model explicitly includes the vibrational motion of the molecule in addition to the rotational and translational ones already considered in the bibliography. Furthermore, we have taken into account explicitly the positions of the carbon atoms in the nanotube for building the potential energy surface.
The results obtained confirm previous observations on the effect that nanoconfinement has on rotational levels. On the one hand, the confinement breaks the j-degeneracy characteristic of the unconfined molecule, giving rise to a structure of broad rotational bands. On the other hand, it causes the appearance of new states whose quantum-mechanical nature is associated to the motion of the center of mass of the H 2 molecule. We have additionally evidenced that the vibration is also affected by the confinement of the nanotube. Indeed, the energy levels found for excited vibrational bands are blue-shifted with respect to the pattern found for the ground vibrational state. It is also remarkable the uniformly decreasing gap between the lowest ortho-and para-energy levels as the vibrational excitation increases, which makes one wonder about a possible exchange in the order of appearance of these levels for higher vibrational bands. Both effects cannot be neglected if one intends to generate an accurate potential energy surface for the H 2 @SWCNT system, especially in the case of the narrower nanotubes.
With respect to the morphology of the nanotube, we have shown that in the wider nanotubes the effective coupling between the translational and the internal degrees of freedom (rotation and vibration) is much smaller than in the narrower nanotubes, producing less distortion in the energy spectrum with respect to the ro-vibrational spectrum of the unconfined molecule. Nevertheless, the presence of a ring-structured minimum in the PES of the (10,0) SWCNT generates a manifold of (orbital) translational states that enhances the density of the spectra obtained for the confined molecule. In such a case, the strongest coupling occurs between the rotational motion of the molecule and the orbital translation of the center of mass about the nanotube axis.
Finally, we have studied the role that the C-H interaction potential plays on the energy spectrum pattern. As expected, the most relevant modifications in the energy spectrum are obtained for the C-H interaction potential with a larger equilibrium distance (the FB potential), yielding an extreme 1D confinement. A milder confining potential such as the one proposed by Bacić and co-workers, produces only slight modifications in the rotational spectrum. The discrepancies are expected to be larger for vibrational excited levels.

ACKNOWLEDGMENTS
This work has been supported by the Spanish Ministerio de Ciencia e Innovación (MICINN) (Project No. CTQ2009-12215) and the Generalitat de Catalunya Agencia de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) (2009SGR17 and XRQTC). Centre de Serveis Científics i Acadèmics de Catalunya (CESCA) is thanked for allocating massive parallel calculations on their supercomputing facilities. The authors indebted to Professor Stavros Farantos and Dr. Stamatis Stamatiadis co-authors of the GridTDSE and who have helped in the implementation of the numerical code.