A nuclear spin and spatial symmetry-adapted full quantum method for light particles inside carbon nanotubes: clusters of 3 He, 4 He, and para -H 2 †

We present a new nuclear spin and spatial symmetry-adapted full quantum method for light fermionic and bosonic particles under cylindrical carbon nanotube confinement. The goal is to address Fermi– Dirac and Bose–Einstein nuclear spin statistics on an equal footing and to deliver excited states with a similar accuracy to that of the ground state, implementing ab initio -derived potential models as well. The method is applied to clusters of up to four (three) 4 He atoms and para -H 2 molecules ( 3 He atoms) inside a single-walled (1 nm diameter) carbon nanotube. Due to spin symmetry effects, the bound states energy landscape as a function of the angular momentum around the tube axis becomes much more complex and rich as the number of 3 He atoms increase compared to the spinless 4 He and para -H 2 counterparts. Four bosonic 4 He and para -H 2 particles form pyramidal-like structures which are more compact as the particle mass and the strength of the inter-particle interaction increases. They feature stabilization of the collective rotational motion as bosonic quantum rings bearing persistent rotational motion and superfluid flow. Our results are brought together with two key experimental findings from the group of Jan-Peter Toennies: (1) the congestion of spectral profiles in doped 3 He droplets as opposed to the case of 4 He droplets (S. Gebenev, J. P. Toennies and A. F. Vilesov, Science , 1998, 279 , 2083); (2) the onset of microscopic superfluidity in small doped clusters of para -H 2 molecules (S. Grebenev, B. G. Sartakov, J. P. Toennies and A. F. Vilesov, Science , 2000, 289 , 1532), but at the reduced dimensionality offered by the confinement inside carbon nanotubes.


