Mini Review: Quantum Confinement of Atomic and Molecular Clusters in Carbon Nanotubes.

We overview our recent developments on a computational approach addressing quantum confinement of light atomic and molecular clusters (made of atomic helium and molecular hydrogen) in carbon nanotubes. We outline a multi-scale first-principles approach, based on density functional theory (DFT)-based symmetry-adapted perturbation theory, allowing an accurate characterization of the dispersion-dominated particle-nanotube interaction. Next, we describe a wave-function-based method, allowing rigorous fully coupled quantum calculations of the pseudo-nuclear bound states. The approach is illustrated by showing the transition from molecular aggregation to quasi-one-dimensional condensed matter systems of molecular deuterium and hydrogen as well as atomic 4He, as case studies. Finally, we present a perspective on future-oriented mixed approaches combining, e.g., orbital-free helium density functional theory (He-DFT), machine-learning parameterizations, with wave-function-based descriptions.

1 Introduction

The cylindrical confinement provided by carbon nanotubes has offered the possibility of studying the pronounced quantum behaviour of 4He atoms and H2 molecules at reduced dimensionality. Recent measurements have demonstrated the formation of two-dimensional (2D) 4He layers on the outer surface of single-walled carbon nanotubes (SWCNTs) (Noury et al., 2019). Also, an experimental study (Ohba, 2016) of gas adsorption at low (2–5 K) temperature revealed a quenched propagation of 4He atoms through carbon nanopores with diameters below 7 Å despite of their small kinetic diameter. The application of orbital-free helium density functional theory (He-DFT) to carbon nanotubes immersed in a helium nanodroplet provided theoretical explication that the experimental observations stem from the exceptionally high zero-point energy of 4He as well as its tendency to form two-dimensional (2D) layers upon adsorption at low temperatures (Hauser and de Lara-Castells, 2016). These conclusions were further confirmed by applying more accurate ab initio potential modelling along with a wave-function (WF)-based approach (Hauser et al., 2017). This study also showed that SWCNTs are filled by more molecules of N2 than 4He atoms, due to the higher zero-point energy of the latter. More generally, the interaction of small atomic and molecular systems with carbon-based nanostructures has attracted a lot of attention recently (see, e.g. Deb et al., 2016; Deb et al., 2019; Paul et al., 2019; Paul et al., 2020 and references therein).

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 (de Lara-Castells and Mitrushchenkov, 2020; de Lara-Castells and Mitrushchenkov, 2021), confirming a previous study using an embedding approach in a broader (ca. 1.4 nm) SWCNT (de Lara-Castells et al., 2017). Experimental evidence, using neutron scattering, on the formation of one-dimensional D2 crystals under carbon nanotube confinement has been reported as well (Cabrillo et al., 2021). Altogether, these studies have confirmed the key role played by the quantum nature of the nuclear degrees of freedom in the confined atomic helium, molecular hydrogen, or molecular deuterium motion.

The quest for the understanding of the aggregation of molecular H2 in carbon nanotubes is also application-oriented as it is being actively used as a clean energy source, substituting fossil fuels. In fact, its combustion produces only heat and water, and it can be efficiently combined with oxygen in a fuel cell to produce electricity. Yet, the usage of hydrogen as a profit fuel requires a substantial development of efficient storage materials (Schlapbach and Züttel, 2001; Züttel, 2003; Zheng et al., 2021) such as metal organic frameworks, Rosi et al. (2003) covalent organic frameworks (Zheng et al., 2021), and carbon-based nanoporous materials (Cheng et al., 2001; Liu et al., 2010). Hydrogen storage methods are aimed to pack hydrogen molecules as close as possible, existing direct experimental evidences for solid-like packing of hydrogen at temperatures of industrial importance (Ting et al., 2015). Actually, the question on the nature of the hydrogen packing in carbon nanotubes is truly fundamental. Thus, the possible existence of either a superfluid (Rossi and Ancilotto, 2016) or a crystalline phase (Del Maestro and Boninsegni, 2017; Ferré et al., 2017) for para-H2 molecules inside carbon nanotubes at zero (Rossi and Ancilotto, 2016; Ferré et al., 2017) or ultra-low temperatures (Del Maestro and Boninsegni, 2017) have been addressed as well (see also (Cazorla and Boronat, 2017) for a recent review). Specially appealing in regards to the impact of quantum effects on the diffusion of H2 and D2 along carbon nanotubes is the experimental finding of reserved trends in their rates upon cooling (Nguyen et al., 2010), with the D2 isotope becoming the faster inspite of its higher mass (Nguyen et al., 2010). The impact of quantum effects in confined H2 and D2 motion has been theoretically confirmed as well (Mondelo-Martell and Huarte-Larrañaga, 2016; Mondelo-Martell et al., 2017; Mondelo-Martell and Huarte-Larrañaga, 2021).

