Quantum dynamics of H 2 in a carbon nanotube : Separation of time scales and resonance enhanced tunneling

Quantum confinement effects are known to affect the behavior of molecules adsorbed in nanostructured materials. In order to study these effects on the transport of a single molecule through a nanotube, we present a quantum dynamics study on the diffusion of H2 in a narrow (8,0) carbon nanotube in the low pressure limit. Transmission coefficients for the elementary step of the transport process are calculated using the flux correlation function approach and diffusion rates are obtained using the single hopping model. The different time scales associated with the motion in the confined coordinates and the motion along the nanotube's axis are utilized to develop an efficient and numerically exact approach, in which a diabatic basis describing the fast motion in the confined coordinate is employed. Furthermore, an adiabatic approximation separating the dynamics of confined and unbound coordinates is studied. The results obtained within the adiabatic approximation agree almost perfectly with the numerically exact ones. The approaches allow us to accurately study the system's dynamics on the picosecond time scale and resolve resonance structures present in the transmission coefficients. Resonance enhanced tunneling is found to be the dominant transport mechanism at low energies. Comparison with results obtained using transition state theory shows that tunneling significantly increases the diffusion rate at T < 120 K.


I. INTRODUCTION
Nanostructured materials stand as one driving force in current research, both for their interesting fundamental properties and for their potential technological applications. 1,2Among the different phenomena appearing in the nanoscale, the structural and dynamical changes experienced by molecules adsorbed in nanometric cavities, known as quantum confinement effects, are of particular interest since they allow to store light gases, achieve chemical and isotopical separation of gaseous mixtures, 3,4 and even tailor the reactivity of chemical species using specific materials. 5The understanding of these effects is of key importance when trying to design such devices.
Since the first prediction of quantum confinement by Beenakker, 6 there has been intensive research both theoretically and experimentally on systems expected to present these effects.As a result of this, the available knowledge on these systems has increased steadily in the last years.][25][26][27][28] One particular feature present in single-walled carbon nanotubes (SWCNTs) that makes them especially interesting materials is their cylindrical shape, which gives rise to the coexistence of a 2D confinement in the space perpendicular to the nanotube's axis and unbound motion along this latter coordinate.Even though 2D confinement has been relatively studied and reviewed in the literature, 12,16,29,30 only few studies have treated both the confined and the quasi-free coordinates in a fully quantum formalism.To the best of our knowledge, only Skouteris and Laganá, 31 treating the motion of an OH radical along a (10,0) SWCNT, and some of the present authors, 32,33 dealing with the H 2 molecule in a (8,0) SWCNT, have performed such studies.The most relevant conclusion drawn from these two studies is that the quasi-free coordinate of the diatomic molecule is only weakly coupled to the set of confined coordinates.This was shown quantitatively in the case of H 2 32 through the analysis of the convolution functions of a set of 6D propagated wave packets and a set of 5D eigenstates calculated at different points along the nanotube's axis.In that study, it was seen that the state mixing between different initial states was small as the diatomic molecule progressed along the nanotube's axis.This low mixing, together with different characteristic times in the sets of confined and unbound coordinates, suggests that the fast motions of the confined degrees of freedom (DOFs) could be adiabatically separated from the slow large amplitude motions in the quasi-free coordinate.
The present work exploits the idea of separation of time scales in different ways.First, a numerically exact scheme is introduced which employs the separation of time scales to reduce the numerical effort.In this scheme, first "diabatic" basis functions describing the dynamics in the five confined coordinates are constructed.In the wave packet dynamics calculation, the total wave function is then expanded in this discrete basis and only the motion in the unbound coordinate is described using a grid representation.Second, an adiabatic approximation which separates the dynamics in the unbound and the constrained coordinates is investigated.
The increased numerical efficiency offered by both approaches facilitates a detailed investigation of the dynamics on a longer time scale not considered by previous work by some of the authors, 32 in which the diffusion rate of H 2 in a (8,0) SWCNT was computed using a full 6D quantum dynamics formalism for the first time.Interestingly, the energy resolution provided by these long time simulations allows one to identify resonances corresponding to metastable states residing in the shallow potential wells along the unbound coordinate.These resonances are found to crucially enhance the transport through the carbon nanotube at low energies.
This paper is organized as follows.In Sec.II, the theoretical modeling of the H 2 @SWCNT system is reviewed first.Then a diabatic representation of the system's Hamiltonian, which allows one to efficiently compute numerically exact results, is introduced (Sec.II B) and an adiabatic approximation separating the confined coordinates and the unbound coordinate describing motion along the nanotube's axis is discussed (Sec.II C).The description of the methodology used to simulate the wave packet dynamics and to compute the cumulative reaction probabilities (CRPs) (Sec.II E) completes Sec.II.Finally, the resulting flux correlation functions and cumulative reaction probabilities are discussed and analyzed and used to compute the diffusion rate of H 2 along the nanotube (Sec.III) and the main conclusions of the work are drawn (Sec.IV).

