Simulating the structural diversity of carbon clusters across the planar to fullerene transition

TL;DR

This study uses advanced simulations to explore the structural diversity of carbon clusters during the transition from planar forms to fullerenes, revealing how energy and hybridization influence their configurations.

Contribution

It introduces a comprehensive simulation approach combining REBO potentials and replica-exchange MD to analyze isomers of C$_n$ clusters, highlighting the role of sp$^2$ hybridization in stability.

Findings

Structural diversity increases with energy and cluster size.

Low-coordinated chains and branched structures become more prevalent at higher energies.

Stability is mainly driven by sp$^2$ hybridization rather than geometry.

Abstract

Together with the second generation REBO reactive potential, replica-exchange molecular dynamics simulations coupled with systematic quenching were used to generate a broad set of isomers for neutral C$_n$ clusters with $n=24$, 42, and 60. All the minima were sorted in energy and analyzed using order parameters to monitor the evolution of their structural and chemical properties. The structural diversity measured by the fluctuations in these various indicators is found to increase significantly with energy, the number of carbon rings, especially 6-membered, exhibiting a monotonic decrease in favor of low-coordinated chains and branched structures. A systematic statistical analysis between the various parameters indicates that energetic stability is mainly driven by the amount of sp$^2$ hybridization, more than any geometrical parameter. The astrophysical relevance of these results is discussed in the light of the recent detection of C$_{60}$ and C$_{60}^+$ fullerenes in the interstellar medium.

Tables

Table 1. Table 1: Point group and absolute binding energy of the reference structures of C 24 , C 42 and C 60 obtained with the REBO potential.

Reference structure	Point group	Binding energy (eV/atom)
C24	$D_{6h}$	6.237
C42	$D_3$	6.614
C60	$I_h$	6.842

Table 2. Table 2: Parameters used for determination of temperatures in the REMD simulations.

Cluster	$T_{\min }$ (K)	$T_{\max }$ (K)	$M$	$R_s$ (Å)
C24	1500	6500	12	4.0707
C42	2500	6500	14	4.9030
C60	2500	6500	16	5.5209

Table 3. Table 3: Linear ( R 𝑅 R ) and partial correlation coefficients (PCC) between various parameters and isomer energy obtained from the databases of minima for C 24 , C 42 and C 60 .

	C24		C42		C60
	$R$	PCC	$R$	PCC	$R$	PCC
$R_g^2$	0.436	-0.217	0.517	0.078	0.623	0.183
$A_3$	0.365	0.056	0.444	0.112	0.433	0.103
$S$	0.313	0.085	0.272	0.066	0.206	0.031
sp2	-0.285	-0.114	-0.906	-0.852	-0.965	-0.933

Equations

$$\text{mbox}acc({\mathbf{R}}_i\rightleftharpoons {\mathbf{R}}_j)=\text{min}\left\{1,\exp \left[({\beta }_i-{\beta }_j)(E({\mathbf{R}}_i)-E({\mathbf{R}}_j))\right]\right\}$$

$$Q^{\alpha \beta }=\frac{1}{n}\sum _{i=1}^n(r_i^{\alpha }-{\bar{r}}^{\alpha })(r_j^{\beta }-{\bar{r}}^{\beta }),\alpha ,\beta =1,2,3,$$

$$R_g^2=\frac{1}{n}\sum _{i=1}^n({\mathbf{r}}_i-\bar{\mathbf{r}}{)}^2=\text{Tr}\mathbf{Q},$$

$$A_3=\frac{3}{2}\frac{\text{Tr}{\mathbf{D}}^2}{(\text{Tr}\mathbf{Q}{)}^2},$$

$$S=27\frac{\text{det}\mathbf{D}}{(\text{Tr}\mathbf{Q}{)}^3},$$

$$\rho (r>0)=\frac{1}{4\pi r^2}\sum _{i=1}^n\delta (r-\left|{\mathbf{r}}_i-\bar{\mathbf{r}}\right|),$$