The high interest on the confinement of clusters of atomic 4He is also related to the superfluid nature of 4He droplets at a temperature close to absolute zero (0.4 K) (Mudrich and Stienkemeier, 2014) and, more specifically, to the onset of such fascinating property in doped clusters made of just a few 4He atoms (Toennies and Vilesov, 2004). Then, the question is whether the confinement provided by SWCNTs can favour the emergence of such feature at reduced dimensionality. The effect of carbon nanotube confinement has been also shown to be remarkable for the 3He isotope inside a SWCNT of 1 nm diameter. Thus, as opposed to the free dimer case, the confined 3He dimer is predicted to be bound, mainly as a consequence of the strong localization along the radial 3He–SWCNT direction (de Lara-Castells and Mitrushchenkov, 2021).

This overview article is organized as follows. The next section presents the computational approach addressing the modeling of intermolecular adsorbate–SWCNT interactions as well as fully coupled quantum calculations of the corresponding pseudo-nuclear wavefunctions. Section 3 presents an illustrative application of the computational approach to the transition from molecular aggregation to quasi-one-dimensional (1D) matter systems composed by either D2/H2 molecules or 4He atoms. Finally, Section 4 closes with a summary of concluding remarks and a few future prospects.

2 Computational Approach

2.1 Ab Initio Modelling of the Adsorbate-Nanotube Interaction

Achieving a correct description of the interaction of the lightest atomic and molecular clusters in nature (i.e., made of He or H2) with its confining SWCNT environment is a challenge even for expensive ab initio methods since they are dominated by long-range dispersion forces (van der Waals forces). Extended dispersion-corrected DFT methods are applicable including confinement effects in extended SWCNT structures but they bear a tendency to overshoot when applied to ultimate cases of dispersion-dominated (e.g., He-surface) interactions. One key idea has consisted in designing a functional [the so-called dlDF functional (Pernal et al., 2009)] which accounts for the dispersionless interaction energy only so that the dispersion corrections can be safely added later. Within a practical implementation of this idea, the incremental method (Stoll, 1992) is applied on non-periodic (small) cluster models of extended systems and combined with periodic dispersionless DFT calculations (de Lara-Castells et al., 2014b; de Lara-Castells et al., 2014a). This approach has been shown to be particularly successful when describing the interaction of atomic helium (de Lara-Castells et al., 2014a) as well as heavier noble gases (de Lara-Castells et al., 2015) with coronene/graphene/graphite surfaces. Very recently, it has been also demonstrated how modern ab initio intermolecular perturbation theory allows a cost-effective and accurate characterization of van der Waals-dominated He-SWCNT and H2-SWCNT interactions using small clusters models of the SWCNTs. The use of small cluster models is justified as complicated, dispersionless contributions are mostly of short-range nature. Dispersion, on the other hand, is long-range, but the corresponding parameters show excellent transferability properties upon increasing the size of the surface cluster models (de Lara-Castells et al., 2014a; de Lara-Castells et al., 2015). Therefore, these parameters can be calculated at high level of ab initio theory on small clusters and then scaled to the actual SWCNT system. As previously emphasized (de Lara-Castells and Hauser, 2020), detailed energy decomposition schemes, an intrinsic feature of methods such as density functional theory (DFT)-based symmetry-adapted perturbation theory [SAPT(DFT)] (Hesselmann and Jansen, 2003; Misquitta et al., 2003; Heßelmann et al., 2005; Misquitta et al., 2005) has been shown to be particularly useful in this respect.