A. The H 2 @SWCNT system
Our system consists of a single H 2 molecule trapped inside the hollow cavity of a single-walled carbon nanotube (SWCNT).The nanotube, aligned along the z axis, is treated within the frozen structure approximation, i.e., the positions of the carbon atoms are fixed and vibrations are not allowed.The geometry of the CNT is obtained from quantum chemistry calculations using the Crystal09 34 software with a B3LYP functional, yielding a unit cell with 8.1 bohrs length and 12 bohrs diameter.On the other hand, we treat the hydrogen molecule inside the nanotube as a six-dimensional system, taking into account all the molecular degrees of freedom (DOFs): vibration of the H-H bond (ρ), orientation with respect to the nanotube's axis (θ), an azimuthal angle of rotation (φ), and translation of the center of mass (c.o.m.) of the molecule in Cartesian coordinates (x, y, and z).This coordinate system is depicted in Fig. 1.Note that, since a spherical coordinate system is used for ρ, θ, and φ, the volume element related to this coordinate system, dV, is dV = ρ 2 sin θdθdφd ρdxdydz.
The full Hamiltonian, Ĥ6D , thus reads (Note that atomic units are used throughout this paper, and therefore = 1.) With the idea in mind of separating the motion of the confined and the unbound coordinates, Eq. ( 1) can be rewritten in a compact form, taking advantage of the separability of the Kinetic Energy Operator (KEO), with Tz being the z-dependent KEO term and Tq gathering the remaining KEO terms related with the set of confined coordinates, ρ, θ, φ, x, and y.For the sake of clarity, we will refer collectively to the set of confined coordinates as q.
The potential energy term, on the other hand, is not separable and couples all the DOFs.Following previous studies, the overall H 2 -SWCNT potential is divided into a Morse term, to account for the H-H bond, and a sum of Lennard-Jones pair interactions that model the van der Waals interactions between C and H, The set of parameters used in V LJ i,j correspond to the Frankland-Brenner (FB) potential 9 with σ = 5.82 bohrs and = 0.0549 hartree.This function has been used previously to define H 2 adsorption and confinement in carbon nanotubes by several authors. 16,19,33Standard parameters are used in the Morse function V H-H : D e = 0.1746 hartree, a = 1.0271 bohr 1 , and R e = 1.4 bohrs.

B. Diabatic representation of the Hamiltonian
In order to introduce a diabatic representation of our system, we start by defining a basis of states {z k } approximately localized at the points z k to represent the wave function along the z coordinate.Applying the resolution of the identity generated by this basis set onto the last two terms of Eq. ( 2), we obtain where we have introduced a reduced 5D Hamiltonian, Ĥ5D (z), which parametrically depends on the value of z and operates on the five coordinates q describing bound motion.We can define an eigenvalue equation for this operator at a given point z k and obtain a set of eigenstates ξ j (q; z k ) and eigenvalues, These ξ j (q; z k ) functions form a complete basis set in the reduced five-dimensional space.If the eigenstates of H 5D (q; z) at a fixed point z 0 , ξ j (q; z 0 ), are used as a basis to described the wave function's dependence on q, a diabatic representation is obtained, Ψ(q, z, t) = j ψ j (z, t)ξ j (q; z 0 ).
Here the vector of z-dependent functions ψ j (z, t) represents the motion along the quasi-free coordinate.Since the basis functions ξ z 0 j = ξ j (q; z 0 ) are z independent, the corresponding matrix representation of the Hamiltonian operator reads It is useful to take a closer look at Eq. ( 8) in order to understand some of its properties.First of all, note that we have not included any additional approximation to go from the general Hamiltonian, Eq. ( 1), to its diabatic form, Eq. ( 5).Therefore, Eq. ( 8) is exact in the complete basis set limit (N z → ∞).The coupling between the confined and unbound coordinates is contained in the 5D Hamiltonian matrix.Since z = z 0 has been used as reference point to obtain the diabatic basis, the matrix representation of the 5D Hamiltonian is diagonal at z = z 0 but non-diagonal elsewhere.However, due to the low coupling between the confined (q) and unbound (z) sets of degrees of freedom, it is expected that the non-diagonal elements of the matrix will be rather small, and thus a relatively low number of ξ z 0 i basis functions will be enough to yield numerically accurate results.