Introduction
Being composed by the lightest atoms and molecules in nature, aggregates of helium and molecular hydrogen are the most paradigmatic cases of finite quantum systems, with the atoms and molecules being held together by weak dispersion forces. The solvation of molecular dopant with helium atoms and hydrogen molecules has allowed to find evidences for the onset of superfluid behaviour at microscopic level and to show the role of the Fermi-Dirac or Bose-Einstein nuclear spin statistics. The molecular spectroscopic measurements in helium droplets by the group of Professor Jan-Peter Toennies demonstrated the existence of very different spectral profiles depending on the considered isotope, 3 He or 4 He: 1,2 As first observed for sulfur hexafluoride, SF 6 3,4 and carbonyl sulfide, OCS, 1 the spectra of molecular impurities inside 4 He droplets attain a very well resolved profile with sharp rotational lines. In contrast, the spectra of OCS inside 3 He droplets exhibited broad (about 1 cm À1 ) spectral profiles, with the temperature of the drop (about 0.1 K) being below the superfluid transition temperature in bulk 3 He (T l = 3 Â 10 À3 K). Moreover, a well-resolved spectral profile was recovered after adding 4 He atoms to 3 He droplets. 1 Since 4 He atoms can be considered as composite spin-less bosonic particles while 3 He atoms are fermionic with a nuclear spin equal to 1/2, the essential role of nuclear spin statistics on microscopic superfluidity was demonstrated. Further experimental observations by the same group 5,6 indicated the onset of molecular superfluidity even in small doped clusters of spin-less bosonic para-H 2 (or ortho-D 2 ) molecules. Further theoretical 7,8 and combined experimentaltheoretical evidences 9 on the superfluidity of small para-H 2 clusters have been found. The cylindrical confinement provided by carbon nanotubes has offered the possibility of studying the pronounced quantum behaviour of 4 He atoms and H 2 molecules at reduced dimensionality. This way, very recent measurements have demonstrated the formation of two-dimensional (2D) 4 He layers on the outer surface of a single-walled carbon nanotubes (SWCNTs). 10 Also, an experimental study 11 of gas adsorption at low (2-5 K) temperature revealed a quenched propagation of 4 He atoms through carbon nanopores with diameters below 7 Å despite of their small kinetic diameter. The application of orbital-free helium density functional theory to carbon nanotubes immersed in a helium nanodroplet provided theoretical explication that the experimental observations stem from the exceptionally high zero-point energy of 4 He as well as its tendency to form 2D layers upon adsorption at low temperatures. These conclusions were further confirmed by applying more accurate ab initio potential modelling along with a wave-function (WF)-based approach. 12 Moreover, the quest for the occurrence or either quasi-superfluid or quasi-crystalline behavior of para-H 2 molecules when confined in SWCNTs has recently motivated an intense theoretical research. [13][14][15] For the case of molecular deuterium, our recent theoretical work has provided conclusive evidence for the transition from molecular aggregation to quantum solid-like packing in a SWCNT of 1 nm diameter, 16 confirming a previous study using an embedding approach in a broader (ca. 1.4 nm) SWCNT. 17 Altogether, these studies [12][13][14][15][16][17][18] have highlighted the key role played by the quantum nature of the nuclear degrees of freedom in the confined helium, hydrogen, or deuterium motion.
The present work is aimed to provide a new symmetry-and spin-adapted full quantum approach capable to describe fermionic particles of spin 1/2 ( 3 He atoms) as well spin-less bosonic particles ( 4 He atoms and para-H 2 molecules) inside SWCNTs. A similar case of 3 He, 4 He and para-H 2 particles around a molecular dopant was successfully treated by the Full-Configuration-Interaction-type Nuclear Orbital (FCI-NO) approach (see, e.g., ref. 8). This FCI-NO treatment was able to accurately deal with the hard-core interaction problem, 19 and to automatically include all the nuclear-statisticinduced (bosonic or fermionic) effects. It handled the many-body problem by expanding the Hamiltonian operator with quantum single-particle basis functions followed by an efficient diagonalization technique. 20 The method was rooted on a previous Hartree-Fock-type implementation 21 which was applied to explain the role of Boson-Fermion statistics on the differences in the spectra obtained by Toennies and collaborators 1 depending on the considered isotope ( 3 He or 4 He). 22 The FCI-NO method resembles the traditional CI method of quantum electronic structure theory and was originally developed for doped fermionic 3 He clusters 19,23,24 and further extended for bosonic systems (for a mini-review see ref. 25) composed by 4 He atoms 26 or para-H 2 molecules. 8 As the method presented in this work, one key advantage of the FCI-NO technique was its capability of providing solvent states with accuracy similar to that for the ground state. This capability was used to characterize collective rotational motion, 8 and to provide insights into the experimental measurements indicating the onset of microscopic superfluid behaviour on doped para-H 2 clusters by Toennies and collaborators. 5,6 The FCI-NO method can not be directly applied to the case of SWCNT as the 'dopant' species. In fact, a heavy dopant in helium (or molecular hydrogen) clusters can be considered as providing a reference coordinate system fixed in space, for the helium (or molecular hydrogen) motion. 8,19,[23][24][25] This ansatz is similar to the adiabatic approach in electronic structure calculations, where the nuclear framework is considered fixed and provides a coordinate system for electrons. Therefore, the electronic WF can be expanded as a sum of products of 'molecular orbitals', with each of them being well localized in space. Similarly, the FCI-NO technique uses an expansion of products of 'nuclear orbitals', being centered at the molecular dopant and fixed in space. In the present case, the SWCNT 'dopant' is also fixed in space but is infinitely long and uniform in one direction. By assuming the coordinate axis Z along the tube axis, the overall motion of the cluster along Z is thus free. Consequently, the spectrum is continuous, with the total WF being delocalized along Z. In order to obtain the bound states, the motion of the cluster center-of-mass (COM) along the coordinate Z needs to be explicitly separated. We know that the total WF, corresponding to a zero momentum along Z and the target bound state, must be independent on Z. Unfortunately, this condition is impossible to satisfy when the WF is taken as a product of one particle 'orbitals'. Thus, unless all orbitals do not depend on Z, the artificial COM kinetic energy will contaminate the total energy. This difficulty is reminiscent to the problem of separating the overall COM in theories where the quantum motion of electrons and nuclei is tentatively described on an equal footing using orbital-like approaches (see, e.g., ref. 27 and references therein). In the present work, this issue is solved by recasting the Hamiltonian in internal coordinates and incorporating an adaptation to the spin-and spatial-symmetry of the confined fermionic and bosonic particles.
A consistently accurate description of the interaction between the confined species and the SWCNT is achieved through the application of ab initio-derived pairwise additive potential models for both helium 12 and hydrogen. 17 For the specific case of van-der-Waals-dominated adsorbate/(carbonbased) surface interactions, 12,17,[28][29][30] its adequacy relies on the excellent transferability properties found for the parameters accounting for the dispersion contribution upon increasing the size of clusters modelling the extended systems. 28 More difficult, dispersionless contributions are mostly of short-range nature and the usage of small cluster models can be also justified. Therefore, different analytical forms are employed to model the dispersionless and dispersion contributions, with the model parameters being computed at ab initio level on short and narrow carbon nanotubes and then scaled to the actual systems. The detailed energy decomposition offered by the ab initio method applied -density functional theory (DFT)based symmetry-adapted perturbation theory (SAPT(DFT)) 31,32has rendered such strategy successful and efficient when describing the adsorption of molecular and atomic species to inner and outer surfaces of carbon nanotubes. 12,17,33 The article is structured as follows. Section 2 presents the method developed in this work. Section 3 illustrates its application to the cluster of quantum particles ( 3,4 He atoms and para-H 2 molecules) inside a SWCNT of 1 nm diameter. We emphasize on the interplay between the mass of the given particles, the strength of the corresponding pair interactions, and their nuclear (fermionic or bosonic) spin statistics, as well as on the onset for the stabilization of collective rotational motion. Finally, Section 4 summarizes our main findings and presents an outline of future perspectives.