In our most recent works (de Lara-Castells and Mitrushchenkov, 2020; de Lara-Castells and Mitrushchenkov, 2021), the SAPT(DFT) method has been used to derive by fitting the dispersionless and dispersion SAPT(DFT) energy contributions to an additive pairwise potential model (PPM) (de Lara-Castells et al., 2016; Hauser and de Lara-Castells, 2017; de Lara-Castells et al., 2017; Hauser et al., 2017), which is a modified version of that proposed by Carlos and Cole (Carlos and Cole, 1980). For the illustrative cases presented in this work, these terms were calculated for the interaction between a single H2 molecule (He atom) and a hydrogen-saturated (unsaturated) nanotube SWCNT(5, 5) tube made of 62 (40) atoms, considering two transverse sections of the tube. The SAPT(DFT) MOLPRO (Heßelmann et al., 2005; Werner et al., 2012) implementation was applied, using the Perdew-Burke-Ernzerhof (PBE) density functional (Perdew et al., 1996) and the augmented polarized correlation-consistent triple-zeta basis (Woon et al., 1994) for all atoms but the hydrogen atoms saturating the dangling bonds, for which the correlation-consistent double-zeta basis was used instead.

The PPM functional form (de Lara-Castells et al., 2016; Hauser and de Lara-Castells, 2017; de Lara-Castells et al., 2017; Hauser et al., 2017) for the dispersionless energy contribution accounts for the typical exponential growth of the dominant dispersionless term, the exchange-repulsion, but also including a Gaussian-type “cushion” to describe weakly attractive tails stemming from other dispersionless terms where R is a cut-off distance, R A-C stands for the distance between the adsorbate center-of-mass and one carbon atom of the SWCNT, and θ C is the angle between the radial vector going from the SWCNT center to one carbon atom and the vector R A-C pointing from the adsorbate center-of-mass to the same C atom. The dimensionless factor γ R in the first term accounts for the anisotropy of the C − C bonds. The sum in Eq. (1) runs over all carbon atoms of the nanotube. For the dispersion part, we apply the typical C6/C8 expansion with the damping functions of Tang and Toennies f (n = 6, 8) (Tang and Toennies, 1984) where γ A is also a dimensionless anisotropy parameter. The inclusion of γ and γ anisotropy terms has been found to be important when modelling potential corrugation effects on both curved carbon and metallic surfaces (de Lara-Castells et al., 2016; Hauser and de Lara-Castells, 2017; Hauser et al., 2017).

2.2 The Pseudo-Nuclear Wave-Function Problem

Our method has been designed to allow bound-state calculations for N identical atoms or molecules (referred to as “pseudo-nuclei”, PNs) inside a SWCNT. Using cylindrical coordinates (x , y , z → ρ , ϕ , z ) for each PN (see left panel of Figure 1), together with a Jacobian transformation of the volume to keep the Hamiltonian “explicitly Hermitian”, , the Hamiltonian describing a motion of N PN’s takes the form where the distance r is explicitly given by

FIGURE 1