$$\text{sp}\text{: }\{\begin{array}{l}{N}_i=1\text{ or }2,\\ {\theta }_k>170^{\circ }\forall k.\end{array}$$

$${\text{sp}}^2\text{: }\{\begin{array}{l}{N}_i=2\text{ or }3,\\ 95^{\circ }<{\theta }_k<135^{\circ }\forall k,\\ \text{Var}({\theta }_k)<20^{\circ },\end{array}$$

$${\text{sp}}^3\text{: }\{\begin{array}{l}{N}_i=4,\\ 85^{\circ }<{\theta }_k<125^{\circ }\forall k,\\ \text{Var}({\theta }_k)<20^{\circ }.\end{array}$$

$$g(r)=\frac{1}{4\pi r^2{\rho }_S}\sum _{i<j}\delta (r-\left|{\mathbf{r}}_i-{\mathbf{r}}_j\right|),$$

$$\text{mbox}PCC(\{X_i\},\{X_j\})=-\frac{P_{ij}}{\sqrt{P_{ii}{P}_{jj}}},$$

$$\tilde{E}$$

Full text

I Introduction

From a few atoms to bulk matter, carbon clusters show a significant ability to hybridize in sp, sp2 or sp3 chemical bonds, reflecting at finite size the wide allotropy of bulk carbon matter. Depending on experimental conditions, carbon clusters can be produced into a very large variety of isomers that have been probed by many groups for more than two decades von Helden et al. (1991, 1993); López et al. (1997); von Orden and Saykally (1998); Handschuh et al. (1995); Lifshitz (2000); Kent et al. (2000); Ueno and Saito (2008); Kosimov et al. (2010); Straatsma et al. (2017); Kono et al. (2018). Below about 20 carbon atoms, (1D) chains and (2D) rings have been identified as the most stable isomers von Helden et al. (1991); von Orden and Saykally (1998); Brito et al. (2018) while (3D) fullerenes were shown to be the most stable form of larger carbon clusters Tománek and Schluter (1991); Kent et al. (2000).

Recently there was a surge of interest in the relaxation dynamics in carbon-based clusters following their brief excitation typically produced by energetic particle collisions Delaunay et al. (2015); Simon et al. (2018) or short laser pulses Marciniak et al. (2015); Ji et al. (2017); Zhen et al. (2016). Such experiments are largely motivated by the increasing evidence from astronomical observations that pure and hydrogenated carbon clusters are indeed present in the interstellar medium (ISM). While only the smallest molecules up to C5 were conclusively detected until 2001 Bernath et al. (1989); Maier et al. (2001), the much more organized fullerenes C60 and probably C70 were observed in the ISM owing to their very characteristic vibrational infrared emission bands Sellgren et al. (2010); Cami et al. (2010). The possible detection of the cation C${}_{60}^{+}$ from its infrared emission bands was also suggested Berné et al. (2013). These fullerene bands accompany the so-called aromatic infrared bands (AIBs), which trace polycyclic aromatic aliphatic mixed hydrocarbons widely observed in the ISM Léger and Puget (1984); Allamandola et al. (1985); Sloan et al. (2007); Pino et al. (2008); Acke et al. (2010); Sloan et al. (2014); Pilleri et al. (2015). Interestingly, infrared spectroscopy was also used to characterize smaller carbon clusters in laboratory experiments by Straatsma and coworkers Straatsma et al. (2017).

The high level of chemical organization of fullerenes necessarily raises questions regarding their formation under the harsh conditions of astrophysical environments. So far, essentially two scenarios have been proposed to explain the presence of such molecules in the ISM. The so-called bottom-up hypothesis, in which fullerenes would be formed by coalescence of smaller entities Yamaguchi and Maruyama (1998); Irle et al. (2003); Zheng et al. (2005), is rather unlikely in the ISM owing to its extremely low densities of matter Berné et al. (2015). In the top-down scenario, fullerenes originate from the decay of larger compounds subject to energetic excitation (cosmic rays, shocks, VUV irradiation…) followed by stepwise isomerization Berné et al. (2015); Chuvilin et al. (2010); Montillaud et al. (2013); Pietrucci and Andreoni (2014); Zhang et al. (2013). In this respect, amorphous carbon clusters were suggested to play a possibly important role Powles et al. (2009); Jones et al. (2017); Sinitsa et al. (2017). Such formation pathways thus question the existence of products that are concomitantly formed and may be stable in space, spanning the range between amorphous and fullerenic carbon clusters.

Experimental observation of the possible formation mechanisms will necessarily rely on spectroscopy, and theoretical interpretation of the measured spectra requires appropriate candidate structures. Unfortunately, the typical approach until now has usually been biased towards certain chemical types such as polyaromatic compounds Bauschlicher et al. (2018), only limited effort being devoted to arbitrary or amorphous conformations despite their known astrophysical relevance Duley and Hu (2012); Jones et al. (2013, 2017); Boutéraon et al. (2019). According to the theoretical study by Kent and coworkers Kent et al. (2000), 24 carbon atoms are needed to produce the first polyaromatic flake, and stable fullerenes arise above the size 26, although not necessarily the most stable conformation. Size 60 is notorious to support fullerenes as their most stable structure. In the absence of known magic numbers in this size range, and as common in clusters physics, 2D and 3D structures are thus likely to coexist von Helden et al. (1991).

The present article aims at exploring the structural diversity of carbon clusters as a function of their most important feature, namely their internal energy, in the size range where they undergo the flake-to-fullerene transition. Although the degree of chemical ordering is obviously expected to vary with increasing excitation energy, the extent of chemical and conformational variety remains undocumented so far, and as far as we are aware our fundamental study represents the first attempt to construct a library of carbon cluster structures in an unbiased way, not focusing on the lowest-energy structures only but giving as much attention to higher-energy isomers that could be formed on the interstellar fullerenic road .

To reach this goal we have systematically explored the conformational landscapes of selected carbon clusters Cn across the planar to fullerene transition and containing $n=24$, 42, or 60 atoms, using advanced molecular dynamics methods and systematic quenching. The numerous isomers thus obtained were sorted and analyzed using a range of structural order parameters, some of them to quantify the nature of the chemical bonds within the cluster. Because this computational investigation is highly statistical, we relied on a simplified but realistic description of the potential energy surface based on the second generation reactive REBO bond-order potential Brenner et al. (2002). Such an approach has already been used in the past to study the structure of carbon clusters in the size range up to 55 atoms by Kosimov and coworkers Kosimov et al. (2010), who predicted a transition from ring structures to graphene flakes occuring above 18 atoms, without reporting any fullerene. In contrast with earlier computational investigations, our approach here is highly statistical and does not focus on the lowest-energy structures. The results thus obtained shed light on the local chemical ordering and global structural arrangement of carbon clusters as their energy varies down to low-energy polyaromatic isomers.

The article is organized as follows. In the next section, we describe the different computational tools used for generating and characterizing the structural diversity of selected carbon clusters at finite internal energy. The results are presented and discussed in Section III in relevance to astrophysical implications. In Section IV we present a systematic statistical analysis in which we correlate the relative energetic stability of the various conformers to their structural and chemical features. Some concluding remarks are finally given in Section V together with ongoing or future extensions.

II Methods

II.1 Sampling the conformational landscape

The energy landscape of carbon clusters Cn is characterized by an exponentially increasing number of local minima and transition states with increasing size $n$, and cannot be sampled exhaustively as soon as this number exceeds a few tens. Furthermore, the barriers separating the various local minima or even more distant funnels on the landscape are likely to be high and make traditional simulation methods based on molecular dynamics or Monte Carlo methods poorly efficient. More importantly, we are not interested in addressing the global optimization problem specifically by focusing on the lowest-energy structures, but rather to explore conformations that are physically relevant also at high energies, as could be produced e.g. upon photonic or collisional excitation, though still under isolated conditions.

We used replica-exchange molecular dynamics (REMD) simulations Swendsen and Wang (1986); Marinari and Parisi (1992) to circumvent this broken ergodicity issue and achieve a broader sampling of the potential energy surface over extended energy ranges.

REMD simulations were performed using the LAMMPS program Plimpton (1995) by propagating $M$ trajectories at fixed temperatures $T_{i}$ ($i=1\dots M$) and occasionally attempting exchanges between the two configurations ${\bf R}_{i}$ and ${\bf R}_{j}$ of neighboring replicas $i$ and $j=i\pm 1$. Such an exchange is accepted using the following acceptance rule Sugita and Okamoto (1999)

[TABLE]

where we have denoted $\beta_{i}=1/(k_{B}T_{i})$ and $E({\bf R})$ the potential energy at configuration ${\bf R}$, discussed below.

The efficiency of this exchange process depends on the overlap between the thermal distributions at temperatures $T_{i}$ and $T_{j}$, which itself is driven by various factors, primarily the two temperatures but also the size $n$ which affects the width of the individual distributions. It is thus important to carefully choose the set of temperatures $\{T_{i}\}$ for each cluster size, taking also into account the need for the upper temperature $T_{M}$ to be high enough to ensure an efficient exploration of configurational space, still below the vaporization temperature. Since the thermal distributions both shift and broaden with increasing temperature, the difference between successive temperatures $T_{i+1}-T_{i}$ must also increase with $i$. Here we employed a geometrical progression in $T_{i}$, namely $T_{i+1}=\alpha T_{i}$, which is optimal in the harmonic limit Mitsutake et al. (2001), to which we have added manually some temperatures to increase the exchange probability.

The MD trajectories were integrated with a $0.1$ fs time step. A Nosé-Hoover thermostat was used to keep temperatures constant with a damping constant of $10$ fs. Exchange between nearest replicas were attempted once every $2\times 10^{4}$ MD steps and the simulation was propagated for 100 ns. Configurations were periodically saved for further structural analysis and quenching every 4 ps, resulting in a total of 25 000 structures per replica.

During the simulations it was also important to forbid fragmentation as it is an irreversible process preventing the correct sampling of size-selected compounds. Here we used a simple spherical harmonic potential container with a 500 eV$/$Å2 spring constant and a radius $R_{s}$ whose value was chosen differently for the three cluster sizes $n$ to reach a common density of 1.7 g/cm3, relevant for disordered polyaromatic materials such as soot Yon et al. (2015). The initial conditions of the REMD simulations were taken as the global energy minima for each cluster, also denoted as reference structures.

II.2 Potential energy surface for carbon clusters

The systematic production of large samples of minima at finite temperature, in which chemical bonds are easily broken and reformed, requires an efficient but chemically realistic method for the potential energy surface $E({\bf R})$. Here methods with an explicit treatment of electronic structure or based on first principles are not practical, especially considering the possible multireference character of small carbon clusters that would make the solution to the electronic problem already cumbersome Montagnon and Spiegelman (2007).

Only an explicit potential energy surface is currently able to handle the tremendous number of configurations gathered with the currently available computational resources, approximate schemes based on tight-binding Wang and Mak (1995); Van-Oanh et al. (2002) or density-based tight-binding Rapacioli et al. (2015) remaining still too expensive for the present rather large clusters. A few realistic potentials are available for carbon, which correctly account for bond breaking and formation and the various hybridization environments displayed by carbon Tersoff (1988); Brenner (1990); Brenner et al. (2002); Stuart et al. (2000); Marks (2000); Ghiringhelli et al. (2005); Zhou et al. (2015). Here we have chosen the adaptive the second-generation reactive bond order (REBO) potential of Brenner Brenner et al. (2002). Brenner type potentials have been notably used to study energetic and mechanical properties of nanotubes Ni et al. (2002), the formation process of fullerenes Yamaguchi and Maruyama (1998); Makino et al. (1997) but also to describe reactivity and formation of pure carbon clusters or hydrocarbons in the astrophysical context  Patra et al. (2014); Delaunay et al. (2015), making it a natural choice also for the present investigation. Moreover it has been shown that the REBO potential gives better results to describe small carbon clusters compared with other bond-order potentials Lai et al. (2016).

II.3 Reference structures

Our REMD exploration requires some initial structures. We have thus conducted a distinct search for the lowest-energy configurations of C24 by comparing various remarkable structures such as polycyclic, chains and rings conformers. In agreement with the earlier study by Kent et al. Kent et al. (2000), the lowest-energy structure found for C24 is the planar, fully dehydrogenated coronene with $D_{6h}$ point group. For C60 the natural reference structure is buckminsterfullerene with point group $I_{h}$, besides a set of 1811 alternative but higher energy fullerenes Sure et al. (2017). Finally, C42 was chosen as an intermediate size cluster with a propensity to form fullerenes. For this cluster, 45 nonequivalent fullerene isomers could be identified Fowler and Manolopoulos (1995) and the REBO potential predicts that the most stable isomer is the only fullerene with $D_{3}$ symmetry, in accordance with DFT calculations Sun et al. (2005); Malolepsza et al. (2009) but at variance with the earlier study by Kosimov et al. Kosimov et al. (2010) who found a graphene flake as for C24.

These reference structures are depicted in Figs. 3-C24-a),  3-C42-a), and 3-C60-a) and their absolute binding energies obtained with REBO are given in Table 1.