C. Adiabatic representation of the Hamiltonian
Although the diabatic approach just presented is rigorous and numerically exact, its physical interpretation neither is straightforward nor offers us a clear picture about the separability of the confined and unbound degrees of freedom beyond the number of basis functions needed to achieve convergence.The adiabatic approximation, on the other hand, can shed light on the physics of the system while, under certain circumstances, providing a reasonably accurate description of the system's dynamics.
An adiabatic representation is obtained by employing the z-dependent eigenstates of H 5D (q; z) as a basis for the expansion of the wave function, The corresponding matrix representation of the Hamiltonian operator contains derivative couplings due to the action of the kinetic energy operator Tz on the z-dependent basis functions ξ j (q; z).To avoid the complications resulting from derivative coupling terms, numerically accurate wave packet calculations typically avoid the use of adiabatic representations.
If the adiabatic approximation is invoked and nonadiabatic transitions between the different adiabatic states ξ j (q; z) are neglected, the motion of the one-dimensional ψj (z, t) can be given (approximately) by where the one-dimensional adiabatic Hamiltonian Ĥ(ad) The above equations provide a simple physical picture: the dynamics in each confined eigenstates, ξ j (q; z) , evolves independently on a confined eigenstate potential energy surface (cePES), ε j (z) onto which a quasi-particle with a single DOF, ψ(z, t), evolves.Furthermore, it should be noted that the adiabatic dynamics described by these equations can be simulated with negligible numerical effort once the cePESs have been determined.The value of cePESs at all grid points z k can efficiently be computed by diagonalizing the 5D Hamiltonian in the diabatic representation at these grid points, The unitary transformation matrix U(z k ) in the above equation defines the transformation from the diabatic to the adiabatic representation.The z-dependent wave packets in diabatic representation, ψ j (z, t), and the z-dependent wave packets ψj (z, t) obtained from a rigorous simulation in the adiabatic basis, which can deviate from the wave packets calculated within the adiabatic approximation, are related via It should be noted that the only difference in the implementation of both approaches is whether the Hamiltonian matrix computed at the different grid points along z is diagonalized (adiabatic) or not (diabatic).Thus, the scheme to carry out a calculation in our system using both approaches consists of three steps: 1. Calculation of the 5D eigenstates at a chosen point along the z dimension, z = z 0 .2. Evaluation (and diagonalization, in the case of the adiabatic approach) of the reduced 5D Hamiltonian operator at the different values of the z coordinate grid, using the basis previously obtained.After this step, we obtain a set of 5D-Hamiltonian matrices (diabatic) or cePESs (adiabatic) at all the grid points z k .3. Propagation of the wave packets using the Hamiltonian matrices or the confined eigenstates potential energy surfaces.
Both the 5D eigenstates and the Hamiltonian matrices (or cePES) can be stored so that steps 1 and 2 need only to be performed just once at the beginning of the study.The stored eigenstates can be used as long as Ĥ5D remains unchanged (i.e., neither the electronic potential energy nor the KEO is modified) and the grid in the unbound coordinate is maintained.The storage of the 5D eigenfunctions is particularly important, since this step is usually the most computationally demanding.