Setting up the problem
The method presented in this work has been developed to allow the calculation of the bound states of N identical atoms or molecules trapped inside a SWCNT. These atoms or molecules will be referred with the term 'pseudo-nuclei' (PNs), having been described as pointlike particles of mass M. The total interaction of the PN cluster with the SWCNT is represented as a sum of one-particle terms, describing the interaction of the SWCNT with a single PN, and the interaction between the PNs. The latter is written as a sum of pairwise terms. We assume that the corrugation is negligible so that the interaction of each PN with the SWCNT (which is assumed to be infinitely long) depends only on the distance of the PN from the nanotube axis, V 1 = V 1 (r). The interaction between each pair of (identical) PNs is assumed to be a function of just the PN-PN distance, V 2 = V 2 (r ij ). Then, the total Hamiltonian operator can be written as a function of the 3 N Cartesian coordinates of the PNs aŝ Fig. 1), together with a Jacobian transformation of the volume to keep the Hamiltonian 'explicitly Hermitian', C ! ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 1 r 2 . . . r N p C; it takes the form where the distance r ij is explicitly given by Considering a SWCNT of 1 nm diameter and helicity index (11,4), the He-CNT and H 2 -CNT one-particle potentials V 1 (r) are presented in the upper panels of Fig. 2 [panels (a) and (b)] along with a few supported bound states. It can be seen that the carbon wall is impenetrable, with the PN-SWCNT being dispersion-dominated. As can be observed in the lower panels of Fig. 2 [panels (c) and (d)], the two-particle interaction V 2 (r ij ) between the PNs is much weaker than the PN-SWCNT interaction, specially in the case of helium. Notice the huge zeropoint energies for the 4 He-SWCNT (ca. 500 cm À1 ) and 4 He-4 He vibrational motions.
Since the center-of-mass coordinate leads to a continuum spectrum, the Hamiltonian has to be recast in internal coordinates which do not contain the overall Z-translation neither the overall rotation coordinate around the Z axis as the PN-SWCNT and PN-PN interactions do not depend on it.
We introduce the following internal coordinates for the z degrees of freedom The angular relative coordinates are defined instead as It should be mentioned that the definition of the overall rotation coordinate is not unique. We use F = f N instead of   with the overall rotation quantum number L given by After separating the coordinates Z and F, a L-dependent Hamiltonian of 3N À 2 internal coordinates is obtained aŝ The kinetic energy operator along the Z axis can be expressed aŝ while the rotational kinetic energy term iŝ : For L 4 0, the angular kinetic energy contribution contains complex terms. Yet, the whole Hamiltonian matrix can be kept real if the total WF is separated into two terms, symmetric and anti-symmetric with respect to a simultaneous change of sign of all w i coordinates, correspondingly, and the anti-symmetric term is taken with the imaginary factor 'i'.