For buckminsterfullerene, the binding energy is reasonably close to the experimentally known value of 6.7 eV/atom Dresselhaus et al. (1996).

II.4 Computational details

The number of temperatures allocated for the REMD trajectories, their lowest and highest values and the container radius employed in the simulations are given for each system size in Table 2.

For all cluster sizes and for each of the $M$ replicas, 25 000 configurations were generated. Given the numbers of replicas employed for each system, a total of 300 000 instantaneous configurations were kept for further analysis for C24, 350 000 instantaneous configurations for C42 and 400 000 configurations for C60.

II.5 Identification of stable structures

From the large sets of instantaneous configurations gathered at finite temperature, the local minima were obtained by systematic quenching using here the Hessian-free truncated Newton method as implemented in LAMMPS Plimpton (1995) and disregarding the hard-wall spherical container. Only connected structures were subsequently kept for further analysis, disconnected configurations being discarded. Here fragments are identified using a maximum nearest neighbor distance of 1.85 Å. In order to speed up the analysis, two locally minimized structures were further considered to be identical if their energy difference lies below 0.01 meV, in which case the highest was discarded as well. After this optimization and screening stage, the numbers of distinct configurations saved for further processing was equal to 51 901, 240 305, and 236 394 for C24, C42, and C60, respectively. We interpret this slightly smaller number for the larger cluster size as reflecting its more magic character, the diversity of conformers at a same upper temperature (here 6500 K) being lower than in a non-magic system such as C42 Wales (2003). However, the energy distributions shown below are essentially robust when considering only the first or second half of the configurations harvested during the REMD trajectories, suggesting that our samples are statistically representative. Yet, the smaller number of distinct configurations for C60 would probably require further scrutiny, and it could be useful to apply alternative approaches employing systematic local minimization such as basin-sampling Bogdan et al. (2006) in future extensions of this work dealing with even more complex systems.