D. Time dependent quantum dynamics: Diffusion rates
Diffusion and adsorption rates stand among the most relevant observables in nanoconfined substance studies, due to the interest in the development of efficient storage devices.In the present work, we will focus on obtaining diffusion rates for the H 2 molecule inside the hollow cavity of a narrow (8,0) nanotube.Following the work of Zhang and Light, 35 diffusion rates are obtained directly from the transition rate between two potential minima through the single hopping approximation, k hop (T ), where l is the distance between adjacent potential minima and d is the dimension of the system, which equals 1 in the present work.This equation implies that thermalization due to coupling to degrees of freedom not explicitly considered in the dynamical simulations is fast compared to the hopping rate k hop (T ).Furthermore, the hopping rate will be approximated by the rate of transition through a dividing surface separating two adjacent minima.
The transition rate through a dividing surface is calculated by the thermal averaging of the corresponding energy dependent Cumulative Reaction Probability (CRP), N(E), where β = 1/k B T and k B is the Boltzmann constant.The partition function, Q(T ), has been computed as the trace of the 5D Boltzmann operator of the system, e β Ĥ5D , and then multiplied by the contribution of the unbound coordinate, z, obtained through the semiclassical expression for a particle in a periodic potential, L mT (here, L stands for the length of the unit cell and m stands for the mass of the diatom).The CRP is the sum of all (energy-dependent) transition coefficients corresponding to open channels describing passage through the barrier separating two adjacent potential minima.
The cumulative reaction probability will be computed following the flux correlation function approach [36][37][38] employing the scheme described in Refs.39 and 40.As a starting point, one defines the flux operator, F = i[ Ĥ, h], where h is a Heaviside dividing surface discriminating reactant and product geometries, and the thermal flux operator, F T = e −β Ĥ/2 Fe −β Ĥ/2 .In the first step of the calculation, the eigenvalues f n and eigenstates |f n of the thermal flux operator, for a conveniently chosen reference temperature T 0 are computed.Then N(E) is calculated via a Fourier transform of the correlation function obtained from the time propagation of the thermal flux eigenstates, It should be pointed out that when the dividing surface is placed at the top of the potential barrier separating both configurations, a straightforward interpretation of the nature of the flux eigenstates can be found: they correspond to the vibrational states of the activated complex of the reaction. 41,42f the system's dynamics shows long-living resonances, the flux correlation function f n |e −i Ĥt |f n converges only very slowly.To avoid excessively large propagation times, then it is preferable to compute the cumulative reaction probability with limited energy resolution.The limited resolution CRP N(E) is defined as the convolution of N(E) with a Gaussian function of width ∆E, N(E) can directly be computed using the flux-correlation function present in Eq. (17), Here an additional damping factor appears in the Fourier integral which determines the required propagation time.If the energy resolution ∆E is significantly smaller than k B T for all temperatures considered, one can substitute N(E) by N(E) in Eq. ( 15) without causing significant errors.

E. The Multiconfigurational Time-dependent Hartree (MCTDH) approach
The Multiconfigurational Time-dependent Hartree (MCTDH) approach 43 has been used to represent the wave function in the different quantum dynamical calculations reported in this work.The efficiency of this method lies to a great extent on the structure of the ansatz.In the MCTDH approach, the f -dimensional wave function, Ψ(Q 1 , . . ., Q f , t), is represented as a sum of Hartree products of time-dependent basis functions called Single Particle Functions (SPFs), where Q i represent the DOFs of the system.The SPFs are in turn represented in a time-independent or primitive basis set, typically consisting of Discrete Variable Representation (DVR) or Fast Fourier Transform (FFT) functions, One can straightforwardly see that, due to the time-dependent nature of the SPF basis set, the number of basis functions, f k=1 n k , required in Eq. ( 20) to achieve numerical convergence will be significantly lower than f k=1 N k , the corresponding number of basis functions used in a standard wave packet propagation scheme.The efficient implementation of the MCTDH approach requires a specific representation of the Hamiltonian.Preferably, the Hamiltonian is written as a sum of products of operators acting only on a single Q k . 43,44hile kinetic energy operator can almost always be written as a sum of products of one-dimensional operators, multidimensional PESs usually are not given in this form.Currently there are two main approaches to address this problem within the MCTDH framework: the fitting of the PES to product form (e.g., using the potfit algorithm) 45,46 or the Correlation DVR (CDVR) approach. 47The latter approach, which uses a time-dependent quadrature corresponding to the SPF basis and allows the use of general potentials without fitting them to product form, has been used in this work for the calculation of the eigenvalues of the reduced 5D Hamiltonian.However, once the diabatic basis is built, CDVR is no longer needed, since the diabatic potential energy matrix depends only on a single coordinate.
Besides the regular real-time propagations, the MCTDH approach can also be used to obtain eigenvalues and eigenfunctions of certain relevant operators which imply imaginary-time propagations.In the present work, block relaxation of stateaveraged MCTDH (SA-MCTDH) wave functions 40 is used to obtain the eigenstates and eigenvalues of the 5D Hamiltonian operator, Ĥ5D (q; z).Furthermore, the eigenvalues and eigenstates of thermal flux operator are calculated by the iterative diagonalization of FT within the state-averaged MCTDH approach. 40