Adaptation to the symmetry of the problem
The total WF of a system composed by N identical particles must be symmetric (for bosonic particles, e.g., 4 He atoms) or anti-symmetric (e.g., 3 He atoms) with respect to the exchange of any two particles. Moreover, the Hamiltonian expressed in both Cartesian and internal coordinates must be invariant under z i -Àz i and, for L = 0, f i -Àf i transformations.
While z and f symmetries are equally present in the internal Hamiltonian, and thus can be easily accounted for, the incorporation of the full permutation symmetry of all N particles (i.e., the S N group) is more complicated. The irreducible representations of this group are described by Young diagrams. The fully symmetric (anti-symmetric) solution for bosons (fermions) correspond to the {N} ({1,1,. . .,1}) Young diagrams. The case of fermions is even more complex as the spin symmetry has to be incorporated as well. Using the FCI-NO approach, 8 all these symmetry properties can be easily implemented via the second quantization technique. The employment of internal coordinates makes the symmetry adaptation far more complex. The internal Hamiltonian is explicitly invariant under any exchange between the first N À 1 particles, and not invariant when the particle labelled as N is involved. It means that we can use only the S NÀ1 symmetry group representations explicitly and it is necessary to further perform a S N -S NÀ1 reduction. This aspect is discussed in detail below.
The case of N = 2 is special and the full permutation symmetry can be explicitly taken into account. Thus, if the angular coordinates are settled as ; the exchange of the two particles 1 2 2 simply corresponds to Notice that the factor 1/2 in the definition of w makes the WF periodic in w, with the associated quantum number being m = m 1 À m 2 . The overall rotational quantum number L = m 1 + m 2 is always integer. Yet, an additional constraint must be imposed on the w dependence as both m 1 and m 2 are integers and thus the parity of m should be the same as the parity of L. In other words, the wave-function must be AEinvariant under the ww + p transformation with a phase +1 for even L and À1 for odd L. This condition is satisfied when using properly adapted basis functions, as explained in Section 2.4. The angular kinetic energy in these coordinates is expressed aŝ which is indeed explicitly invariant under the 1 2 2 permutation.
To conclude this subsection, let us describe how a given explicit symmetry is imposed into the Hamiltonian. As iterative Davidson-type techniques are used for the Hamiltonian diagonalization (see Section 2.4), the employment of symmetryadapted basis functions is not efficient. Instead, we use a projection technique and, specifically, character-based projectors. The projector to a given irreducible representation a of the symmetry group is expressed as, 34 where the sum runs over all group operators, f a stands for the dimension of representation, |G| denotes the group order, and w (a) are the characters of the representation. This projector is applied all along the iterative procedure to ensure that the calculated solutions belong to the required representation.
2.3 Reduction S N -S NÀ1 and incorporation of the spin symmetry As mentioned above, the Hamiltonian expressed in internal coordinates does not explicitly include the symmetry group S N . Yet, the solution must belong to the S N representation as the original Hamiltonian features this symmetry. In order to identify the specific S N representation, the symmetries of the solutions are tested a posteriori. An acceptable solution for the case of spinless bosons must belong to the {N} Young diagram. By removing the coordinate of the particle labelled as N, this solution might only correspond to the {N À 1} diagram of the S NÀ1 group as it must obviously be symmetric with respect to any permutation of first N À 1 PNs. However, the solution that transforms as the {N À 1} representation of the S NÀ1 symmetry group may correspond not only to the {N} but also to the {N À 1,1} diagram of S N . To verify that the good {N} solution has been obtained, we calculate a 'symmetry factor Q' as: Then, from general theory of the S N group representations, 35 it can be shown that the factor Q for the solution corresponding : Therefore, the factor Q can be used to select the correct solutions. The details of its calculation when applying the DVR technique are given in Section 2.4. It involves a modification in the grid, and thus this factor may also serve as an indicator of the accuracy of the numerical treatment (grid quality) itself. A similar approach can be used for the fermionic case. Here we may have degenerate solutions belonging to different representations of the S NÀ1 subgroup so that several Q factors need to be calculated. We provide full details for bosonic and fermionic cases with explicit expressions for N = 2, 3, and 4 in Section S1 of the ESI. † It is also stressed that the differences between the energies of correct (pure fermionic or bosonic) solutions and those corresponding to solutions which do not satisfy the associated full S N symmetry serve as an estimate of the magnitude of exchange effects. In the calculations presented in this work, these energy differences have been of the order of a few wave-numbers.