II.6 Structural analysis

Both local and global parameters were used to characterize the various isomers obtained for the three clusters sizes. Global parameters provide information about the overall shape and atomic distribution around the center of mass, while local parameters give insight into the chemical arrangement at the atomic level.

For a $n$-atom cluster Cn with equilibrium configuration ${\bf R}=\{{\bf r}_{i}\}$ for $i=1\dots n$, the $3\times 3$ gyration tensor ${\mathbf{Q}}$ is defined from its components $Q^{\alpha\beta}$ as Šolc and Stockmayer (1971); Theodorou and Suter (1985); Blavatska and Janke (2010); Calvo et al. (2012)

[TABLE]

where $r_{i}^{\alpha}$ are the Cartesian coordinates ($x$, $y$ or $z$) of atom $i$, $\bar{r}^{\alpha}$ are the corresponding coordinates for the cluster center of mass Blavatska and Janke (2010). Three rotationally invariant quantities can be defined from the tensor ${\bf Q}$ that respectively measure the geometrical extension, the asphericity, and the prolateness of the atomic distribution. The squared radius of gyration $R_{g}^{2}$ is first defined by

[TABLE]

where $\text{Tr}\,\mathbf{Q}$ stands for the trace of the gyration tensor. The other two quantities are defined from the traceless tensor $\mathbf{D}=\mathbf{Q}-\frac{1}{3}\mathbf{I}\text{Tr}\mathbf{Q}$, where $\mathbf{I}$ is the 3$\times$3 identity matrix. The asphericity parameter $A_{3}$ is defined by Calvo et al. (2012)

[TABLE]

The asphericity varies from 0 for a purely spherical system to 1 for a perfectly linear structure. Finally, the prolateness $S$ is given for configuration $\mathbf{R}$ by

[TABLE]

and varies from $S=-1/4$ for a perfectly planar disk (oblate structure) to $S=2$ for a perfectly linear chain (prolate structure).

The tendency to form hollow, planar or close-packed structures was investigated also from the radial density $\rho(r)$ around the center of mass,

[TABLE]

where $\delta$ refers to the Dirac function.

Another useful information for pure carbon clusters having a tendency for sp and sp2 hybridization is the number of rings of each size. Here we counted the numbers of 3-, 5-, 6-, and 7-membered rings from the analysis of nearest neighbor atoms connectivity and denote them as Rℓ where $\ell$ is the length of the ring.

Turning now to local order parameters, the hybridization state of each atom was quantified for all stable configurations ${\bf R}$ of the databases. In many computational studies that do not explicitly compute the electronic structure Galli et al. (1989, 1990); Wang et al. (1993); Bhattarai and Drabold (2017); Dozhdikov et al. (2017) hybridization is defined based on coordination numbers only. However, such a definition cannot account for the chemical complexity and possibly reactive atoms that are under coordinated. Here we use both coordination and geometric information to assign hybridization states. More precisely, we define $N_{i}$ the coordination number of atom $i$, neighbors lying within 1.85 Å from atom $i$, and we evaluate all angles in which atom $i$ is a vertex. There are $M_{i}=N_{i}(N_{i}-1)/2$ such angles which we denote by $\theta_{k}$ for $k=1\dots M_{i}$, leaving the dependence on atom $i$ as implicit. Atom $i$ is then said to be in sp hybridization state provided that it is not overcoordinated and the angles in which it is central are close enough to 180∘:

[TABLE]

Likewise, the atom is assigned sp2 hybridization with the following criteria

[TABLE]

where $\text{Var}(\theta_{k})$ denotes the variance of all angles for which $i$ is a vertex. Here again an upper limit for the angular boundary of 20∘ is chosen in such a way as to include atoms involved in hexagonal and pentagonal rings as properly sp2. The variance limit ensures that for $N_{i}=3$ all 4 atoms remain close to a common plane. It was notably chosen to include fullerene structures (for buckminsterfullerene Var$(\theta_{k})=6^{\circ}$). Finally atom $i$ is considered to be in a sp3 hybridization state accordingly with

[TABLE]

Here the variance condition on the angles is required to ensure the atom has a 3D tetrahedral environment. With the above definitions, it may occur that a given atom is neither sp, sp2 or sp3, in which case it will be referred to as ambiguous. This notably occurs for ringlike structures containing too few atoms, hence producing angles that exceed 10∘, one such example being illustrated in Fig. 3-C24-d). With the present criterion on angles, at least 36 atoms would be needed in the perfect ring for the sp2 assignment to be recovered.

Finally, as our last local structural quantity we have determined the pair correlation function $g(r)$ from

[TABLE]

where $\rho_{S}=3n/4\pi R_{S}^{3}$ is the original density used to constrain the REMD simulations. The pair correlation function ignores the overall distribution of atoms within the container and thus provides a complementary information to the radial distribution. Furthermore it is indirectly sensitive to the chemical nature of the bonds, the various hybridization states being related to different equilibrium distances in carbon-carbon bonds. For example, with the REBO potential the CC distance of acetylene (C2H2) is 1.21 Å, increasing to 1.32 Å in ethylene (C2H4) and to 1.54 Å in ethane (C2H6). These values are in very satisfactory agreement with quantum chemical calculations at the DFT/M06-2X/6-311++G** level, which for the same molecules yield 1.1981, 1.3260, and 1.5254 Å, respectively.

III Results and discussions

The distributions of local minimum energies obtained for the three clusters C24, C42 and C60 with our computational scheme and the REBO potential are represented in Fig. 1. Here the energies were shifted relative to the reference structure energy to highlight the ranges explored by the REMD simulations, these isomers being depicted at the top of Fig. 3 and the ranges themselves highlighted by arrows.

For C24, the distribution exhibit a single broad ($\sim$0.2 eV/atom) and asymmetric peak located around 0.31 eV/atom. For C42 the distribution exhibits two peaks at 0.11 and 0.56 eV/atom, showing structures on a large energy scale up to $\sim$0.85 eV/atom. Similarly the distribution for C60 exhibits two peaks located at 0.11 and 0.77 eV/atom and the overall energies of the structures vary up to $\sim$1 eV/atom. To understand the difference between C24 and the two larger cluster sizes and the nature of the two peaks in C42 and C60, the distributions of quenched structures were determined as a function of the squared radius of gyration and the proportion of sp2 carbon atoms (Fig. 2).

As for the energy distribution, Fig. 2 shows a single peak for C24 and two peaks for C42 and C60. The lowest energy peaks in C42 and C60 can be identified to the peak with larger proportion of sp2 and smallest squared gyration radius. The highest energy peaks in C42 and C60 can be identified to the peak with lowest proportion of sp2 and largest squared gyration radius. A direct examination of the structures allows the lowest-energy peak to be assigned to cage-like structures while the highest-energy peak corresponds to pretzel-like conformations, as already identified by Kim and Tománek Kim and Tománek (1994) who simulated the melting of fullerenes using a tight-binding approach Kim and Tománek (1994).

In Fig. 3, typical structures are shown for different internal energies. At low energy ($\leq$0.1 eV/atom), C24 tends to be planar as in the reference structure. In this energy range, some deformed caged-structures and various fullerene isomers are mainly found for C42 and C60. A large proportion of atoms are found as sp2 with only a small number of atoms with ambiguous hybridization. As the internal energy increases, a clear increase in the number of ambiguous hybridization also occurs, concomitant with the appearance of pretzel-like structures around $\sim$0.2 eV/atom for C24 and $\sim$0.4 eV/atom for both C42 and C60. At higher energies, i.e. $\sim$0.5 eV/atom for C24 and $\sim$0.7 eV/atoms for C42 and C60, a trend to form branched atomic chains with a larger number of sp carbons is seen as well.

After this first qualitative analysis, we now proceed on a quantitative statistical footing based on the different parameters described in the previous section. The energetic stability of these isomers can be correlated with their structural properties, starting with the global structural order parameters provided by the gyration tensor. The energy-resolved distributions of parameters $R_{g}^{2}$, $A_{3}$ and $S$ are shown in Fig. 4 for the three cluster sizes. At a given energy the distribution is normalized.