III. RESULTS AND DISCUSSION
In Sec.II B it was argued that the diabatic approach is formally exact in the limit of a complete basis set.In the case of a quasi-separable system, with a very small coupling between confined and unbound coordinates, convergence will be achieved with few eigenfunctions.Furthermore, an adiabatic approximation which separates the motion of the fast confined degrees of freedom from that of the slower unbound ones will yield good results.In Sec.III, we will first discuss the validity of the model by studying the cePESs obtained with the adiabatic approach and the 5D Hamiltonian matrix in the diabatic representation.Then the approach will be used to rigorously compute CRPs using the diabatic representation.
The accuracy of the adiabatic approximation will be studied by comparison with full-dimensional Hamiltonian results.Finally, the diffusion rates of H 2 into the hollow cavity of the nanotube will be computed.

A. Diabatic Hamiltonian matrices and adiabatic cePESs
In the present work, a total of 50 confined 5D eigenstates were used as a basis for the diabatic (and adiabatic) representation of q.To ensure convergence of the 50 ξ j (z 0 ) basis functions, a higher number of states was calculated in the SA-MCTDH diagonalization, namely, 80.The numerical parameters for the MCTDH wave function representation used to calculate the confined eigenstates are given in Table I: SPF basis size and primitive's grid type, number of points, and the span of the representation.Note that the primitive's grid type is a FFT grid, except for the angular coordinate θ, for which a Cotangent DVR (cot-DVR) 48 was chosen.This representation was constructed specifically to avoid the singularity appearing due to the term 1 sin θ in the Hamiltonian, Eq. ( 1).See Ref.48 for more details.The calculation of the 5D eigenstates at a value of z = −1.36 a 0 , corresponding to a minimum in the electronic PES, took 14 h in a 16 core Intel Xeon E5-4620 0@2.20 GHz machine.The subsequent calculation of the diabatic Hamiltonian matrix and the adiabatic cePESs for 512 equidistant values in the z ∈ [−56.066,56.066] a 0 range-corresponding to the length of 14 (8,0)-SWCNT unit cells-took only few seconds in the same machine.
Through the analysis of the adiabatic surfaces, one can detect regions where avoided crossings appear, and therefore the coupling between states could be stronger.The adiabatic surfaces ε j (z) of the H 2 molecule inside the CNT are plotted in Fig. 2 as a function of the diffusion coordinate, z.It should be noted that only states transforming according to the same irreducible representation of the system's symmetry group can interact.cePESs corresponding to states of different symmetry can thus cross without giving rising to any coupling.In Fig. 2, cePESs corresponding to states which are symmetric with respect to the permutation of the two hydrogen atoms are displayed by full lines while cePESs of antisymmetric states are shown as dashed lines.
In Fig. 2, it is readily seen that the lower cePESs are wellbehaved, smooth functions and therefore the coupling between them is expected to be negligible.However, as one goes up in energy, the density of states is largely increased and avoided crossings appear.This is seen in the higher energy surfaces in Fig. 2. Considering the z-dependence, one finds that near the minima of the cePESs all surfaces are well distinguishable.
No avoided crossings are found in this region and the coupling between the different cePESs can be expected to be negligible.However, in the steep regions of the PES and near the maxima, avoided crossings are visible.The coupling between different confined eigenstates is potentially significant in these areas.
Following the analysis of the cePESs, we can graphically represent the 5D Hamiltonian matrix of Eq. ( 8) at different points of z, paying special attention to the areas where avoided crossings appear in the cePESs.This representation will provide further information about which confined states, if any, are significantly coupled: large off-diagonal matrix elements would indicate significant couplings and a breakdown of the adiabatic approximation.On the other hand, if the couplings are small and H 5D is quasi-diagonal, the adiabatic approximation can be expected to yield good results.
Furthermore, in this case just a few basis functions will be required in the diabatic representation in order to obtain accurate results.
Figure 3 shows a representation of the H 5D matrix for H 2 , using fifty 5D eigenstates as a basis for q, at two points along the PES: a minimum of the PES (z = 1.36 bohr), and at z = 0 bohr (indicated by the vertical dotted line in Fig. 2), a point where many avoided crossings of the cePESs appear and the largest couplings can be expected.The matrix elements are represented in logarithmic scale for the sake of clarity.As it was expected, the off-diagonal terms at z = 1.36 bohr, the reference point used for the definition of the diabatic basis, effectively vanish.More interesting are the results obtained at a point in the strong coupling region at z = 0 bohr.Even here, the largest off-diagonal terms are roughly 50-100 times smaller than the diagonal elements.This finding indicates a very low coupling between the unbound motion in z and the dynamics in the confined DOFs q.It strongly suggests the suitability of an adiabatic separation of the time scales in studies of the quantum dynamics of nanoconfined species.