Implementation: discrete variable representation (DVR) approach
To evaluate the matrix elements of the internal Hamiltonian, we use the Discrete Variable Representation (DVR) approach 36 as in our previous works, 16,17,33,37 with the basis set obtained as a direct product of functions for the different coordinates. The main advantage is that the potential term V 2 is diagonal in the DVR basis so that it needs to be evaluated just on the set of DVR grid points.
To fully account for the stronger V 1 interaction, we use the potential-optimized (using V 1 ) PO-DVR method 38 (see also ref. 16). The underlying DVR functions are 2D harmonic oscillator functions as they 'absorb' the 1 4r 2 singular term at r = 0 into the kinetic energy contribution. As the N = 1 excitation energies are much larger than the clusters formation energies (because of the V 2 contributions are much smaller than the V 1 contributions), a few PO-DVR functions are sufficient. The good performance of such a basis has been assessed by the fast convergence of the total energies with the number of PO-DVR grid points along the r coordinates, n r .
Sinc-DVR functions are employed for both t-and w-coordinates. The number of functions on the t-coordinate and the grid interval is chosen to ensure the convergence of the bound-state WFs and energies. Sinc-DVR functions are defined for the angular coordinates w on the interval [0. . .2p], with the kinetic energy adjusted to the periodicity of the WFs. These functions are symmetrized with respect to the w -2p À w transformation and the anti-symmetric functions are taken with phase factor i so that the Hamiltonian matrix is kept as real. For the particular N = 2 case, the Sinc-DVR functions are also symmetrized with respect to the ww + p transformation so that even-m and odd-m functions can be distinguished.
The diagonalization of the resulting Hamiltonian matrix is performed using the Jacobi-Davidson (JD) algorithm. 39 In contrast with standard Davidson algorithms, the JD technique has been found to be very robust and capable of providing converged energies even for ill-behaved 'hard-core' interactions such as the short-range He-He interaction in doped helium clusters. 19,20,24 As explained above, the factors Q are evaluated to assign correctly the permutation symmetry. Their explicit expression is: where . .N À 1. Equivalent expressions hold for the angular coordinates w. To derive these relations, we must consider that the overall rotation is accounted for via the coordinate F = f N , with exp(iLf N ) as the corresponding contribution to the total WF. For L 4 0, the factors exp(ÀiLw 1 ) are calculated via the DVR approach. The Q factor evaluation obviously requires the interpolation of the DVR functions to the transformed grid points. Its deviation from the 'exact' Q values is used as a numerical test for the quality of the DVR basis set expansion.