At low energy, the three structural parameters smoothly converge to the corresponding values in the reference structure. Being shown on the same scale, the square gyration radii are rather similar for C24 and C42, despite different topologies (planar versus hollow cage). Such differences are better manifested on the asphericity, which vanishes for both cage isomers C42 and C60. From these 2D distributions the behavior of C24 stands out as different from the two larger clusters, with non monotonic variations in the average $R_{g}^{2}$ and $A_{3}$ with increasing isomer energy while the trends are all monotonically increasing in the latter case. A decrease in the square gyration radius near 0.2 eV/atom concomitant with a decrease in $A_{3}$ indicates spherical structures less extended than fully dehydrogenated coronene and corresponding to the 3D compact structures, as shown in Fig. 3-C24-d). Above this energy range, the distribution of squared gyration radius $R_{g}^{2}$ becomes much broader and reaches $\sim$20 Å2.

For the two larger clusters, the global structural indices display more regular variations with increasing isomer energy, metastable configurations exhibiting larger gyration radii, a greater asphericity and the prolateness remaining low in magnitude but with increasing fluctuations extending mostly to positive values. These fluctuations are most prominent above 0.5 eV/atom and again convey the greater structural diversity. In contrast, C24 remains planar until the energy reaches about 0.2 eV/atom, at which stage the rather sharp variations in the three structural indicators are consistent with the appearance of more compact configurations. In particular, the nearly spherical cage-like structures are found with very low $A_{3}$ in this energy range.

Above 0.4 eV/atom for C24, 0.6 eV/atom for C42 and 0.7 eV/atom for C60, the asphericity and prolateness indices both explore larger values much deviating from the reference structures. Visual inspection indicates that these high energy isomers are usually branched with several chains and a few rings only, as depicted in Fig. 3-C24-f), Fig. 3-C42-f), Fig. 3-C60-f). The linear chain isomers, for which $A_{3}=1$ and $S=2$, were indeed found for C24 and for C42 but not for C60 as they do not fit into the spherical container.

The radial densities sorted with increasing configuration energy are represented in Fig. 5 for the three cluster sizes.

The ranges of variations in the radial density match those exhibited by the global parameters originating from the gyration tensor. In particular, the formation of more compact structures in C24 below 0.2 eV/atom, the clear cagelike character in the two larger clusters at low energies manifested by a main peak in the radial density, and the loss of these fullerenic structures above 0.5 eV/atom for C42 and 0.6 eV/atom for C60 are all reflected in Fig. 5. Interestingly, carbon atoms also much more likely occupy the central regions when the energy exceeds 0.2–0.6 eV/atom depending on system size, which is consistent with the loss of hollow structures and the increasing occurrence of chains and branched configurations.

Turning to hybridization states, we first show in Fig. 6 and on the example of C60 only how the angles between connected triplets of carbon atoms are distributed when the central atom has two or three neighbors, which we anticipate to be potential candidates as sp and sp2 states, respectively.

The angles corresponding to regular rings R3, R5, R6 and R7 are also shown to highlight the occurrence of such regular polygons in the carbon clusters. At low energies, the angular distribution for carbon atoms having two neighbors only is peaked around 110–120∘, indicating a majority of pentagonal and hexagonal rings. The angles open at energies higher than 0.3 eV/atom, consistently with the formation or chains or larger rings. Five- and six-membered rings, which are predominant at low energies and in fullerene structures, concomitantly decrease above this same approximate energy threshold. Chains themselves show a signature at an angle of 180∘. In 3-coordinated atoms, the distribution is strongly peaked around 108∘ and 120∘ which are the expected values for buckminsterfullerene composed of perfect pentagonal and hexagonal rings. For both coordination numbers, three-membered rings arise above 0.8 eV/atom, but mostly as traces in 2-coordinated atoms and much more significantly in 3-coordinated atoms as seen through the increasing occurence of 60∘ angles. This difference suggests that in clusters that are compact enough, the three-membered rings lie in their inner regions. This is corroborated by examining their average distance to the center of mass, which in average approximately equates 83% of the gyration radius.

Based on the above analysis, the effects of excitation energy on the relative proportion of the various hybridation states can now be discussed. The fractions of sp, sp2, sp3, and ambiguous atoms were thus evaluated for all structures in our databases and for the three cluster sizes, the results being shown in Fig. 7.

For the three systems, the reference structures exhibit pure sp2 hybridization state, as expected for the present polyaromatic isomers. As energy increases, the extent of sp2 hybridization drops around energies where the most salient structural changes were noted earlier, that is approximately at 0.2 eV/atom for C24, 0.4 eV/atom for C42 and 0.6 eV/atom for C60. However, sp hybridization becomes significant only at energies higher than these thresholds, while no signature of sp3 is seen whatsoever. Ambiguous hybridization states thus populate the intermediate energy range where configurations become less compact, with a rather high amount of curved linear chains or large rings that do not fall in either of the sp or sp2 categories.

The steady increase in sp hybridization indicates the greater importance of linear chains in high energy structures. At intermediate energies, many configurations show fewer or shorter such chains, at the expense of rings or curved chains, whose atoms are interpreted as ambiguous until the angle becomes small enough and compatible with sp2 hybridization.

The same qualitative trends are noted in larger clusters, the rise in sp atoms and the maximum in the amount of ambiguous atoms being both shifted to higher energies. Extrapolating these trends, we speculate that in even larger clusters the propensity for sp2 hybridization would become even stronger and more robust against energy excitations, the proportion of linear chains being concomitantly lower and delayed to higher energies.

A complementary quantity is the number and size of the rings contained in the configurations as function of their energy. Figure 8 illustrates this specific property, in average, for the 3-, 5-, 6-, and 7-membered rings. 4-membered rings, which are much scarcer, were not considered in this study. Such rings are occasionally found but deviate significantly from the perfect square and exhibit angles closer to 80∘ (and a more standard angle near 120∘). In the case of C60 they can be seen in Fig. 6 for three-coordinated atoms as an horizontal spot with very low magnitude near 80∘.