Top panel. Left-hand side: Radial scan of the interaction energies between a single H2 molecule and a short (single-walled) carbon nanotube (sCNT) of helicity index (5, 5). The pairwise potential model is compared with reference ab initio calculations using the SAPT(DFT) approach. Right-hand side: Radial scan of the total interaction energies (full lines) and dispersion contributions (dashed lines) between a single H2 molecule and carbon nanotubes with helicity indexes of increasing value. Source: (de Lara-Castells et al., 2017) Bottom panel. Left-hand side: Cylindrical coordinates describing a single particle in a tube. Right-hand side: interaction potentials PN-CNT (A, B) and PN-PN (C, D) together with the corresponding low-lying bound states of a single confined PN. (A, C): PN = H2, (B, D): PN = 4He. Source: (de Lara-Castells and Mitrushchenkov, 2021).

We note that, in order to preserve the cylindrical symmetry of the system, the very small corrugation appearing along the azimuthal degree of freedom ϕ is not accounted for. Considering long nanotubes, it is thus assumed that the interaction potential depends only on ρ. The atomic structure of the SWCNT, however, is implicitly considered through the SAPT(DFT) calculations used to fit the parameters of our PPM (see Section 3.1). These potentials are used within the Discrete Variable Representation (DVR) approach (Bačic̀ and Light, 1989) upon fitting them to polynomials of ρ:

To separate the overall Z-translation and overall rotational motions, we introduce relative coordinates as t = z − z and χ = ϕ − ϕ , for i = 1 … N − 1. The overall Z-translation and overall rotation coordinates are conveniently defined as and Φ = ϕ . The overall rotation quantum number reads Λ = m 1 + m 2 + ⋯ + m , where m is the integer standing for the projection of the angular momentum of the i-th particle onto the nanotube axis. Since both PN–SWCNT and PN–PN interactions do not depend on Φ, the coordinate Φ can be separated via the factor . Finally, the full kinetic energy is written as where and

We note that the internal Hamiltonian does not account explicitly for the bosonic permutation symmetry of all N particles. The symmetry is automatically included for the particles labelled as 1 … N − 1 since it is equivalent to a simple exchange of the corresponding coordinates t , ρ , χ . The exchange of the particles labelled as i and N, however, results in linear transformations of the coordinates t and χ . The symmetry with respect to the i ↔ N exchange is analyzed a posteriori by defining a “bosonic symmetry factor” Q as the matrix element:

The Q factor is unity for true bosonic solutions. This test is used as a criteria to evaluate the wavefunction accuracy.

To calculate the eigenvalues of the internal Hamiltonian, we use the DVR approach (Bačic̀ and Light, 1989; de Lara-Castells and Mitrushchenkov, 2020; de Lara-Castells and Mitrushchenkov, 2021). This approach allows an easy and quick evaluation of the pair interaction potential V 2, with the basis set being obtained as a direct product of functions for the different coordinates. Sinc-DVR functions are conveniently employed for the t and χ coordinates. When dealing with the polar radii ρ, however, it is necessary to explicitly treat the singular kinetic energy term , which stems from the Jacobian transformation to cylindrical coordinates. For this purpose, we use the DVR basis obtained for the finite basis set representation (FBR) built from the radial functions of two-dimensional (2D) Harmonic oscillator functions. However, this basis renders the grid size too large for clusters with, e.g., 4 PNs since the internal Hamiltonian becomes 10-dimensional (10 D). This problem is efficiently solved by using potential-optimized DVR (PO-DVRs) functions (Echave and Clary, 1992), allowing for a very fast energy convergence as well. The diagonalization of the resulting Hamiltonian matrix is carried out using the Jacobi-Davidson algorithm (Sleijpen and van der Vorst, 1996). This technique has been very successful even for characterizing ill-behaved interactions such as, e.g., the “hard-core” interaction problem of doped helium clusters (de Lara-Castells et al., 2009a; de Lara-Castells et al., 2009b).

3 Illustrative Applications

In the first subsection, we will illustrate how the interaction of molecular hydrogen with various SWCNTs can be modelled at ab initio level. In the second subsection, we will show how our computational approach has allowed to reveal the transition from van-der-Waals-type molecular aggregation to quasi-1D condensed matter systems for atomic 4He and molecular H2 and D2 inside a SWCNT of 1 nm diameter.

