(Ro)vibrational Spectroscopic Constants, Lifetime and QTAIM Evaluation of Fullerene Dimers Stability

The iconic caged shape of fullerenes gives rise to a series of unique chemical and physical properties; hence a deeper understanding of the attractive and repulsive forces between two buckyballs can bring detrimental information about the structural stability of such complexes, providing significant data applicable for several studies. The potential energy curves for the interaction of multiple van der Waals buckyball complexes with increasing mass were theoretically obtained within the DFT framework at ωB97xD/6−31G(d) compound model. These potential energy curves were employed to estimate the spectroscopic constants and the lifetime of the fullerene complexes with the Discrete Variable Representation and with the Dunham approaches. It was revealed that both methods are compatible in determining the rovibrational structure of the dimers and that they are genuinely stable, i.e., long-lived complexes. To further inquire into the nature of such interaction, Bader’s QTAIM approach was applied. QTAIM descriptors indicate that the interactions of these closed-shell systems are dominated by weak van der Waals forces. This non-covalent interaction character was confirmed by the RDG analysis scheme. Indirectly, QTAIM also allowed us to confirm the stability of the non-covalent bonded fullerene dimers. Our lifetime calculations have shown that the studied dimers are stable for more than 1 ps, which increases accordingly with the number of carbon atoms.