For the three clusters, the reference structures only contain 5- and 6-membered rings. In all cases, the number of hexagonal rings steadily decreases with increasing energy. In both C24 and C60 this decrease benefits to pentagonal and heptagonal rings, consistently with the appearance of topological defects such as Stone-Wales pairs.

As the internal energy reaches the threshold values where global structural changes take place (0.2 eV/atom for C24, 0.5 eV/atom for C42, 0.7 eV/atom for C60), the numbers of rings having 5 or more atoms reaches a minimum and only residual 3-rings are found although with an increasing propensity. The loss of large rings is consistent with more atoms being present as linear chains in such structures, as depicted in Fig. 3. Three-membered rings being energetically rather disfavored, they only appear — occasionally — at the highest excitation energy as a way to connect linear chains into branched structures.

We have finally considered the pair distribution function as another structural indicator sensitive to the chemical arrangement and in particular hybridization state. Such quantities are shown in Fig. 9 for the three cluster sizes and at selected energies.

At low energy of 0.1 eV/atom and for all sizes, the pair distribution function is dominated by sp2 carbons. For C42 and C60, in the range 1.4–1.5 Å  carbon atoms are involved in hexagonal and heptagonal rings with the two corresponding peaks located at 1.43 and 1.47 Å, respectively. Note that for buckminsterfullerene the CC distance for a bond located between two hexagonal ring is 1.40 Å  while bonds connecting hexagonal and pentagonal rings together are closer to 1.45 Å. For the C42 reference structure with $D_{3}$ point group, the CC distance between two hexagonal rings is 1.41 Å, the CC distance between an hexagonal and a pentagonal ring varies between 1.45 and 1.46 Å. Finally the CC distance between two pentagonal rings varies between 1.48 and 1.49 Å. For C24, the CC distance of the $D_{6h}$ reference structure for the inner bonds varies between 1.44 Å and 1.45 Å. For the outer bond the distance varies between 1.37 Å and 1.41 Å.

As energy increases, ambiguous and sp hybridization states become increasingly important and are manifested by a narrower peak near 1.34 Å, and eventually 1.31 Å for singly coordinated atoms. A small residual peak above 1.6 Å is also found at very high energy. This peak results from the few carbon atoms involved in three-carbon rings.

The fullerene cage C60 and its cation were both recently observed in the ISM. Identification of the neutral buckminsterfullerene could be achieved owing to its few and very specific active bands Omont (2016). The observation of other clusters, even of the fullerene form, should be far more challenging. One important issue is to determine how such highly organized molecules could be formed in extremely dilute media, notably in the top-down assumption where C60 would originate from dehydrogenation and subsequent rearrangement of larger polycyclic aromatic hydrocarbons. Under such a model, our results indicate that the route to buckminsterfullerene is far from straightforward and encompasses multiple branched isomers along the way. It is noteworthy that the high energy configurations found in this work resemble the pretzel phase previously identified by Kim and Tománek Kim and Tománek (1994) in their simulation of fullerene melting. With the threshold energy reported here, the pretzel phase of C60 might be present as well in the ISM.

IV Statistical data analysis

In the previous section the relation between isomer energy and various structural parameters was highlighted, emphasizing the large configurational diversity arising as energy increases. Here we adopt a different but complementary point of view by questioning to which extent energetic stability is statistically related to these geometric parameters. More precisely, we have performed a systematic correlation analysis between isomer energy and the square gyration radius $R_{g}^{2}$, the asphericity $A_{3}$, the prolateness $S$ and the sp2 fraction, for the three databases of minima obtained for C24, C42, and C60.

The most straightforward way to quantify linear correlation between two sets of data is based on the traditional Pearson’s coefficient that we denote $R_{ij}$ for two sets of variables $X_{i}$ and $X_{j}$ among the five quantities of isomer energy, $R_{g}^{2}$, $A_{3}$, $S$ and sp2 fraction. We are mostly interested in correlations between energy (acting as $X_{i}$, our output parameter) and any of the four other quantities (acting as $X_{j}$ for $j\neq i$ and treated as input parameters).

One known issue with Pearson correlation coefficients is that they incorporate possibly strong correlations between the various input parameters, and an efficient way of removing these correlations consists of considering partial correlation coefficients (PCCs) Hamby (1994) instead. The PCC between variables $X_{i}$ and $X_{j}$ involves the full linear correlation matrix ${\bf R}=R_{ij}$ (including also those elements $R_{jk}$ between geometrical parameters $X_{j}$ and $X_{k}$ for $k\neq i$) and reads

[TABLE]

where $P_{ij}$ is the cofactor of the element $R_{ij}$ in the determinant $|{\bf R}|$ of matrix $\mathbf{R}$.

The linear ($R_{ij}$) and partial correlation coefficients obtained between the four geometrical parameters and the isomer energy are given in Table 3 for the three cluster sizes.

From these data, the smaller cluster C24 appears to behave somewhat differently from the two larger clusters C42 and C60. For these sizes that support fullerene-type isomers, linear correlations seen from $R$ quantities are highest with the sp2 fraction, the negative sign obtained for $R$ being consistent with the intuitive observation that isomers lowest in energy have the highest sp2 fraction. The three shape parameters always correlate positively but with Pearson coefficients never exceeding 0.65, such a poor linear correlation appears consistent with the scatter plots reported in Fig. 4.