Interaction potentials
As mentioned in the introduction, we have used previously developed ab initio-derived potential models for the PN-SWCNT interaction (PN = He, H 2 ), using different analytical forms for the dispersionless and dispersion energy contributions. To preserve the cylindrical symmetry of the system, the very small corrugation appearing along the azimuthal degree of freedom f is ignored. Moreover (see e.g., ref. 12), the adsorbate-SWCNT interaction is almost constant along the Z axis, decaying slowly at the borders. Then, considering long nanotubes, we can assume that the interaction potential depends only on r. The atomic structure of the SWCNT is implicitly considered through the SAPT(DFT) calculations used to fit the parameters of the employed pairwise potential model. 12,17 For convenience, the details are given in Section S2 of the ESI. † The adequacy of the potential model has been attested in previous works on both the He-SWCNT (see Fig. 2 of ref. 12) and the H 2 -SWCNT (see Fig. 3 of ref. 17) interactions. In order to employ these potentials within the DVR approach, they were fit to the polynomials of r: The parameters of these polynomial fits are provided in Section S2 of the ESI. † To model the He-He interaction, we have employed the very accurate potential of Cencek et al., 40 which includes nonadiabatic and relativistic corrections. For the H 2 -H 2 interaction, the rotationally-averaged potential energy surface reported by Hinde 41 has been employed instead, as in our previous works. 16,17 3.2 Bound states of ( 3,4 He) N and (para-H 2 ) N clusters (N o 4) Table 1 summarizes the results for the bound states of ( 3,4 He) N and (para-H 2 ) N clusters inside the SWCNT(11,4) tube for N r 3. To better understand the spatial structure of the WFs, Table 1 Bound states of ( 3,4 He) N and (para-H 2 ) N clusters inside a SWCNT (11,4) nanotube, with N r 3. For L 4 0, the energy is given with respect to that corresponding to the L = 0 state. The binding energy E bind is defined as the negative of the difference between the energy of a given bound-state and that corresponding to the dissociation channel (which is also indicated). Under dissociation, the total L value (i.e., the sum of the L's values of the monomers) is preserved. 42 For N = 1, the L value is also denoted as m. Using our basis sets (see appendix), the energy values are converged to better than 0.01 cm À1 for N r 3. The labels '+' and 'À' indicate that the wave function is symmetric (+) or anti-symmetric (À) with respect to the transformation z -Àz (referred to as z-symm)  Considering the lowest-energy bound states, we can observe that the values of the binding energies of (ortho-H 2 ) N complexes are at least 10 times larger than those associated to ( 4 He) N clusters. Comparing the two helium isotopes, we notice that the   binding energy of ( 4 He) N complexes is twice as large as for ( 3 He) N clusters. It is remarkable that, in contrast with the free dimer, the confined ( 3 He) 2 system is bound (although by just À0.3 cm À1 ), despite of the free Z motion. In fact, according to a previously developed one-dimensional model (see Section S3 of the ESI of ref. 12), the effective 3 He-3 He interaction along the tube axis would be very similar to that of the free dimer in such a narrow nanotube of 1 nm. This finding can be explained as follows: the confinement causes the single 3 He distribution to be much more localized along the r direction. By accounting for the z symmetry and according to the obtained density distributions, the two 3 He atoms can place themselves at opposite sides of the nanotube diameter and not along the tube axis. Then, the attractive V 2 interaction makes possible that such configuration is bound, even at the price of a slightly increased kinetic energy (f localization).
The differences between the bound states of confined 3 He and 4 He atoms for N = 1 originates from the higher zero-point energy (the mass) of the former. Due to the spin symmetry, however, the levels structure for N 4 1 and the 3 He isotope (as well as the dissociation channels for N = 3) differs much from the 4 He counterpart. While there is only one bound state for L = 0, the number of states of opposite t-symmetry for L 4 0 increases sharply as the number of 3 He atoms grow. In particular, notice that a significant number of bound states have very close energies already for the ( 3 He) 3 cluster. Of course, this trend is expected to be even more pronounced for larger N values. This outcome is reminiscent to that reported for doped 3 He clusters at both Hartree-Fock 22 and FCI-NO levels. 25 In both studies, 22,25 the congested profile of the experimental spectra in doped 3 He droplet 1 was attributed to the high energy degeneracy exhibited by fermionic clusters.
We can also observe from Table 1 that collective rotational states for ( 4 He) N clusters with L = N are associated to larger binding energies, as opposed to those of ( 3 He) N clusters. However, no energy minima are obtained at L = N for either 4 He or para-H 2 bosonic particles. This feature would signal the stabilization of collective rotational motion which, importantly, is found upon the aggregation of four bosonic particles (see Section 3.3).
It worth stressing that the clusters dissociation often breaks the symmetry in t, which is present in all bound states (see Table 1). For example, the ( 4 He) 2 cluster with L = 1 dissociate to two 4 He atoms with different overall rotation quantum numbers (m = 0 and m = 1). It indicates the onset of t-symmetry breaking at infinite distance, which would have to be accounted for in dynamical studies of the helium motion in SWCNTs. This feature also correlates (at least for the 4 He isotope) with the fact of having (nearly degenerate) states that are anti-symmetric in t for N = 2 and odd L values. Indeed, for even L values, the dissociation channel is symmetric and the anti-symmetry in t is not present. States with odd L values are not symmetric in t and the dissociative WF is in fact a combination of AE states. For N = 3, all the dissociation channels are obviously non symmetric as combinations of the channels of the helium dimer and a single helium atom. However, antisymmetric in t states are not bound in this case.
As mentioned above, WFs bearing the same L value can be either symmetric or anti-symmetric with respect to the z -Àz exchange. The relatively large values for the splittings between AE states (see Table 1) come from the fact that the fragments wave functions (WFs) are very diffused. As a result, their overlaps are significant and contribute to the splittings. For instance, antisymmetric WFs should be zero at z = 0 in contrast with symmetric WFs. This additional restriction results in noticeable energy differences. The reason why the values of the splittings are large can be also understood as follows: the nanotube is large enough to hold two PNs with the same z value (i.e., t = 0). As the repulsive interaction is due to V 2 , which has the minimum (of just a few cm À1 ) at about 3 Å, such a configuration is just slightly penalized (if ever). It is additionally slightly penalized by the rotational kinetic energy term, due to an increased localization along the f coordinates. Thus, the barrier corresponding to the PN's exchange (or the end-over-end rotation) is found to be quite low. Actually, for N = 2, such t = 0 configuration corresponds to the ground state. This feature can be viewed making an analogy to a 1D model potential with two symmetric wells and a barrier. When the barrier is low enough, the splittings between symmetric and anti-symmetric solutions become large. In our case, it is difficult to plot an effective 1D V 2 potential along the z coordinate as all internal coordinates (t, w and, to a lesser extent, r) are strongly coupled.
Focusing now on the 2D densities ( Fig. 3-6), we observe that all systems present long tails in the t coordinates, though more compact than those of the free 4 He 2 dimer. All the 2D densities for the (PN) 2 CSWCNT(11,4) systems (see Fig. 3) indicate that the most probable structure is a dimer with the inter-nuclear axis perpendicular to the tube axis. The delocalization trend as the particle mass decreases and the strength of the interparticle interaction decreases is very pronounced. The analysis of the 2D densities of the (PN) 3 CSWCNT(11,4) systems is even more appealing (see Fig. 4). 43 For instance, it can be seen that, for PN = 4 He and ortho-D 2 , the WF bears a maximum at t 1,2 = 0, which can be attributed to a triangular structure. In stark contrast, for PN = 3 He and para-H 2 , the density along the t coordinate has an hexagonal-like shape which would correspond to an structure lying aligned along the tube axis. Actually, 4 He atoms also execute wide amplitude motion between triangular and aligned-type structures so that the corresponding 2D density in the w coordinates deviates noticeably from that of D 2 molecules. Hence, three D 2 molecules and 4 He atoms form a very flexible triangular structure (see also Fig. 5) while H 2 and 3 He atoms are delocalized along the tube axis. It is clear that the mass of the quantum particle is the dominant factor followed by the strength of the inter-particle interaction.

3.3
The special case of ( 4 He) 4 16. It can be observed that, whatever the bosonic particle be, the most probable structure is pyramidal-like. However, it features an increasing spatial delocalization when going from ortho-D 2 through para-H 2 to 4 He particles. On the one hand, although para-H 2 molecules are twice as light as 4 He atoms, the three times deeper attractive well of the pair potential (see Fig. 1) causes the structure of the (para-H 2 ) 4 cluster being more compact that the 4 He counterpart. In fact, the extended profile of the 2D densities in the w coordinates reflects a quasiindependent relative rotational motion of two ( 4 He) 2 dimers lying orthogonal to the tube axis. On the other hand, since ortho-D 2 molecules are twice as heavy as para-H 2 molecules, the solidlike nature of the (ortho-D 2 ) 4 structure (a regular tetrahedron) becomes even more apparent than the (para-H 2 ) 4 counterpart. Interestingly, (para-H 2 ) 4 clusters also adopt a (vibrationally averaged) tetrahedral geometry under spherical nanoconfinement in clathrate hydrate cages. 44,45 These similarities indicate that the aggregation of multiple H 2 molecules is driven by quantum nuclear effects and very weak dispersive H 2 -H 2 interactions rather than the much more stronger H 2 -SWCNT (or H 2 -cage) energy contribution.

and (para-H 2 ) 4 clusters
We have also investigated the stabilization of collective rotational motion. To this end, the lowest-state energies of the PN 4 CSWCNT (11,4) systems are represented at each L value (the so-called yrast states) in Fig. 6 (see also Table 2). As found for doped para-H 2 8 and confined ortho-D 2 16 clusters, four 4 He atoms or para-H 2 molecules are necessary for the appearance of the collective rotational states as energy minima such as that presented at L = N (see Fig. 7). The occurrence of energy minima in yrast lines at L = nN (n being an integer) has been connected with the onset of a superfluid behaviour and persistent flow in quantum rings made of 4 He atoms by Bloch. 46 Actually, bosonic symmetry is a necessary but not a sufficient condition.
For instance, the stabilization of collective rotational motion is not featured by (PN) N clusters with N o 4 (PN = 4 He, para-H 2 , see Table 1). As previously discussed, 8 such singularities appear when the hard-core H 2 -H 2 (or 4 He-4 He) repulsive interaction is intensified. This situation occurs, e.g., when a belt around a dopant molecule is close to completion with 4 He atoms or H 2 molecules. Such behaviour has been derived for bosonic particles featuring an infinity repulsive interaction when confined in a 1D ring. 47 However, we remark that our outcome is issued from highly accurate full quantum calculations and not a theoretical model. As for the case of ortho-D 2 molecules, 16 a classical analysis of the potential minima landscape leads to the prediction of a quasi-1D chain of pyramidal-like structures as the number of 4 He atoms and para-H 2 molecules increase (see Fig. 5). Because the well-depth of the V 1 potential is much larger than the V 2 counterpart, the 4 He atoms and the H 2 molecules can be considered as fixed at a distance r min from the nanotube axis, where r min is the position of the V 1 potential minimum. Then, the rest of the geometry is defined by V 2 . By requiring that all r ij distances have the values at which the V 2 potential minimum is located, the numbers of variables and conditions are the same for N = 4 so that the potential minimum is just a single point in Table 2 Energies of confined ( 4 He) 4 and (para-H 2 ) 4 clusters for different L values. The number of PO-DVR radial functions (n r ) is indicated as well as those of sinc-DVR functions in t and f coordinates (n t and n w ). For L 4 0, DE is the difference between the total energy and that corresponding to the ground-state with L = 0. The labels '+' and 'À' indicate that the wave function is symmetric (+) or anti-symmetric (À) with respect to the transformation z -Àz (referred to as z-symm) the parameter space. Considering the z-extension of the pyramidal structure for N = 4, we can conclude that the structure for N c 4 is a 1D chain of pyramidal-shaped structures (see Fig. 5), exhibiting a high delocalization when being formed by 4 He atoms.

Conclusions
This work addresses the characterization of multiple quantum particles inside carbon nanotubes and, in particular, in identifying the role of the nuclear spin statistics of the considered particles. It has been motivated by the experimental evidence of microscopic superfluidity (and its lack) in doped 4 He ( 3 He) droplets and para-H 2 clusters by the group of Professor Jan-Peter Toennies (see, e.g., ref. 1, 2, 5 and 6), but at the reduced dimensionality offered by confinement inside carbon nanotubes. For this purpose, we have developed a highly accurate full quantum method which is capable of dealing with Fermi-Dirac and Bose-Einstein nuclear spin statistics of the confined quantum particles on an equal footing, delivering excited and disentangled states with similar accuracy to the ground state, and using ab initio-derived potential models. In contrast with orbital-based techniques, we numerically solve the Schrödinger differential equation for the wave function with the full set of independent variables represented via a DVR grid so that the energy convergence has just to be ensure vs. the grid parameters and the grid size. In order to illustrate its performance, we have applied our new method to clusters of up to four (three) 4 He atoms and para-H 2 molecules ( 3 He atoms) inside a single-walled carbon nanotube of 1 nm diameter. The energy landscape of the bound states as a function of the angular momentum around the tube axis reveals very different profiles for clusters of 3 He and 4 He atoms. Due to the spin of the fermionic 3 He atom, the increasingly congested structure of the energy spectrum is evident as the number of 3 He atoms increases, in contrast with the spin-less (bosonic) 4 He counterpart. This finding is in line to that reported for doped 3 He clusters: 22,25 the calculated spectral profile of the dopant molecule becomes increasingly more congested upon augmentation of the number of 3 He atoms, explaining the experimental observations in doped 3 He droplets. 1 The wavefunction and density analysis of (PN) 3 clusters (PN = 3 He, 4 He, and para-H 2 ) indicates the occurrence of wide amplitude motions, specially for the 3 He isotope due its lighter mass and weakness of the He-He interaction. As for ortho-D 2 , 16 evidence of stabilization of collective rotational motion and persistent flow, 46 is found in clusters of four bosonic para-H 2 molecules and 4 He atoms, forming tetrahedral pyramidal-like structures. Remarkably, the main findings (arrangement in an unique (vibrationally averaged) tetrahedral geometry and exhibition of collective rotational motion) are shared with those uncovered for the case of spherical confinement in clathrate hydrate by Bačić and collaborators. 44,45 In our case of cylindrical confinement, the analysis of the potential minima landscape serves to predict a 1D chain of pyramidal-like structures along the tube axis as the resulting condensed matter system as the number of bosonic particles increases.
Pushing our understanding beyond common intuition, the present study along with previous works 16,17 point out that the transition from van der Waals-type molecular aggregation to quasi-1D condensed matter systems under cylindrical carbon nanotube confinement is driven by purely dispersive interadsorbate interactions together with quantum nuclear effects and not the much more stronger adsorbate-tube interaction forces. As a future perspective, we are aimed to adapt our Full Configuration Interaction Nuclear Orbital (FCI-NO) approach, originally developed for doped helium clusters, to the present case in which the carbon nanotube acts as the 'dopant' species. The essential advantage of this ansatz is that the WF is expanded using products of nuclear orbitals, allowing the employment of second quantization techniques and an automatic inclusion of Fermi-Dirac and Bose-Einstein nuclear spin statistical effects. This would imply the choice of an orbital representation such that the equality hZi = 0 is satisfied. For this purpose, one path is the imposition of the symmetry on the z coordinate for each orbital, and the follow-up subtraction of the center-of-mass kinetic energy. Whatever the specific strategy for further developments be, we hope that the highly accurate method presented in this work will serve to produce the necessary benchmark results.

Conflicts of interest
There are no conflicts to declare. as well as Andreas W. Hauser for his collaboration on previous works addressing confinement in carbon nanotubes.