B. Dynamics of H 2 diffusion
The two approaches described in Secs.II B and II C were used to run quantum dynamics simulations of H 2 diffusion inside the SWCNT.Cumulative reaction probabilities were obtained at a reference temperature of 100 K and the diffusion rate was computed by thermal averaging of N(E).Additional calculations were preformed at a reference temperature of 1000 K to confirm the accuracy of the results.
In order to obtain converged results in the diabatic and adiabatic schemes, several parameters need to be optimized, namely: the primitive basis set used in both q and z, the propagation time, the complex absorbing potential (CAP), and the number of flux eigenstates contributing to N(E) in the chosen temperature range.Using the relation between flux eigenstates and eigenstates of the confined Hamiltonian, we can get an approximate threshold for the number of relevant flux eigenstates at a given temperature: they will correspond to approximately twice the number of confined eigenstates relevantly populated at that temperature.In our case, for the highest temperature considered, 1000 K, up to 11 confined eigenstates would be significantly populated.For the sake of numerical accuracy, 26 flux eigenstates were obtained by the iterative diagonalization of the thermal flux operator.These eigenstates were propagated in real time to finally obtain N(E).Regarding the primitive basis, converged results were obtained using 50 diabatic basis functions in the q coordinate and 20 SPFs in the unbound DOF z.The different basis sets used in the diabatic and adiabatic calculations are given in Table II.
Finally, both time convergence and CAP optimization required a special treatment due to the very low energies involved in the diffusion process.Regarding the CAP, it was seen that regular polynomial forms of the imaginary potential failed to avoid both reflexions and transmissions of the wave packets in the energy range involved.To overcome this issue, a transmission-free absorbing potential 49,50 was used.This potential has the advantage that no transmission of the wave packet is possible.Hence, one could in principle reduce the amount of reflexion arbitrarily by increasing the length of the imaginary potential.In our case, to avoid all possible reflections, a length of 20 Å was required.The thermal flux eigenstates were propagated for 10 ps which allows one to obtain the symmetric flux correlation in Eqs. ( 17) and (19) for times up to 20 ps. 41s seen in Eq. ( 17), the cumulative reaction probability is obtained as the Fourier transform of a flux autocorrelation function.Although an infinite propagation time is in principle needed to obtain N(E), generally convergence is achieved at relatively short times, as the wave packet leaves the interaction region.A measure to qualitatively check the convergence of a quantum dynamics calculation is given by the flux-flux correlation function, an its time integral, the flux-position correlation function, Generally, this function increases and potentially oscillates while the wave packet remains in the interaction region, but asymptotically reaches a constant value.The total C fp for both approaches developed in the present work, together with the full dimensional data obtained from a previous work, 33 is shown in Fig. 4 for the reference temperature of 100 K.There are several conclusions to be drawn from the exploration of this quantity.First, focusing in the short-time region (see the inset of Fig. 4), it is seen that the results obtained using the three different approaches agree almost perfectly.The agreement between the data from our previous 6D work and the new results obtained using the diabatic approach described in Sec.II B confirms that both numerically exact calculations (6D and diabatic representation) are numerically converged.Furthermore, it is found that the adiabatic approximation yields very accurate results at this time scale.This finding confirms the separation of time scales of the dynamics in the unbound and constrained DOFs.Second, one notes the complex, long-living structures in the flux-position correlation function which does not reach to a constant value even at 20 ps.This behavior indicates the existence of long-living resonances.Due to the long propagation times required, these structures were not studied in a previous work.Their investigation in the present study is facilitated by the numerical efficiency of the diabatic approach described in Sec.II B.
Finally, it must be noted that C fp resulting from the diabatic and adiabatic propagations are very similar until about 10 ps but differ significantly beyond this value.The differences indicate that the coupling between the motion in unbound and constrained DOFs becomes relevant at about 10 ps.This time scale corresponds to coupling constants of about 0.06 meV.It should be noted that the gradual decay of accuracy resulting from the assumption of separability of time scales in the adiabatic calculation is consistent with the results reported in a previous work by some of the authors. 32Here it was shown that the coupling between unbound and confined degrees of freedom increased each time a wave packet crossed a maximum of the PES.
To investigate the resonances giving rise to the long-living structures in the flux-correlation function in more detail, the limited resolution CRP was computed as defined in Eq. (18), N(E), with a resolution of ∆E = 0.12 meV, and is displayed in Fig. 5. Comparing the rigorous results obtained by the diabatic approach with the results of the adiabatic approximation, one immediately finds that the resonance structures are well described within the adiabatic approach.The agreement between the results of the diabatic and adiabatic calculations is perfect at low energies, almost perfect at energies up to about 0.4 eV, and reasonably good at higher energies.It should be noted that the results of diabatic calculation show a slightly more noisy behavior at higher energies indicating remaining numerical inaccuracies due to imperfect convergence.FIG. 5. Comparison of the total N(E) for the diffusion of H 2 in an (8,0) CNT obtained through the diabatic (solid, black) and adiabatic (dashed, red) approaches.
As it was hinted by the complex shape of the flux-position correlation function, N(E) presents a number of structures, specially in the lowest energy region.These can be nicely explained studying Fig. 6, where the CRP obtained in adiabatic approximation is plotted together with the individual contributions arising from the first five flux eigenstates.In order to investigate the resonant character of the peaks found in the CRP, we have calculated the eigenstate spectrum of H 2 confined to a single unit cell of the carbon nanotube, i.e., applying the periodic boundary condition at the limits of the box displayed in Fig. 2. The computed eigenenergies are plotted in Fig. 6 as vertical lines.These eigenstates appear in pairs corresponding to which are either symmetric or antisymmetric with respect to reflection in z-planes located at the top of the potential barriers.In Fig. 6, solid vertical lines correspond to symmetric eigenstates while dotted ones are used for anti-symmetric ones.
The connection between the eigenenergies and the resonances is clearly visible in the figure: the contribution of a given flux eigenstate to the CRP, N i , rises to unity whenever the total energy coincides with the value of a symmetric eigenstate.On the other hand, it decays again to zero when the energy FIG. 6.Total N(E) and the first five individual contributions to the CRP obtained with the adiabatic approach at a reference temperature T 0 of 100 K. N 3 and N 4 come from degenerate states.Vertical lines correspond to the eigenstates of the 6D Hamiltonian computed in a unit cell with periodic boundary conditions (see text for details).FIG. 7. Diffusion rate for H 2 in an (8,0) CNT computed using the diabatic approach (solid line) and the adiabatic approximation (dashed).TST results (dotted) are also given for comparison.
corresponding to an anti-symmetric eigenstate is reached.In the energy interval between the two eigenenergies, resonance enhanced tunneling through the barrier is possible.Comparing the resonance positions in Fig. 6 with the cePESs in Fig. 2, the low-lying resonances in the CRP can be straightforwardly assigned.The two prominent resonances below 0.355 eV result from resonance enhanced tunneling on the two lowest cePESs.The corresponding resonance wave function is mainly localized in the potential wells and shows no nodes in this region.The corresponding tunneling splitting is in the meV range.The next eigenstates appear about 5 meV above the first ones.Comparing this excitation energy with the barrier height of about 10 meV on the lowest cePESs, one would expect a much larger tunneling splitting.Actually the corresponding tunneling splittings seen in Fig. 6 are about 4 meV, which confirms this expectation.The two overlapping CRP contributions resulting from these states give rise to the interesting feature seen at about 0.36 eV.At even higher energies, the available energy exceeds the barrier height and both channels corresponding to the two lowest cePESs are completely open.Therefore the CRP settles at a value of two until the first resonance states located on the higher cePES appear between 0.38 and 0.39 eV.
After Boltzmann averaging and integration of the CRP at different temperatures, we first obtain k(T ) and then, using Eq. ( 14), the diffusion rate.This quantity, computed by both approaches developed in this work, is shown in Fig. 7 for the H 2 molecule in an (8,0) CNT.Transition state theory results are also plotted for comparison.As it was expected from the study of the different N(E), the results obtained with the diabatic and adiabatic formalisms agree perfectly, both showing a significant tunneling contribution below 125 K.