Removing the correlations between the various parameters in the PCC quantities confirms this trend for the two larger clusters, and highlights the sp2 fraction as the most sensitive parameter causing relative energetic stability among conformers. However, for C24 this quantity performs not as well, and none of the quantities shows a strong partial correlation with energy. This can be understood because for C24, within the main peak of the energy distribution (between 0.2 and 0.6 eV/atom) the sp2 fraction remains relatively constant. The negative sign obtained with the square gyration radius is related with the markedly different distributions displayed in Fig. 4, where for C24 structures near 0.2 eV/atom are found to be more compact than the planar global minimum, at variance with the behavior noted for the two larger clusters.

The present statistical analysis thus suggests that, for fullerenes, the fraction of sp2 atoms is the main factor responsible for the energetic stability of the various conformers, geometric shape parameters appearing of lesser importance. Proceeding further, we have attempted to represent the energy $X_{0}$ as a simple function of the other parameters $\{X_{i}\}$, $1\leq i\leq 4$, using a simple linear functional form as

[TABLE]

where the 16 parameters $a$ are obtained by minimizing exactly the least square penalty function $|\widetilde{E}-X_{0}|^{2}$ summed over all minima.

The optimized fitted form for the isomer energy is shown in Fig. 10 against the reference value, for the three systems and as a 2D distribution.

These plots generally reflect the shapes of the distributions in Fig. 1, which are bimodal for the two larger clusters but only unimodal for C24. Even after optimizing the parameters, the representation of the isomer energy as a function of the three shape indicators $R_{g}^{2}$, $A_{3}$, $S$ together with the sp2 fraction appears rather approximate, especially for the fullerenes subsets. However, the trends are satisfactory with Pearson coefficients between the REBO and the fitted energies exceeding 0.95 for both C42 and C60. The smaller cluster C24 actually exhibits less chemical ordering and weaker correlations between energy and all other quantities analyzed, hence our representation of the energy as a linear form of these parameter unsurprisingly performs poorly.

V Concluding remarks

Owing to their additional free surface, carbon clusters are anticipated to exhibit a structural diversity that matches or even exceeds that of the known allotropy of bulk carbon. The present article was aimed at quantifying this diversity from a computational perspective, by performing a broad unbiased sampling of low-energy structures of carbon clusters containing a few tens of atoms. For this purpose the reactive semi-empirical REBO potential was employed for its ability to describe all hybridization states of carbon atoms together with REMD simulations as our main exploration tool. The simulations were further processed by systematically performing local optimizations in order to uncover the inherent structures and sort them as a function of energy. Various structural indicators of global or local character were evaluated to relate the structural and energetic features of these configurations to one another.

The lowest-energy structures obtained for the three sizes of interest here with 24, 42, or 60 atoms were found to be mostly polyaromatic and either planar as in the case of the fully dehydrogenated coronene flake for C24, or cage-like as in C42, buckminsterfullerene providing the most symmetric case of a spherical structure for C60. Excess energy appears as the main driving force causing the structural diversity of other isomers, and we evaluated 0.2 eV/atom as a threshold energy above which this diversity explodes, leading to a mixture of isomers that are neither fully planar nor perfectly cage-like but contain an increasing number of rings and connecting chains. Along this transformation the atoms initially keep a mostly sp2 character but evolve toward a greater proportion of sp hybridization (through an intermediate state where hybridization is ambiguous), no significant sp3 character being noted. It should be noticed that the spherical container used here to restrict to prevent excessively dissociated structures also disfavors elongated but connected chains as well as large planar structures for the bigger clusters, which could be entropically favored at high energies and thus be of significance as reaction intermediates. Periodic boundary conditions at fixed density or pressure or a Monte Carlo framework restricted to sample connected configurations only could complement the present approach and possibly generate new relevant configurations.

A systematic statistical analysis of the databases of minima gathered for the three systems was performed using linear and partial correlation coefficients in order to unravel a possible sensitivity of the relative isomer energy to the various indicators evaluated. For the two clusters supporting fullerenes, the sp2 fraction was found as the main driving force correlating with isomer stability, while none of the shape parameters was found particularly relevant. A linear representation of the energy as a function of the remaining parameters was proposed as an approximate description.

In view of astrophysical implications, one extension of the present work could consider the effects of a single charge on the structure of such carbon clusters. Ionized clusters are much more convenient to study in laboratory experiments owing to mass spectrometry selection, and their structures can be indirectly measured by techniques such as ion mobility von Helden et al. (1991, 1993). Cationic clusters could be modeled either by incorporating appropriate electrostatic corrections to the present REBO model, e.g. through fluctuating charges Rappé and Goddard (1991), or by approaches explicitly accounting for electronic structure such as tight-binding Van-Oanh et al. (2002); Kim and Tománek (1994) or DFTB Elstner and Seifert (2014).

Besides their ionic character, the pure carbon clusters studied here are an oversimplified description of the chemical ISM. The presence of hydrogen should be considered, especially in relation with the presence of polycyclic aromatic hydrocarbons as the starting point toward fullerene formation. To cope with the additional chemical complexity resulting from the presence of unlike atoms, improvements of the present REMD method or alternative approaches such as basin-sampling could be beneficial. Such investigations about the possible contribution of hydrogen at selected fractions will be the addressed in the near future.

VI Acknowledgments

The authors gratefully acknowledge financial support by the Agence Nationale de la Recherche (ANR) grant ANR-16-CE29-0025. This work was supported by grants from Région Ile-de-France and by the GDR EMIE 3533.