3.1 H2-Nanotube Interaction Potentials for SWCNTs of Increasing Diameter

As an illustrative case, the top (left-hand) panel of Figure 1 shows the total H2/CNT interaction potential as a function of the radial distance r between the molecule center-of-mass and the SWCNT(5,5) center, along with the electrostatic Eelec, exchange-repulsion Eexch-rep, induction Eind, and dispersion contributions Edisp. Upon applying the pairwise PPM presented in Section 3.1, the radial scans of interaction energies shown in the upper (right-hand) panel of Figure 1 are obtained for SWCNTs of increasing diameter.

As can be observed in the top panel of Figure 1 (left-hand side), the H2/CNT attractive interaction is dispersion-dominated, with the exchange-repulsion dominating the whole dispersionless term. This repulsive energy contribution grows exponentially as the molecule-surface distance decreases although such behaviour is somewhat smoothed out by the attractive electrostatic contribution. The induction term contributes very little at the potential minimum but considerably at the SWCNT cage. The potential minimum is located at the center of the narrowest nanotube (diameter of 6.74 Å), where the adsorbates benefit from the dispersion interaction with carbon atoms at both sides of the carbon cage. However, upon increasing the nanotube diameter (see right-hand panel at the top of Figure 1), the dispersion becomes very small at the nanotube center and the potential minima shift towards the carbon cage. For a nanotube of diameter of 27.21 Å, the value of the well-depth (−56 meV) is already close to that obtained for a graphene sheet and measured for the H2 molecule adsorbed onto graphite (−51.7 meV from Ref. 29). From the top (left-hand side) panel of Figure 1 it can also be observed that the pairwise potential model closely reproduces the results obtained with the SAPT(DFT) method. Similar conclusions hold for the case of the He/SWCNT interaction (Hauser et al., 2017), also emphasizing the usefulness of the SAPT(DFT)-derived PPM (see Section 3.1).

3.2 From van-der-Waals Aggregation to Quasi-One-Dimensional Chains

Considering a SWCNT of 1 nm diameter and helicity index (11, 4), the He-SWCNT and H2-SWCNT one-particle potentials V 1(ρ) are shown in the bottom (right-side) panels of Figure 1 along with a few supported bound states of a single 4He atom or H2 molecule. Figure 2 presents the plots of the ground-state two-dimensional (2D) densities of PN4 clusters (PN = 4He, para-H2, and ortho-D2).

FIGURE 2

Plot of 2D densities for the ground state of PN4 ⊂ SWCNT(11, 4) complexes. Left to right: PN = 4He, para-H2, ortho-D2. (A) t 1 and t 2 coordinates; (B) t 1 and χ 1 coordinates; (C) t 1 and χ 2 coordinates; (D) χ 1 and χ 2 coordinates. Source: (de Lara-Castells and Mitrushchenkov, 2021).