IV. SUMMARY AND CONCLUSIONS
The differences in the time scales of the fast motion in the confined DOFs and the slow motion along the nanotube's axis in the H 2 @(8,0) SWCNT system were exploited to devise efficient numerical methods to study the system's quantum dynamics.In the rigorous, numerically exact approach, first a diabatic basis and a corresponding Hamiltonian matrix describing the motion in the confined coordinates is constructed.The resulting diabatic Hamiltonian matrix only depends on the z-coordinate describing the unbound motion of H 2 along the nanotube.Thus, the 6D Schrödinger equation is turned into a 1D equation describing the motion on a set of coupled diabatic potential energy surfaces.Due to the small size of the couplings, converged results can be obtained with a small number of diabatic basis states.The resulting scheme is extremely numerically efficient and allows one to study the system's dynamics on a time scale which could not be accessed in previous studies. 32,33nvoking an adiabatic approximation to decouple the slow and fast motions, the description of the system's dynamics can be further simplified.In the adiabatic approximation, the dynamics precede on uncoupled adiabatic potential energy curves.Comparison with rigorous results shows that the adiabatic approximation works very well for the present system.The adiabatic picture is then successfully utilized to interpret the system's dynamics.
The investigation of the long-time dynamics in the H 2 @(8,0) SWCNT system facilitated by this methodological development yielded interesting new insights.It was found that the motion of H 2 along the nanotube axis at low energies can occur via resonance enhanced tunneling.Corresponding structures are clearly visible in the computed energy-dependent transmission coefficients.The resonance can be assigned to quantum states localized in the wells of the adiabatic PESs.This result is particularly interesting as it suggests that tunneling might be the dominant mechanism for the motion of individual H 2 molecules trough narrow CNTs a low temperature.The consequences of this finding should be investigated more deeply, since it might affect potential applications of carbon nanotubes, such as their usefulness as quantum sieves.

FIG. 1 .
FIG. 1. Scheme of the coordinate system used in the present work.

FIG. 3 .
FIG. 3. Graphic representation of the Hamiltonian matrix in the basis of SPFs at the point of minimum (right) and maximum (left) coupling.Energy scale is in eV.

TABLE I .
Basis set and grid MCTDH representation used in the calculation of the 5D eigenstate basis for H 2 .
3.5-3.5 a 0 FIG.2.Adiabatic cePESs generated by the 30 lowest energy eigenstates of the H 2 @SWCNT system.cePESs corresponding to confined eigenstates which are symmetric with respect to permutation of the two hydrogen atoms are displayed by black solid lines while cePESs corresponding to antisymmetric states are shown as red dashed lines.

TABLE II .
Basis set and grid MCTDH representation used in the propagations.