Introduction
Fullerenes, also known as buckyballs, are hollow polyhedral molecules constituted of sp2-hybridized covalently bonded carbon atoms, with the faces made up of a combination between pentagons and hexagons [1][2][3]. They were originally predicted by Osawa in 1970 [1] and first synthesized in 1985 by Kroto and co-workers [2]. The iconic caged shape of buckyballs gives rise to interesting chemical properties, such as trapping small molecules inside fullerenes [4,5], organometallic fullerenes acting as hydrogen adsorbents [6], and endohedral doping with metals resulting in high electrical conductivity [7]. The unique physical and chemical properties of fullerenes attracted a lot of attention aimed at their possible applications to organic electronics. For example, as components of photovoltaic cells [8,9], for the development of sensors employed for the capture of pollutants [10,11], or for electron transport through self-assembled monolayers of functionalized fullerenes [12,13]. In this context, information about the intermolecular potential between fullerene molecules, the equilibrium distance, and energy minima of the van der Waals contacts are of great importance to understanding these phenomena [14].
Dynamic force microscopy associated with dispersion-corrected density functional theory (DFT) was used to map the intermolecular potential of the C 60 fullerene by Sweetman and co-workers, showing that the different orientations of the C 60 introduced a variation from~50 to 60 meV in the potential minimum [15]. These authors also showed that there is a subtle interplay between attractive and repulsive forces governing the intermolecular potential, in which variations of the dispersion interaction are screened by repulsive interaction even at large distances.
In one of our previous manuscripts, we investigated the influence that the different orientations of dimers of C 70 have on the potential energy curve and in the spectroscopic constants [16]. It was shown that different configurations might be accessed by a fullerene dimer, which may change the spectroscopic and stability characteristics of those molecules. Furthermore, in another of our recent research [11], the importance of non-covalent interactions (NCI) was investigated to better understand the characteristics of the intermolecular interaction between the C 59 X doped fullerene (X = B, Al, Ga, Si, Ge, N, and P) and the emerging pollutant Carbamazepine CBZ. In this research, it was shown that intermediate stable states could also be accessed due to a small change in the configuration of the interaction, allowing the discovery of new active sites for CBZ adsorption on C 59 X fullerenes.
Therefore, a deeper understanding of these features concerning NCI, as well as the spectroscopic properties and structural stability of the van der Waals bonds provided by theoretical data, are of great significance when studying supramolecular fullerenes for several applications. Baggio and co-workers [17] have reported a theoretical study of the binding energies between ammonia and metallo-phthalocyanines, in which the rovibrational spectroscopic constants were obtained by approximating both species as a huge "monoatomic" compound. Such an approach neglects all intramolecular interactions, thus greatly simplifying the solution of the nuclear Schrödinger equation while reporting good correlation with experimental data, characterizing the spontaneity of ammonia/metallophthalocyanines complexes at temperatures around 500 K. The same protocol was applied, by the same authors, to predict the thermodynamics and stability of methanol noble gasses complexes [18].
This work aims at calculating the potential energy curves (PECs) of numerous fullerene dimers of different sizes. Specifically, we are interested in studying the stability and determining the spectroscopic properties of (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 . All dimers employed in this research are presented in Figure 1. Spectroscopic properties constitute an interesting tool to assess the stability of the dimer form because long-lived complexes tend to oscillate around the equilibrium distance, and therefore, quantized (ro)vibrational levels are expected to exist in the potential well. To achieve this goal, we employed the same methodology applied in the calculation of spectroscopic properties of diatomic systems [16,[19][20][21] to the case of fullerene dimers, as will be further discussed in this article. The stability of these fullerene dimers was also investigated through lifetime calculations, while the quantum theory of atoms in molecules (QTAIM) [22] was employed to determine the dimers' type of intermolecular interaction.

Results and Discussion
Let us begin our discussion by presenting the PECs obtained for the (C84)2 fullerene dimer in Figure 2 and for the remaining dimers in Figures S1-S6. As discussed in the Methodology section, the ab initio ωB97xD/6−31G(d) points (blue dots in Figure 2) were fitted through the Ryd6 analytical function (blue line in Figure 2). The quantitative pointby-point difference between ab initio energies and fitted Ryd6 PEC accounts for the error in the fitting procedure (red line in Figure 2). The optimized Ryd6 parameters for each fullerene dimer are reserved in Tables S1 and S2 to avoid the proliferation of the results. Notice from Figure 2 that the fitting errors do not surpass 0.05 kcal mol −1 for the (C84)2 dimer. The same fitting accuracy was transferable among the other fullerene dimers, indicating that the optimization of the Ryd6 parameters is consistent and reliable, as confirmed by the global accuracy of the χ 2 values in the range of 10 −5 kcal mol −1 to 10 −4 kcal mol −1 (see Tables S1 and S2). is also plotted to magnify the overall fit accuracy in the entire range of internuclear separation.

Results and Discussion
Let us begin our discussion by presenting the PECs obtained for the (C 84 ) 2 fullerene dimer in Figure 2 and for the remaining dimers in Figures S1-S6. As discussed in the Methodology section, the ab initio ωB97xD/6−31G(d) points (blue dots in Figure 2) were fitted through the Ryd6 analytical function (blue line in Figure 2). The quantitative point-bypoint difference between ab initio energies and fitted Ryd6 PEC accounts for the error in the fitting procedure (red line in Figure 2). The optimized Ryd6 parameters for each fullerene dimer are reserved in Tables S1 and S2 to avoid the proliferation of the results. Notice from Figure 2 that the fitting errors do not surpass 0.05 kcal mol −1 for the (C 84 ) 2 dimer. The same fitting accuracy was transferable among the other fullerene dimers, indicating that the optimization of the Ryd6 parameters is consistent and reliable, as confirmed by the global accuracy of the χ 2 values in the range of 10 −5 kcal mol −1 to 10 −4 kcal mol −1 (see Tables S1 and S2).

Results and Discussion
Let us begin our discussion by presenting the PECs obtained for the (C84)2 fullerene dimer in Figure 2 and for the remaining dimers in Figures S1-S6. As discussed in the Methodology section, the ab initio ωB97xD/6−31G(d) points (blue dots in Figure 2) were fitted through the Ryd6 analytical function (blue line in Figure 2). The quantitative pointby-point difference between ab initio energies and fitted Ryd6 PEC accounts for the error in the fitting procedure (red line in Figure 2). The optimized Ryd6 parameters for each fullerene dimer are reserved in Tables S1 and S2 to avoid the proliferation of the results. Notice from Figure 2 that the fitting errors do not surpass 0.05 kcal mol −1 for the (C84)2 dimer. The same fitting accuracy was transferable among the other fullerene dimers, indicating that the optimization of the Ryd6 parameters is consistent and reliable, as confirmed by the global accuracy of the χ 2 values in the range of 10 −5 kcal mol −1 to 10 −4 kcal mol −1 (see Tables S1 and S2). is also plotted to magnify the overall fit accuracy in the entire range of internuclear separation. is also plotted to magnify the overall fit accuracy in the entire range of internuclear separation.
Other published analytical potentials were also evaluated for the (C 60 ) 2 and compared to the Ryd6 potential. The intermolecular potential forms of Girifalco, Smith-Thakkar, and Lim were fitted to the ωB97xD/6−31G(d) points (Figure 3), and their corresponding expressions were reserved in Table S5 along with their optimized parameters and squared correlation effects. Originally, these analytic forms were proposed to fit bulk data obtained from crystals, neglecting deformation and orientation effects [23]. Although deformation can be neglected in such types of van der Waals complexes, orientation effects should play an important role in the intermolecular interaction of carbon nanostructures [16]. Furthermore, we see in Figure 3 that the potential forms of Girifalco, Smith-Thakkar, and Lim were unable to describe the PEC at the equilibrium region. Since the NSE depends on a good description of the PEC, these potential forms were not effectively suitable to describe the spectroscopic constants for the fullerene dimers studied herein. Other published analytical potentials were also evaluated for the (C60)2 and compared to the Ryd6 potential. The intermolecular potential forms of Girifalco, Smith-Thakkar, and Lim were fitted to the ωB97xD/6−31G(d) points (Figure 3), and their corresponding expressions were reserved in Table S5 along with their optimized parameters and squared correlation effects. Originally, these analytic forms were proposed to fit bulk data obtained from crystals, neglecting deformation and orientation effects [23]. Although deformation can be neglected in such types of van der Waals complexes, orientation effects should play an important role in the intermolecular interaction of carbon nanostructures [16]. Furthermore, we see in Figure 3 that the potential forms of Girifalco, Smith-Thakkar, and Lim were unable to describe the PEC at the equilibrium region. Since the NSE depends on a good description of the PEC, these potential forms were not effectively suitable to describe the spectroscopic constants for the fullerene dimers studied herein. To better appreciate the overall aspects of the obtained PECs, we present in Figure 4 the superimposed intermolecular potential for all dimers studied in this work. Here, we are eager to explore the size effects on the potential energy minimum (De) and its separation (Re) explicated in Table 1, along with other results obtained with different methods for the (C60)2 and the (C70)2 dimers. From Figure 4 and Table 1, it is noticed that as the number of carbon atoms in each monomer increases, the intermolecular interaction is intensified, and the potential well becomes more profound. To better appreciate the overall aspects of the obtained PECs, we present in Figure 4 the superimposed intermolecular potential for all dimers studied in this work. Here, we are eager to explore the size effects on the potential energy minimum (D e ) and its separation (R e ) explicated in Table 1, along with other results obtained with different methods for the (C 60 ) 2 and the (C 70 ) 2 dimers. From Figure 4 and Table 1, it is noticed that as the number of carbon atoms in each monomer increases, the intermolecular interaction is intensified, and the potential well becomes more profound.     [24] had the values estimated for covalently bonded dimers, [25,26] for free dimers, [27,28] obtained values for dimers in a cluster of fullerenes.

Fullerene Dimers
We now focus on the results obtained by solving the NSE (Equation (4)). The values of the spectroscopic constants obtained using the Dunham and DVR methods are presented in Table 2. There was a satisfactory agreement between the values of the spectroscopic constants through both methods. The rotational constant Be and the rotorvibrating coupling constant αe calculated for the dimers were negligibly small due to the relatively large reduced mass for the studied systems. The reduced mass μ values were taken to be 120.11 , 144.13 u, 216.19 u, 360.33 u, 420.38 u, and 504.46 u for the (C20)2, (C24)2, (C36)2, (C60)2, (C70)2, and (C84)2 dimer, respectively. Concerning the fundamental vibrational frequency ωe, we can see from Table 2 that there was a small blue shift (less than ~10 cm −1 ) when increasing the size of the dimers. By analyzing the anharmonicity effects measured by the ωexe constant, we notice that it becomes more relevant for the smallest dimers, representing one-tenth of the harmonic frequency for the (C20)2 and (C24)2 dimers.  We now focus on the results obtained by solving the NSE (Equation (4)). The values of the spectroscopic constants obtained using the Dunham and DVR methods are presented in Table 2. There was a satisfactory agreement between the values of the spectroscopic constants through both methods. The rotational constant B e and the rotor-vibrating coupling constant α e calculated for the dimers were negligibly small due to the relatively large reduced mass for the studied systems. The reduced mass µ values were taken to be 120.11 u, 144.13 u, 216.19 u, 360.33 u, 420.38 u, and 504.46 u for the (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 dimer, respectively. Concerning the fundamental vibrational frequency ω e , we can see from Table 2 that there was a small blue shift (less than~10 cm −1 ) when increasing the size of the dimers. By analyzing the anharmonicity effects measured by the ω e x e constant, we notice that it becomes more relevant for the smallest dimers, representing one-tenth of the harmonic frequency for the (C 20 ) 2 and (C 24 ) 2 dimers. The vibrational eigenvalues (with J = 0) of the NSE obtained through the DVR method are presented in Table 3 for the lowest-lying vibrational states of the fullerene dimers.
Because of the small rotational effects on the interaction mode (B e~1 0 −4 cm −1 and α e~1 0 −7 cm −1 ), the dimers are essentially a rigid rotator system; therefore, the values for ε υ,J with J = 0 were suppressed in Table 3. Additionally, we estimated the maximum number of vibrational levels harbored by the PECs for each fullerene dimer, i.e., we determined the largest value of υ in ε υ,J satisfying the inequality ε υ,J < D e . It turns out that the number of bound vibrational states for each dimer was 60, 70, 146, 222, 237, and 265 for the (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 systems, respectively. Table 3. Quantized energies of the first ten vibrational levels (in cm −1 ) for the (C 20 With the rovibrational constants in hand, we can also calculate the vibrational population for each state and the temperature effects of the Boltzmann distribution of the quantized levels, which are displayed in Figure 5. The calculated spectroscopic constants were inserted in the rovibrational partition function that can be found in ref. [17]. To confirm the stability of these dimers, a lifetime calculation was performed using the De and ωe parameters occurring in Equation (11). The lifetimes of the van der Waals complexes as a function of temperature are presented in Figure 6. As a general rule, when a given complex presents lifetimes larger than 1.0 ps, it is deemed to be a stable one; on In Figure 5, we first focused on the cases of the strongest interacting systems, i.e., (C 70 ) 2 and (C 84 ) 2 where we can see the anharmonicity of the vibrational structure. For instance, at room temperature (300 K), we see that, even for those strongest interacting systems, less than 1 4 of the dimers lie on the vibrational ground state and that the inclusion of higher vibrational excited states is mandatory for a good description of the normal modes. The harmonic oscillator accurately models these systems only at very low temperatures of 10 K. Therefore, once again, it highlighted the importance of using adjustable potential functions, such as Ryd6, to adjust the potential energy curve employed in NSE to achieve a suitable description of the system's spectroscopic constants.
To confirm the stability of these dimers, a lifetime calculation was performed using the D e and ω e parameters occurring in Equation (11). The lifetimes of the van der Waals complexes as a function of temperature are presented in Figure 6. As a general rule, when a given complex presents lifetimes larger than 1.0 ps, it is deemed to be a stable one; on the other hand, if the value is inferior to this threshold value, the complex is considered transitory [29]. Figure 6 reveals that all fullerene dimers studied here present lifetimes that are larger than 1.0 ps, even when the temperature is raised to normal conditions. The numerical values are reserved in Table S4 for temperatures ranging from 200 to 500 K. The inset plot in Figure 6 shows that even the weakest complexes (C 20 ) 2 and (C 24 ) 2 are of the order of τ~1 ps, thus relatively stable. Figure 6 also confirms that the lifetime values have the same trends observed for the interaction energies, D e , and therefore (C 84 ) 2 is the most stable system, followed by (C 70 ) 2 , (C 60 ) 2 , (C 36 ) 2 , (C 24 ) 2, and (C 20 ) 2 . Understanding the nature of the NCIs occurring in van der Waals complexes is of primary interest to assess the stability of the dimers and to gain insight into the physical nature of the interaction. To guarantee accurate NCI description, the fullerene dimers were optimized at the equilibrium distance obtained from the PECs calculations. Although the rigid scan showed acceptable energies, a dimer structure optimized at the minimum energy should be considered, especially for the analysis of the interaction type between the monomers. The QTAIM approach developed by Bader [22] is a qualitative quantum chemical tool to evaluate the NCIs from an analysis of the topological properties of the electron density. Within the QTAIM methodology, the so-called bond paths describe a line of local maximum electron density connecting two nuclear critical points (NCPs, i.e., the actual position of the atoms) that interacts either covalently or noncovalently. For this purpose, the electron density ( ) and its Hessian matrix (∇ = + + ) were analyzed at the region between two interacting molecules. The algebraic signs of the eigenvalues ( = 1, 2,3) of the Hessian matrix distinguish the nature of the critical points. At the nuclei, all eigenvalues are negatives. At this region, the density exhibits local maximum, the so-called nuclear critical point (NCP, i.e., the atoms' actual Understanding the nature of the NCIs occurring in van der Waals complexes is of primary interest to assess the stability of the dimers and to gain insight into the physical nature of the interaction. To guarantee accurate NCI description, the fullerene dimers were optimized at the equilibrium distance obtained from the PECs calculations. Although the rigid scan showed acceptable energies, a dimer structure optimized at the minimum energy should be considered, especially for the analysis of the interaction type between the monomers. The QTAIM approach developed by Bader [22] is a qualitative quantum chemical tool to evaluate the NCIs from an analysis of the topological properties of the electron density. Within the QTAIM methodology, the so-called bond paths describe a line of local maximum electron density connecting two nuclear critical points (NCPs, i.e., the actual position of the atoms) that interacts either covalently or non-covalently. For this purpose, the electron density (ρ ) and its Hessian matrix (∇ 2 ρ = λ 1 + λ 2 + λ 3 ) were analyzed at the region between two interacting molecules. The algebraic signs of the eigenvalues λ i (i = 1, 2, 3) of the Hessian matrix distinguish the nature of the critical points. At the nuclei, all eigenvalues are negatives. At this region, the density exhibits local maximum, the so-called nuclear critical point (NCP, i.e., the atoms' actual position). When two eigenvalues of the Hessian matrix are negative ( λ 1 , λ 2 ), we observed the formation of a bond critical point (BCP). When only one value of the Hessian matrix is negative ( λ 1 ), it is observed as a ring-critical point (RCP). Generally, RCP are points of steric effects. When none of the eigenvalues are negative, the formation of a local minimum in the density, named cage critical point (CCP), is observed. At the BCP, the values of ρ and ∇ 2 ρ are closely related to the nature of the interaction. Figure 7 portrays the molecular graphs of the interacting fullerene dimers separated at the minimum of the well showing the BCPs, RCPs, CCPs, and bond paths. Interestingly, for those fullerenes that were more spherically symmetric, we observed fewer CPs and bond paths connecting the two (C n ) 2 monomers.  Table S6 presents the values of the Electron density ρcp (in e/a0 3 ), the Laplacian of the electron density ∇ 2 ρcp (in e/a0 5 ) and the total energy HBCP (in Hartree).
To facilitate the view and comparation with the results shown the Table S6, the CP where numbered. The QTAIM analysis was performed for the optimized dimers at the equilibrium position.
To facilitate the analyses and interpretation of the NCIs, QTAIM offers a relatively simple and straightforward interpretation of the intermolecular forces by calculating parameters at the BCP, such as the electron density (ρBCP), its Laplacian (∇ 2 ρBCP) and the total electron energy density (HBCP), which are summarized in Table S6. The nature of the interacting BCP can be categorized depending on the signal and magnitude of ∇ 2 ρBCP and HBCP. When ∇ 2 ρBCP < 0 and HBCP < 0, it is a good indicator of a very strong interaction typical of open-shell interactions (covalent bonds). On the other hand, when ∇ 2 ρBCP > 0 and HBCP Figure 7. QTAIM critical points for the fullerene dimers. The orange dots represent the BCPs between the molecules, the yellow dots represent the critical points of the RCP ring, and the green dots represent the cage critical points of CCP. Table S6 presents the values of the Electron density ρ cp (in e/a 0 3 ), the Laplacian of the electron density ∇ 2 ρ cp (in e/a 0 5 ) and the total energy H BCP (in Hartree).
To facilitate the view and comparation with the results shown the Table S6, the CP where numbered. The QTAIM analysis was performed for the optimized dimers at the equilibrium position.
To facilitate the analyses and interpretation of the NCIs, QTAIM offers a relatively simple and straightforward interpretation of the intermolecular forces by calculating parameters at the BCP, such as the electron density (ρ BCP ), its Laplacian (∇ 2 ρ BCP ) and the total electron energy density (H BCP ), which are summarized in Table S6. The nature of the interacting BCP can be categorized depending on the signal and magnitude of ∇ 2 ρ BCP and H BCP . When ∇ 2 ρ BCP < 0 and H BCP < 0, it is a good indicator of a very strong interaction typical of open-shell interactions (covalent bonds). On the other hand, when ∇ 2 ρ BCP > 0 and H BCP > 0, the BCP along the bond path is commonly associated with closed-shell interactions (NCIs) when accompanied by a small value of ρ BCP [30][31][32].
From Table S6, we observe small values of ρ BCP at the BCPs between the monomers of the order of 10 −3 e a 0 −3 , and associated with the positive values of ∇ 2 ρ BCP and H BCP. This confirms, from a QTAIM standpoint, that the intermolecular interactions can be weak closed-shell interactions, i.e., they attract each other through van der Waals forces.
In association with QTAIM, the analysis of the reduced density gradient (RDG) aids the interpretation of the nature of the NCIs. RDG introduces a generalization of the approach focused on critical points. This allows the production of isosurfaces that yield additional evidence for non-covalent interactions. Figure 8 presents the diagrams of isosurfaces for NCIs present in the fullerene dimers. In Figure 8, green-colored isosurfaces show non-covalent interactions, while red-colored isosurfaces present steric effects.
From Table S6, we observe small values of ρBCP at the BCPs between the monomers of the order of 10 , and associated with the positive values of ∇ 2 ρBCP and HBCP. This confirms, from a QTAIM standpoint, that the intermolecular interactions can be weak closed-shell interactions, i.e., they attract each other through van der Waals forces.
In association with QTAIM, the analysis of the reduced density gradient (RDG) aids the interpretation of the nature of the NCIs. RDG introduces a generalization of the approach focused on critical points. This allows the production of isosurfaces that yield additional evidence for non-covalent interactions. Figure 8 presents the diagrams of isosurfaces for NCIs present in the fullerene dimers. In Figure 8, green-colored isosurfaces show non-covalent interactions, while red-colored isosurfaces present steric effects. In the RDG analysis, while the behavior of indicates the critical points in a 2D graph, the combination of electron density and the signal of provides a tool to distinguish between repulsive or attractive interactions. For < 0 it is reported as an attractive interaction, while for repulsive interaction, > 0 [32]. For this reason, a real space function, ( ) , namely the product of with is employed to help with this task of identifying the type of the interaction. In summary, large negative values of ( ) are indicative of strong attractive interactions, such as hydrogen bounds or dipole-dipole interactions. Conversely, large positive values of ( ) represents non-bonding interactions, and when ( ) ≈ 0, it is indicative of van der Waals interactions [30,32]. Figure 9 presents the scatter graph of ( ) versus RDG of non-covalent interaction for fullerenes dimers. In the graph, the regions where the values of the function ( ) are between −0.010 and 0.010, reveal van der Walls interactions. Those regions are the spikes located on the X-axis in the graph and are represented by the green and In the RDG analysis, while the behavior of s indicates the critical points in a 2D graph, the combination of electron density ρ and the signal of λ 2 provides a tool to distinguish between repulsive or attractive interactions. For λ 2 < 0 it is reported as an attractive interaction, while for repulsive interaction, λ 2 > 0 [32]. For this reason, a real space function, sign(λ 2 )ρ, namely the product of λ 2 with ρ is employed to help with this task of identifying the type of the interaction. In summary, large negative values of sign(λ 2 )ρ are indicative of strong attractive interactions, such as hydrogen bounds or dipole-dipole interactions. Conversely, large positive values of sign(λ 2 )ρ represents non-bonding interactions, and when sign(λ 2 )ρ ≈ 0, it is indicative of van der Waals interactions [30,32]. Figure 9 presents the scatter graph of sign(λ 2 )ρ versus RDG of non-covalent interaction for fullerenes dimers. In the graph, the regions where the values of the function sign(λ 2 )ρ are between −0.010 and 0.010, reveal van der Walls interactions. Those regions are the spikes located on the X-axis in the graph and are represented by the green and light brown isosurfaces shown in Figure 8. Regions where sign(λ 2 )ρ assumes values over 0.020 are regions of steric effects. Those regions are represented by the red surfaces in Figure 8.  QTAIM results and RDG index showed electron density depletion in the intermolecular region. The analyses of the electronic density, the Laplacian of the electronic density, and the total energy at the critical points in association with the RDG isosurfaces and with the spikes observed in the scattering graphs, presented in Figures 8  and 9 confirm the non-covalent interaction as the responsible for the fullerene dimer QTAIM results and RDG index showed electron density depletion in the intermolecular region. The analyses of the electronic density, the Laplacian of the electronic density, and the total energy at the critical points in association with the RDG isosurfaces and with the spikes observed in the scattering graphs, presented in Figures 8 and 9 confirm the non-covalent interaction as the responsible for the fullerene dimer stabilization. This van der Waals nature of the NCIs correlates very well with the values of the potential depth obtained by our DFT estimates.
In addition, QTAIM can be employed as an indirect tool to evaluate the stability of weak interacting molecules [33][34][35][36]. When an unstable interaction is analyzed under the QTAIM schemes, the presence of critical degenerate points can be observed in the intermolecular region. The closer the RCP and the BCP are, the more unstable the interaction tends to be. In unstable intermolecular interactions, a small energy disturbance can cause the breaking of an RCP to form a BCP [22]. As can be seen in the QTAIM parameters, shown in Table S6, no degenerated critical point was observed. This way, the QTAIM analyses deal with energetic and lifetime evaluation, confirming the fullerene dimers' stability, in particular for the (C 70 ) 2 and (C 84 ) 2 systems which showed a well-defined face-to-face interaction, presented for stable and direct bond paths in the intermolecular region.
Finally, the NCI analyses provided by QTAIM and RDG approach allow us to verify the interacting more favorable face for each fullerene dimer as well the trend of the bond. It was observed by QTAIM analysis that the dimers of C 36 , C 70, and C 84 fullerenes are prone to interact through their hexagonal-hexagonal faces. The RDG isosurfaces confirm this trend through its hexahedron form. The (C 60 ) 2 dimers presented an interaction via one pentagonal face oriented to an inter-pentagon bonding (see Table S8 in supporting information and Figure S7). Although this configuration differs from the more energetically stable configuration observed by Sharapa and co-authors' work [23], it is notable that our calculations returned that these pentagonal faces are oriented towards the configuration of inter-pentagon connections. Similar inter-pentagon relative orientation has been observed for bulk C 60 fullerenes at low temperatures [18].

Methods
Intermolecular Interaction Potential. To determine the PECs for all systems, we performed a rigid scan on the interaction coordinate, i.e., keeping one of the monomers stationary and then moving away the second monomer with incremental steps of 0.1 Å with respect to the center-to-center distance, as shown in Figure 1A. The binding energies (E Binding ) were estimated by means of the supermolecular approach given by [37], where E dimer is the energy of the fullerene dimer, E mon,A and E mon,B represent the energy of monomers A and B, respectively. Figure 1B shows the relative orientations of each fullerene dimer used to build the targeted PECs. The structures for each fullerene shown in Figure 1B were retrieved from the C n fullerenes database [38]. We selected the C 20 (I h ), C 24 (D 6d ), C 36 (D 6h ), C 60 (I h ), C 70 (D 5h ), and C 84 (D 2 ) fullerenes because they presented the lowest energy in the ground state and the lowest strain energy being considered as the more stable isomers [39,40]. Electronic structure calculations were performed using DFT at the ωB97xD/6−31G(d) chemical model [23,24], as implemented in the Gaussian 09 package [41]. The rangeseparated hybrid functional ωB97xD was chosen because it is dispersion-corrected through Grimme's GD2 empirical dispersion model [42] and has shown to be very reliable in capturing the binding energy of the C 60 dimer when Grimme's empirical dispersion model GD3 is used [23]. Grimme's GD2 model has an extra energy term that includes the pair-wise London dispersion interaction between the monomers and uses a similar damping function to that used by the GD3 model. Basis Set Superposition Error (BSSE) correction was also included through the counterpoise method of Boys and Bernardi [43]. PECs were then modeled by fitting the ab initio points using the 6th degree Rydberg analytical function (termed here as Ryd6) [44], which is written as, where a j are adjustable coefficients, and R e represents the equilibrium position determined from quantum-chemical calculations. The parameters were optimized employing the Generalized Simulated Annealing (GSA) [45,46] method, required to minimize the cost function given by where E ab,i corresponds to the ab initio energies, E theo,i is the fitted energy and N p is the number of points. Rovibrational Spectroscopic Constants. After the complete description of the PECs of the interacting dimers, we solved the nuclear Schrödinger equation (NSE) to obtain its dynamic properties. For the fullerene dimers studied here, the radial NSE, under the Born-Oppenheimer approximation, is given by where F(R) = RΨ υ,J (R) and µ is the reduced mass of the system. The term U e f (R) in Equation (4) is the effective potential: which includes the intermolecular potential U(R) and centrifugal distortions due to the coupling with the rotational levels. By solving the NSE given by Equation (4), we are considering the dimers as two "large atoms", i.e., we are neglecting all intramolecular effects in an analogy to the methodology of Baggio and co-workers [17]. The spectroscopic constants of the fullerene dimers were obtained employing two different approaches. The first approach consists in solving the NSE to obtain its eigenvalues and then inserting them in a set of equations for the spectroscopic constants as described in refs. [20,21,47].
In this work, the NSE was solved using the Discrete Variable Representation (DVR) method [47]. The spectroscopic constants given in Equation (6) are the harmonic frequency ω e , the anharmonicity constants ω e x e and ω e y e α e and γ e are the rovibrational coupling constants.
The second approach is the method developed by Dunham [48]. In Dunham's method, the spectroscopic constants are calculated by taking derivatives of higher orders of the PEC expanded as a Taylor series around the equilibrium distance R e . The second (d 2 ), third (d 3 ), and fourth (d 4 ) order derivatives of U(R) equate to the following spectroscopic constants [15]. and Rovibrational population analyses as a function of the temperature were conducted using Boltzmann statistics, implementing a rovibrational partition function to account for the non-rigid rotator and anharmonic effects. More details about the explicit form of the partition function are provided in ref. [17]. The rovibrational level occupations were assessed in the range of 10-400 K.
Lifetime of the van der Waals complexes. Slater's theory [49] was employed to measure the lifetime of the dissociative process of the fullerene dimers. The unimolecular dissociation of the complex is supposed to occur when the interaction coordinate reaches the dissociation energy threshold, D e [50]. If the complex has enough energy to reach D e , the frequency of dissociation becomes simply the vibrational frequency, ω e , of the interacting molecules, and the dissociation rate constant becomes: where R = N A k b is the ideal gas constant, and ε 0,0 is the rovibrational ground state of the fullerene dimer. The reciprocal of Equation (10) renders the dissociation lifetime of the complexes.
According to Slater's theory, lifetime estimations are only valid for systems under high/intermediate pressures.
QTAIM Characterization of the NCIs. To investigate the type of intermolecular interaction, we take advantage of the topological analysis provided by Bader's theory, QTAIM [51][52][53]. In the QTAIM approach, every topological idiosyncrasy of ρ(r) (whether a maximum, minimum, or saddle point) is associated with a specific point in space, referred to as the critical point, where the gradient of ρ(r) vanishes [51]. Currently, the QTAIM approach has found large applicability in the study of many electronic properties; specifically, it has been employed in the study of weak interacting systems [10,11,54]. The topological characterization of the NCIs was studied in the case of fullerene dimers. The NCIs were classified according to the topological properties at the Bond Critical Point (BCP): the electron density (ρ(r)), Laplacian of electron density (∇ 2 ρ(r)), and energy density (H CP )). All QTAIM properties were computed using the wavefunction analysis program Multiwfn.
Reduced density gradient characterization of the NCIs. In association with QTAIM, reduced density gradient (RDG) [30,32] was employed to characterize the NCIs. The RDG analysis provides a dimensionless index, presented in Equation (12), called reduced density gradient, s, which is based in the electron density ρ(r) and its derivate (∇ρ(r)) [32].
For values of ρ close to zero, s tends to diverge, typically in intra-or intermolecular regions. To characterize the type of interaction at these critical points, a real space function, sign(λ 2 )ρ, is used [32]. Isosurfaces for RDG can also be drawn, which allows visualizing the different types of non-covalent interactions. RDG analysis has been applied to help categorize non-covalent interactions in different systems [31,55]. RDG properties were computed using the wavefunction analysis free program Multiwfn [56]. The draw of the isosurfaces and molecules for both QTAIM and RDG analysis was made with VMD software version 1.9.3 [57]. To ensure that the NCIs would be well analyzed, the fullerenes dimers were fully optimized at the equilibrium distance obtained from the PEC calculations.

Conclusions
The present work was devoted to obtaining the spectroscopic properties and characterizing the intermolecular potential and hence, the stability of the fullerene dimers (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 by means of the solution of the Nuclear Schrodinger Equation. The fitting procedure of the PECs obtained proved to be satisfactory, resulting in errors smaller than 0.05 kcal mol −1 in the equilibrium region. The interaction between the buckyballs increases with the system size. Lifetime calculations showed that the studied dimers are stable for more than 1 ps, which increases accordingly with the number of carbon atoms. Both Dunham's and DVR's methods returned similar results, ensuring the reliability of the methodology. The spectroscopic constants ω e , ω e x e and B e were determined and revealed that the anharmonicity effects associated with the ω e x e constant become more operative for the smallest dimers representing one-tenth of the harmonic frequency for the (C 20 ) 2 and (C 24 ) 2 dimers. The vibrational population for each state and the temperature effects of the Boltzmann distribution of the quantized levels were calculated. At room temperature (300 K), less than 25% of the dimers lie on the vibrational ground state, which indicates that, for a good description of the normal modes, the inclusion of higher excited states is mandatory. On the other hand, at a temperature of~10 K, almost 100% of the dimers are in the ground state.
QTAIM and RDG analyses indicate that the fullerenes dimers are stabilized due to van der Waals interactions. Investigation of the bond path and BCP descriptors enabled the determination of the preferential interaction configuration of the (C 60 ) 2 dimer and compared it with experimental results obtained with C 60 in-bulk. We expect that the theoretical data provided here should be useful for future studies and dimers applications in different areas.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules28135023/s1, Figure S1: Potential energy curves (top) fitted with analytic Rydberg function for the (C 20 ) 2 dimer calculated at ωB97xD/6−31G(d) level of theory. Point-by-point error (bottom) is also plotted to magnify the overall fitting accuracy in the entire range of internuclear separation; Figure S2: Potential energy curves (top) fitted with analytic Rydberg function for the (C 24 ) 2 dimer calculated at ωB97xD/6−31G(d) level of theory. Point-by-point error (bottom) is also plotted to magnify the overall fitting accuracy in the entire range of internuclear separation; Figure S3: Potential energy curves (top) fitted with analytic Rydberg function for the (C 36 ) 2 dimer calculated at ωB97xD/6−31G(d) level of theory. Point-by-point error (bottom) is also plotted to magnify the overall fitting accuracy in the entire range of internuclear separation; Figure S4: Potential energy curves (top) fitted with analytic Rydberg function for the (C 60 ) 2 dimer calculated at ωB97xD/6−31G(d) level of theory. Point-by-point error (bottom) is also plotted to magnify the overall fitting accuracy in the entire range of internuclear separation; Figure S5: Potential energy curves (top) fitted with analytic Rydberg function for the (C 70 ) 2 dimer calculated at ωB97xD/6−31G(d) level of theory. Point-by-point error (bottom) is also plotted to magnify the overall fitting accuracy in the entire range of internuclear separation; Figure S6: Potential energy curves (top) fitted with analytic Rydberg function for the (C 84 ) 2 dimer calculated at ωB97xD/6−31G(d) level of theory. Point-by-point error (bottom) is also plotted to magnify the overall fitting accuracy in the entire range of internuclear separation; Figure S7: (A) Optimized C 60 -fullerene dimer as determined from ωB97xD/6−31G(d) calculation depicting the most stable orientation, where the pentagon of C 60 faces the inter-pentagon bond of the adjacent monomer. (B) side view and (C) superimposed view to highlight the relative orientation between the monomers; Table S1: Parameters utilized for the potential fit through the Rydberg function for the (C 20 ) 2 , (C 24 ) 2 , and (C 36 ) 2 dimers.; Table S2: Parameters utilized for the potential fit through the Rydberg function for the (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 dimers.; Table S3: Quantized energies of the first ten vibrational levels (in cm −1 ) for the (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C70) 2 , and (C 84 ) 2 dimers; Table S4: Lifetime values for the (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C 70 ) 2 and (C 84 ) 2 dimers from a temperature of 200 K until 500 K; Table S5: Coefficients of fitted expressions of the ωB97XD/6−31G(d)-derived intermolecular potential for the (C 60 ) 2 non-covalent dimer. Fitting parameters for the Ryd6 analytical potential are shown in Table S2. Table S6: Electron density ρcp (in e/a 0 3 ), the Laplacian of the electron density ∇ 2 pcp (in e/a 0 5 ), and the total energy HBCP (in Hartree) associated with the bond-critical point occurring along the bond path connecting each monomer for the (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 , (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 dimers, separated by a distance R e . Table S7: Coordinates for the C 20 , C 24 , C 36 , C 60 , C 70 , and C 84 optimized fullerenes; Table S8: Coordinates for the (C 20 ) 2 , (C 24 ) 2 , (C 36 ) 2 (C 60 ) 2 , (C 70 ) 2 , and (C 84 ) 2 dimers at the equilibrium point. References [58][59][60] are cited in the supplementary materials.