As previously reported (de Lara-Castells and Mitrushchenkov, 2020; de Lara-Castells and Mitrushchenkov, 2021) the structuring of PN clusters, with N = 1 − 3, is characterized by a high pseudo-nuclear delocalization. Thus, the triangular-like structure exhibited by confined (4He)3 clusters (see graphical abstract) freely rotates inside the nanotube. However, as can be observed in Figure 2, whatever the bosonic particle be, the most probable structure is pyramidal-like. This confinement effect is specially remarkable for the case of 4He, due to the high delocalization of the bare PN–PN dimer (see bottom panel (d) of Figure 1). Yet, it features an increasing spatial delocalization when going from molecular ortho-D2 through molecular para-H2 to atomic 4He. Inspite of para-H2 molecules being twice as light as 4He atoms, the three times deeper attractive well of the pair potential (see bottom panels (a) and (b) of Figure 1) causes a much more compact structuring of the (para-H2)4 cluster when compared with the 4He counterpart. Particularly, we note that the extended profile of the 2D densities in the χ coordinates reflects a quasi-independent relative rotational motion of two (4He)2 dimers lying orthogonal to the tube axis. Since ortho-D2 molecules are twice as heavy as para-H2 molecules, the solid-like nature of the (ortho-D2)4 structure (a regular tetrahedron) becomes even more apparent than for the (para-H2)4 counterpart. The formation of (vibrationally averaged) tetrahedral structures of (para-H2)4 clusters is not a special effect from the SWCNT confinement. In fact, the same feature has also been revealed in clathrate hydrate cages (Sebastianelli et al., 2008; Witt et al., 2010). These similarities indicate that the aggregation of multiple light PNs are driven by quantum nuclear effects and very weak dispersive PN–PN interaction forces rather than the much more stronger PN–SWCNT (or PNT-cage) energy contribution (see bottom panels (a) and (b) of Figure 1). As discussed in detail in our previous work (de Lara-Castells and Mitrushchenkov, 2020), the onset of solid-like packing for molecular deuterium is explained by analyzing the potential minima landscape, allowing also to predict the formation of a one-dimensional chain of tetrahedral structures along the tube axis when the number of D2 molecules increase. Similar conclusions holds for the cases of 4He atoms and H2 molecules (de Lara-Castells and Mitrushchenkov, 2021). Moreover, these special structures feature stabilization of collective rotational motion resembling the behaviour of quantum rings exhibiting persistent current (charged particles) or persistent flow (neutral particles). This characteristic has been connected with the onset of superfluid motion and persistent flow in quantum rings made of 4He atoms by Bloch (Bloch, 1973). It is remarkable that such feature occurs already in clusters made of four 4He atoms at the reduced dimensionality offered by the confining medium.

4 Concluding Remarks and Future Directions

Summarizing, this mini review shows how a challenging case of confined molecular system can be accurately characterized using multi-scale first-principles modelling. In order to emphasize the key role of quantum effects, we have chosen the lightest atomic and molecular species in nature as the confined object and SWCNTs as the confining medium. In the first part, we have overviewed how an ab initio-derived pairwise potential model can accurately characterize the intermolecular interaction between the confining medium and the confined object. In the second part, we have illustrated how a computational approach, allowing for rigorous fully coupled quantum calculations of the pseudo-nuclear bound states, has been capable of revealing the transition from van der Waals-type molecular aggregation to quasi-1D condensed matter systems under cylindrical carbon nanotube confinement. Remarkably, it has been also shown that the structuring is driven by purely dispersive PN-PN interactions together with quantum pseudo-nuclear effects and not the much more stronger PN-tube interaction forces. In this way, it can be also understood why, inspite of its weak cohesive PN-PN interaction, neutron scattering experimental measurements have very recently shown the formation of quasi-one-dimensional structures of molecular D2 corresponding to cylindrical sections of the hexagonal close-packed bulk crystal (Cabrillo et al., 2021).

As a future perspective, we are aimed to adapt our Full Configuration Interaction Nuclear Orbital (FCI-NO) approach (de Lara-Castells et al., 2009b), 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 wave function is expanded using products of nuclear orbitals, allowing the employment of second quantization techniques, an automatic inclusion of Fermi–Dirac and Bose–Einstein nuclear spin statistical effects, and, then, a characterization of 3He and 4He isotopes on an equal footing. The second future prospect is the creation of embedding approaches linking our wave-function-based description with orbital free He-DFT (Dalfovo et al., 1995; Ancilotto et al., 2017), following previous efforts (Hauser et al., 2017). In this sense, the highly accurate method presented in our latest work (de Lara-Castells and Mitrushchenkov, 2021) is expected to provide benchmark results guiding new, machine-learning driven, parameterizations in He-DFT, as already illustrated in electronic structure theory (Meyer et al., 2020).