Simulating the structural diversity of carbon clusters across the planar to fullerene transition
Maëlle A. Bonnin, Cyril Falvo, Florent Calvo, Thomas Pino, Pascal Parneix
SSimulating the structural diversity of carbon clusters across theplanar to fullerene transition
Ma¨elle A. Bonnin, Cyril Falvo,
1, 2
Florent Calvo, Thomas Pino, and Pascal Parneix ∗ Institut des Sciences Mol´eculaires d’Orsay, Univ. Paris-Sud,Universit´e Paris-Saclay, 91405 Orsay, France. Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France
Abstract
Together with the second generation REBO reactive potential, replica-exchange molecular dy-namics simulations coupled with systematic quenching were used to generate a broad set of isomersfor neutral C n clusters with n = 24, 42, and 60. All the minima were sorted in energy and analyzedusing order parameters to monitor the evolution of their structural and chemical properties. Thestructural diversity measured by the fluctuations in these various indicators is found to increasesignificantly with energy, the number of carbon rings, especially 6-membered, exhibiting a mono-tonic decrease in favor of low-coordinated chains and branched structures. A systematic statisticalanalysis between the various parameters indicates that energetic stability is mainly driven by theamount of sp hybridization, more than any geometrical parameter. The astrophysical relevanceof these results is discussed in the light of the recent detection of C and C +60 fullerenes in theinterstellar medium. ∗ [email protected] a r X i v : . [ phy s i c s . a t m - c l u s ] M a r . INTRODUCTION From a few atoms to bulk matter, carbon clusters show a significant ability to hybridizein sp, sp or sp chemical bonds, reflecting at finite size the wide allotropy of bulk carbonmatter. Depending on experimental conditions, carbon clusters can be produced into a verylarge variety of isomers that have been probed by many groups for more than two decades [1–11]. Below about 20 carbon atoms, (1D) chains and (2D) rings have been identified as themost stable isomers [1, 4, 12] while (3D) fullerenes were shown to be the most stable formof larger carbon clusters [7, 13].Recently there was a surge of interest in the relaxation dynamics in carbon-based clustersfollowing their brief excitation typically produced by energetic particle collisions [14, 15]or short laser pulses [16–18]. Such experiments are largely motivated by the increasingevidence from astronomical observations that pure and hydrogenated carbon clusters areindeed present in the interstellar medium (ISM). While only the smallest molecules up toC were conclusively detected until 2001 [19, 20], the much more organized fullerenes C and probably C were observed in the ISM owing to their very characteristic vibrationalinfrared emission bands [21, 22]. The possible detection of the cation C +60 from its infraredemission bands was also suggested [23]. These fullerene bands accompany the so-called aro-matic infrared bands (AIBs), which trace polycyclic aromatic aliphatic mixed hydrocarbonswidely observed in the ISM [24–30]. Interestingly, infrared spectroscopy was also used tocharacterize smaller carbon clusters in laboratory experiments by Straatsma and coworkers[10].The high level of chemical organization of fullerenes necessarily raises questions regardingtheir formation under the harsh conditions of astrophysical environments. So far, essentiallytwo 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 ofsmaller entities [31–33], is rather unlikely in the ISM owing to its extremely low densitiesof matter [34]. In the top-down scenario, fullerenes originate from the decay of larger com-pounds subject to energetic excitation (cosmic rays, shocks, VUV irradiation...) followed bystepwise isomerization [34–38]. In this respect, amorphous carbon clusters were suggested toplay a possibly important role [39–41]. Such formation pathways thus question the existenceof products that are concomitantly formed and may be stable in space, spanning the range2etween amorphous and fullerenic carbon clusters.Experimental observation of the possible formation mechanisms will necessarily rely onspectroscopy, and theoretical interpretation of the measured spectra requires appropriatecandidate structures. Unfortunately, the typical approach until now has usually been biasedtowards certain chemical types such as polyaromatic compounds [42], only limited effortbeing devoted to arbitrary or amorphous conformations despite their known astrophysicalrelevance [40, 43–45]. According to the theoretical study by Kent and coworkers [7], 24carbon atoms are needed to produce the first polyaromatic flake, and stable fullerenes ariseabove the size 26, although not necessarily the most stable conformation. Size 60 is notoriousto support fullerenes as their most stable structure. In the absence of known magic numbersin this size range, and as common in clusters physics, 2D and 3D structures are thus likelyto coexist [1].The present article aims at exploring the structural diversity of carbon clusters as a func-tion of their most important feature, namely their internal energy, in the size range wherethey undergo the flake-to-fullerene transition. Although the degree of chemical ordering isobviously expected to vary with increasing excitation energy, the extent of chemical andconformational variety remains undocumented so far, and as far as we are aware our funda-mental study represents the first attempt to construct a library of carbon cluster structuresin an unbiased way, not focusing on the lowest-energy structures only but giving as muchattention 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 se-lected carbon clusters C n across the planar to fullerene transition and containing n = 24, 42,or 60 atoms, using advanced molecular dynamics methods and systematic quenching. Thenumerous isomers thus obtained were sorted and analyzed using a range of structural orderparameters, 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 butrealistic description of the potential energy surface based on the second generation reactiveREBO bond-order potential [46]. Such an approach has already been used in the past tostudy the structure of carbon clusters in the size range up to 55 atoms by Kosimov andcoworkers [9], who predicted a transition from ring structures to graphene flakes occuringabove 18 atoms, without reporting any fullerene. In contrast with earlier computationalinvestigations, our approach here is highly statistical and does not focus on the lowest-3nergy structures. The results thus obtained shed light on the local chemical ordering andglobal structural arrangement of carbon clusters as their energy varies down to low-energypolyaromatic isomers.The article is organized as follows. In the next section, we describe the different computa-tional tools used for generating and characterizing the structural diversity of selected carbonclusters at finite internal energy. The results are presented and discussed in Section III inrelevance to astrophysical implications. In Section IV we present a systematic statisticalanalysis in which we correlate the relative energetic stability of the various conformers totheir structural and chemical features. Some concluding remarks are finally given in SectionV together with ongoing or future extensions. II. METHODSA. Sampling the conformational landscape
The energy landscape of carbon clusters C n is characterized by an exponentially increasingnumber of local minima and transition states with increasing size n , and cannot be sampledexhaustively as soon as this number exceeds a few tens. Furthermore, the barriers separatingthe various local minima or even more distant funnels on the landscape are likely to be highand make traditional simulation methods based on molecular dynamics or Monte Carlomethods poorly efficient. More importantly, we are not interested in addressing the globaloptimization problem specifically by focusing on the lowest-energy structures, but rather toexplore conformations that are physically relevant also at high energies, as could be producede.g. upon photonic or collisional excitation, though still under isolated conditions.We used replica-exchange molecular dynamics (REMD) simulations [47, 48] to circumventthis broken ergodicity issue and achieve a broader sampling of the potential energy surfaceover extended energy ranges.REMD simulations were performed using the LAMMPS program [49] by propagating M trajectories at fixed temperatures T i ( i = 1 . . . M ) and occasionally attempting exchangesbetween the two configurations R i and R j of neighboring replicas i and j = i ±
1. Such anexchange is accepted using the following acceptance rule [50]acc( R i (cid:42)(cid:41) R j ) = min { , exp [( β i − β j )( E ( R i ) − E ( R j ))] } (1)4here we have denoted β i = 1 / ( k B T i ) and E ( R ) the potential energy at configuration R ,discussed below.The efficiency of this exchange process depends on the overlap between the thermaldistributions at temperatures T i and T j , which itself is driven by various factors, primarily thetwo 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 ensurean efficient exploration of configurational space, still below the vaporization temperature.Since the thermal distributions both shift and broaden with increasing temperature, thedifference between successive temperatures T i +1 − T i must also increase with i . Here weemployed a geometrical progression in T i , namely T i +1 = αT i , which is optimal in theharmonic limit [51], to which we have added manually some temperatures to increase theexchange probability.The MD trajectories were integrated with a 0 . × MD steps and the simulation waspropagated for 100 ns. Configurations were periodically saved for further structural analysisand 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 irreversibleprocess preventing the correct sampling of size-selected compounds. Here we used a simplespherical harmonic potential container with a 500 eV / ˚A spring constant and a radius R s whose value was chosen differently for the three cluster sizes n to reach a common densityof 1.7 g/cm , relevant for disordered polyaromatic materials such as soot [52]. The initialconditions of the REMD simulations were taken as the global energy minima for each cluster,also denoted as reference structures. B. Potential energy surface for carbon clusters
The systematic production of large samples of minima at finite temperature, in whichchemical bonds are easily broken and reformed, requires an efficient but chemically realisticmethod for the potential energy surface E ( R ). Here methods with an explicit treatment ofelectronic structure or based on first principles are not practical, especially considering the5ossible multireference character of small carbon clusters that would make the solution tothe electronic problem already cumbersome [53].Only an explicit potential energy surface is currently able to handle the tremendousnumber of configurations gathered with the currently available computational resources,approximate schemes based on tight-binding [54, 55] or density-based tight-binding [56]remaining still too expensive for the present rather large clusters. A few realistic potentialsare available for carbon, which correctly account for bond breaking and formation and thevarious hybridization environments displayed by carbon [46, 57–62]. Here we have chosenthe adaptive the second-generation reactive bond order (REBO) potential of Brenner [46].Brenner type potentials have been notably used to study energetic and mechanical propertiesof nanotubes [63], the formation process of fullerenes [31, 64] but also to describe reactivityand formation of pure carbon clusters or hydrocarbons in the astrophysical context [14, 65],making it a natural choice also for the present investigation. Moreover it has been shownthat the REBO potential gives better results to describe small carbon clusters comparedwith other bond-order potentials [66]. C. Reference structures
Our REMD exploration requires some initial structures. We have thus conducted adistinct search for the lowest-energy configurations of C by comparing various remarkablestructures such as polycyclic, chains and rings conformers. In agreement with the earlierstudy by Kent et al. [7], the lowest-energy structure found for C is the planar, fullydehydrogenated coronene with D h point group. For C the natural reference structure isbuckminsterfullerene with point group I h , besides a set of 1811 alternative but higher energyfullerenes [67]. Finally, C was chosen as an intermediate size cluster with a propensity toform fullerenes. For this cluster, 45 nonequivalent fullerene isomers could be identified [68]and the REBO potential predicts that the most stable isomer is the only fullerene with D symmetry, in accordance with DFT calculations [69, 70] but at variance with the earlierstudy by Kosimov et al. [9] who found a graphene flake as for C .These reference structures are depicted in Figs. 3-C -a), 3-C -a), and 3-C -a) andtheir absolute binding energies obtained with REBO are given in Table I. For buckmin-sterfullerene, the binding energy is reasonably close to the experimentally known value of6 eference structure Point group Binding energy (eV/atom)C D h D I h , C and C obtained with the REBO potential.Cluster T min (K) T max (K) M R s (˚A)C D. Computational details
The number of temperatures allocated for the REMD trajectories, their lowest and highestvalues and the container radius employed in the simulations are given for each system sizein Table II. For all cluster sizes and for each of the M replicas, 25 000 configurations weregenerated. Given the numbers of replicas employed for each system, a total of 300 000instantaneous configurations were kept for further analysis for C , 350 000 instantaneousconfigurations for C and 400 000 configurations for C . E. Identification of stable structures
From the large sets of instantaneous configurations gathered at finite temperature, thelocal minima were obtained by systematic quenching using here the Hessian-free truncatedNewton method as implemented in LAMMPS [49] and disregarding the hard-wall sphericalcontainer. Only connected structures were subsequently kept for further analysis, discon-nected configurations being discarded. Here fragments are identified using a maximumnearest neighbor distance of 1.85 ˚A. In order to speed up the analysis, two locally mini-7ized structures were further considered to be identical if their energy difference lies below0.01 meV, in which case the highest was discarded as well. After this optimization andscreening stage, the numbers of distinct configurations saved for further processing wasequal to 51 901, 240 305, and 236 394 for C , C , and C , respectively. We interpret thisslightly smaller number for the larger cluster size as reflecting its more magic character, thediversity of conformers at a same upper temperature (here 6500 K) being lower than in anon-magic system such as C [72]. However, the energy distributions shown below are es-sentially robust when considering only the first or second half of the configurations harvestedduring the REMD trajectories, suggesting that our samples are statistically representative.Yet, the smaller number of distinct configurations for C would probably require furtherscrutiny, and it could be useful to apply alternative approaches employing systematic localminimization such as basin-sampling [73] in future extensions of this work dealing with evenmore complex systems. F. Structural analysis
Both local and global parameters were used to characterize the various isomers obtainedfor the three clusters sizes. Global parameters provide information about the overall shapeand atomic distribution around the center of mass, while local parameters give insight intothe chemical arrangement at the atomic level.For a n -atom cluster C n with equilibrium configuration R = { r i } for i = 1 . . . n , the 3 × Q is defined from its components Q αβ as [74–77] Q αβ = 1 n n (cid:88) i =1 ( r αi − ¯ r α )( r βj − ¯ r β ) , α, β = 1 , , , (2)where r αi are the Cartesian coordinates ( x , y or z ) of atom i , ¯ r α are the correspondingcoordinates for the cluster center of mass [76]. Three rotationally invariant quantities canbe defined from the tensor Q that respectively measure the geometrical extension, the as-phericity, and the prolateness of the atomic distribution. The squared radius of gyration R g is first defined by R g = 1 n n (cid:88) i =1 ( r i − ¯r ) = Tr Q , (3)where Tr Q stands for the trace of the gyration tensor. The other two quantities are definedfrom the traceless tensor D = Q − I Tr Q , where I is the 3 × A is defined by [77] A = 32 Tr D (Tr Q ) , (4)The asphericity varies from 0 for a purely spherical system to 1 for a perfectly linear struc-ture. Finally, the prolateness S is given for configuration R by S = 27 det D (Tr Q ) , (5)and varies from S = − / oblate structure) to S = 2 for aperfectly linear chain ( prolate structure).The tendency to form hollow, planar or close-packed structures was investigated also fromthe radial density ρ ( r ) around the center of mass, ρ ( r >
0) = 14 πr n (cid:88) i =1 δ ( r − | r i − ¯ r | ) , (6)where δ refers to the Dirac function.Another useful information for pure carbon clusters having a tendency for sp and sp 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 denotethem as R (cid:96) where (cid:96) is the length of the ring.Turning now to local order parameters, the hybridization state of each atom was quan-tified for all stable configurations R of the databases. In many computational studies thatdo not explicitly compute the electronic structure [78–82] hybridization is defined basedon coordination numbers only. However, such a definition cannot account for the chemicalcomplexity and possibly reactive atoms that are under coordinated. Here we use both coor-dination and geometric information to assign hybridization states. More precisely, we define N i the coordination number of atom i , neighbors lying within 1.85 ˚A from atom i , and weevaluate all angles in which atom i is a vertex. There are M i = N i ( N i − / θ k for k = 1 . . . 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 theangles in which it is central are close enough to 180 ◦ :sp: N i = 1 or 2 ,θ k > ◦ ∀ k. (7)9ikewise, the atom is assigned sp hybridization with the following criteriasp : N i = 2 or 3 , ◦ < θ k < ◦ ∀ k, Var( θ k ) < ◦ , (8)where Var( θ k ) denotes the variance of all angles for which i is a vertex. Here again an upperlimit for the angular boundary of 20 ◦ is chosen in such a way as to include atoms involved inhexagonal and pentagonal rings as properly sp . The variance limit ensures that for N i = 3all 4 atoms remain close to a common plane. It was notably chosen to include fullerenestructures (for buckminsterfullerene Var( θ k ) = 6 ◦ ). Finally atom i is considered to be in asp hybridization state accordingly withsp : N i = 4 , ◦ < θ k < ◦ ∀ k, Var( θ k ) < ◦ . (9)Here the variance condition on the angles is required to ensure the atom has a 3D tetrahedralenvironment. With the above definitions, it may occur that a given atom is neither sp, sp or sp , in which case it will be referred to as ambiguous. This notably occurs for ringlikestructures containing too few atoms, hence producing angles that exceed 10 ◦ , one suchexample being illustrated in Fig. 3-C -d). With the present criterion on angles, at least 36atoms would be needed in the perfect ring for the sp assignment to be recovered.Finally, as our last local structural quantity we have determined the pair correlationfunction g ( r ) from g ( r ) = 14 πr ρ S (cid:88) i 0 0.2 0.4 0.6 0.8 1 1.2C C o un t 0 0.2 0.4 0.6 0.8 1 1.2C C o un t 0 0.2 0.4 0.6 0.8 1 1.2C C o un t Energy (eV/atom) FIG. 1. (Color online) Energy distributions of quenched structures for C , C , and C relative tothe corresponding reference structure (lowest isomer). The vertical blue arrows locate the highestenergy isomers found in the REMD simulations. 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 ˚A, respectively. III. RESULTS AND DISCUSSIONS The distributions of local minimum energies obtained for the three clusters C , C andC 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 theranges explored by the REMD simulations, these isomers being depicted at the top of Fig. 3and the ranges themselves highlighted by arrows. For C , the distribution exhibit a singlebroad ( ∼ thedistribution exhibits two peaks at 0.11 and 0.56 eV/atom, showing structures on a largeenergy scale up to ∼ exhibits two peakslocated at 0.11 and 0.77 eV/atom and the overall energies of the structures vary up to ∼ and the two larger cluster sizes and11 % s p c o un t C % s p c o un t C g2 ( Å ) 0 0.2 0.4 0.6 0.8 1 % s p c o un t C FIG. 2. (Color online) Distribution of quenched structures for C , C , and C as a function ofthe squared gyration radius and the fraction of sp carbon atoms. the nature of the two peaks in C and C , the distributions of quenched structures weredetermined as a function of the squared radius of gyration and the proportion of sp carbonatoms (Fig. 2). As for the energy distribution, Fig. 2 shows a single peak for C and twopeaks for C and C . The lowest energy peaks in C and C can be identified to thepeak with larger proportion of sp and smallest squared gyration radius. The highest energypeaks in C and C can be identified to the peak with lowest proportion of sp and largestsquared gyration radius. A direct examination of the structures allows the lowest-energypeak to be assigned to cage-like structures while the highest-energy peak corresponds topretzel-like conformations, as already identified by Kim and Tom´anek [83] who simulatedthe melting of fullerenes using a tight-binding approach [83].In Fig. 3, typical structures are shown for different internal energies. At low energy12 ≤ 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 C andC . A large proportion of atoms are found as sp with only a small number of atoms withambiguous hybridization. As the internal energy increases, a clear increase in the numberof ambiguous hybridization also occurs, concomitant with the appearance of pretzel-likestructures around ∼ and ∼ and C . Athigher energies, i.e. ∼ and ∼ and C , a trend toform 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 footingbased on the different parameters described in the previous section. The energetic stabilityof these isomers can be correlated with their structural properties, starting with the globalstructural order parameters provided by the gyration tensor. The energy-resolved distribu-tions of parameters R g , A and S are shown in Fig. 4 for the three cluster sizes. At a givenenergy the distribution is normalized.At low energy, the three structural parameters smoothly converge to the correspondingvalues in the reference structure. Being shown on the same scale, the square gyrationradii are rather similar for C and C , despite different topologies (planar versus hollowcage). Such differences are better manifested on the asphericity, which vanishes for both cageisomers C and C . From these 2D distributions the behavior of C stands out as differentfrom the two larger clusters, with non monotonic variations in the average R g and A withincreasing 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 decreasein A indicates spherical structures less extended than fully dehydrogenated coronene andcorresponding to the 3D compact structures, as shown in Fig. 3-C -d). Above this energyrange, the distribution of squared gyration radius R g becomes much broader and reaches ∼ 20 ˚A .For the two larger clusters, the global structural indices display more regular variationswith increasing isomer energy, metastable configurations exhibiting larger gyration radii,a greater asphericity and the prolateness remaining low in magnitude but with increasingfluctuations extending mostly to positive values. These fluctuations are most prominentabove 0.5 eV/atom and again convey the greater structural diversity. In contrast, C remains planar until the energy reaches about 0.2 eV/atom, at which stage the rather sharp13 -a) 0 eVC -b) 0.037 eVC -c) 0.174 eVC -d) 0.200 eVC -e) 0.225 eVC -f) 0.521 eV C -a) 0 eVC -b) 0.029 eVC -c) 0.400 eVC -d) 0.536 eVC -e) 0.638 eVC -f) 0.718 eV C -a) 0 eVC -b) 0.033 eVC -c) 0.206 eVC -d) 0.661 eVC -e) 0.731 eVC -f) 0.878 eV FIG. 3. (Color online) Selection of quenched structures obtained for C , C , and C . Carbonatoms are colored accordingly with their identified hybridization state, sp and sp atoms beingshown in red and blue, respectively, and ambiguous atoms in green. No sp atom is present inthese structures. variations in the three structural indicators are consistent with the appearance of morecompact configurations. In particular, the nearly spherical cage-like structures are foundwith very low A in this energy range. 14 R g A S -3 -2 -1 C C C -3 -2 -1 C C C FIG. 4. (Color online) Distributions of the squared radius of gyration R g , asphericity A andprolateness S as a function of the energy per carbon atom for C , C , and C . The blue circlescorrespond to reference structures for C , to the 45 fullerene isomers of C and to the 1812fullerene isomers of C . Above 0.4 eV/atom for C , 0.6 eV/atom for C and 0.7 eV/atom for C , the asphericityand prolateness indices both explore larger values much deviating from the reference struc-tures. Visual inspection indicates that these high energy isomers are usually branched withseveral chains and a few rings only, as depicted in Fig. 3-C -f), Fig. 3-C -f), Fig. 3-C -f).The linear chain isomers, for which A = 1 and S = 2, were indeed found for C and forC but not for C as they do not fit into the spherical container.The radial densities sorted with increasing configuration energy are represented in Fig. 515 r ( Å ) C r ( Å ) -3 -2 -1 C r ( Å ) Energy (eV/atom) FIG. 5. (Color online) Radial densities of the structures as a function of the distance from thecenter of mass, for increasing isomer energy of C , C and C . The blue circles locate thecorresponding peaks in the radial densities of the highly symmetric reference structures. for the three cluster sizes. The ranges of variations in the radial density match those ex-hibited by the global parameters originating from the gyration tensor. In particular, theformation of more compact structures in C below 0.2 eV/atom, the clear cagelike characterin 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 C and 0.6 eV/atom forC are all reflected in Fig. 5. Interestingly, carbon atoms also much more likely occupy thecentral regions when the energy exceeds 0.2–0.6 eV/atom depending on system size, whichis consistent with the loss of hollow structures and the increasing occurrence of chains and16 0 0.2 0.4 0.6 0.8 1energy (eV/atom) 40 60 80 100 120 140 160 180 a n g l e ( d e g r ee s ) N=2 0 0.2 0.4 0.6 0.8 1energy (eV/atom) 10 -4 -3 -2 -1 N=3R R R R FIG. 6. (Color online) Angular distribution for carbon atom with two neighbours (left panel) orthree neighbours (right panel) in C . The horizontal dashed lines represent the angles for idealC rings (60 ◦ ), R rings (108 ◦ ), R (120 ◦ ) and R rings (128.6 ◦ ). branched configurations.Turning to hybridization states, we first show in Fig. 6 and on the example of C onlyhow the angles between connected triplets of carbon atoms are distributed when the centralatom has two or three neighbors, which we anticipate to be potential candidates as sp andsp states, respectively. The angles corresponding to regular rings R , R , R and R arealso shown to highlight the occurrence of such regular polygons in the carbon clusters. Atlow energies, the angular distribution for carbon atoms having two neighbors only is peakedaround 110–120 ◦ , indicating a majority of pentagonal and hexagonal rings. The angles openat energies higher than 0.3 eV/atom, consistently with the formation or chains or largerrings. Five- and six-membered rings, which are predominant at low energies and in fullerenestructures, concomitantly decrease above this same approximate energy threshold. Chainsthemselves show a signature at an angle of 180 ◦ . In 3-coordinated atoms, the distribution isstrongly peaked around 108 ◦ and 120 ◦ which are the expected values for buckminsterfullerenecomposed 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 andmuch more significantly in 3-coordinated atoms as seen through the increasing occurenceof 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 average17 p i ( E ) spsp ambiguous 0 0.2 0.4 0.6 0.8 1 C p i ( E ) p i ( E ) Energy (eV/atom) FIG. 7. (Color online) Probabilities of carbon atoms showing sp, sp , sp or ambiguous hybridiza-tion states in the various configurations of C , C and C , as a function of their energy. distance to the center of mass, which in average approximately equates 83% of the gyrationradius.Based on the above analysis, the effects of excitation energy on the relative proportionof the various hybridation states can now be discussed. The fractions of sp, sp , sp , andambiguous atoms were thus evaluated for all structures in our databases and for the threecluster sizes, the results being shown in Fig. 7. For the three systems, the reference structuresexhibit pure sp hybridization state, as expected for the present polyaromatic isomers. Asenergy increases, the extent of sp hybridization drops around energies where the mostsalient structural changes were noted earlier, that is approximately at 0.2 eV/atom forC , 0.4 eV/atom for C and 0.6 eV/atom for C . However, sp hybridization becomessignificant only at energies higher than these thresholds, while no signature of sp is seenwhatsoever. Ambiguous hybridization states thus populate the intermediate energy rangewhere configurations become less compact, with a rather high amount of curved linear chainsor large rings that do not fall in either of the sp or sp categories.The steady increase in sp hybridization indicates the greater importance of linear chains18n high energy structures. At intermediate energies, many configurations show fewer orshorter such chains, at the expense of rings or curved chains, whose atoms are interpretedas ambiguous until the angle becomes small enough and compatible with sp hybridization.The same qualitative trends are noted in larger clusters, the rise in sp atoms and themaximum in the amount of ambiguous atoms being both shifted to higher energies. Ex-trapolating these trends, we speculate that in even larger clusters the propensity for sp hybridization would become even stronger and more robust against energy excitations, theproportion 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 configu-rations as function of their energy. Figure 8 illustrates this specific property, in average, forthe 3-, 5-, 6-, and 7-membered rings. 4-membered rings, which are much scarcer, were notconsidered in this study. Such rings are occasionally found but deviate significantly from theperfect square and exhibit angles closer to 80 ◦ (and a more standard angle near 120 ◦ ). In thecase of C they can be seen in Fig. 6 for three-coordinated atoms as an horizontal spot withvery 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 withincreasing energy. In both C and C this decrease benefits to pentagonal and heptagonalrings, consistently with the appearance of topological defects such as Stone-Wales pairs.As the internal energy reaches the threshold values where global structural changes takeplace (0.2 eV/atom for C , 0.5 eV/atom for C , 0.7 eV/atom for C ), the numbers of ringshaving 5 or more atoms reaches a minimum and only residual 3-rings are found althoughwith an increasing propensity. The loss of large rings is consistent with more atoms beingpresent as linear chains in such structures, as depicted in Fig. 3. Three-membered ringsbeing energetically rather disfavored, they only appear — occasionally — at the highestexcitation energy as a way to connect linear chains into branched structures.We have finally considered the pair distribution function as another structural indicatorsensitive to the chemical arrangement and in particular hybridization state. Such quantitiesare shown in Fig. 9 for the three cluster sizes and at selected energies. At low energy of0.1 eV/atom and for all sizes, the pair distribution function is dominated by sp carbons. ForC and C , in the range 1.4–1.5 ˚A carbon atoms are involved in hexagonal and heptagonalrings with the two corresponding peaks located at 1.43 and 1.47 ˚A, respectively. Note thatfor buckminsterfullerene the CC distance for a bond located between two hexagonal ring is19 o f c y c l e s R ringsR ringsR ringsR rings 0 2 4 6 8 10 12 C o f c y c l e s o f c y c l e s Energy (eV/atom) FIG. 8. (Color online) Average numbers of R , R , R , and R rings as a function of isomer energyin C , C , and C . The solid circles at zero energy correspond to the values in the referencestructures. reference structure with D point group, the CC distance between two hexagonalrings is 1.41 ˚A, the CC distance between an hexagonal and a pentagonal ring varies between1.45 and 1.46 ˚A. Finally the CC distance between two pentagonal rings varies between 1.48and 1.49 ˚A. For C , the CC distance of the D h reference structure for the inner bondsvaries between 1.44 ˚A and 1.45 ˚A. For the outer bond the distance varies between 1.37 ˚Aand 1.41 ˚A.As energy increases, ambiguous and sp hybridization states become increasingly impor-tant and are manifested by a narrower peak near 1.34 ˚A, and eventually 1.31 ˚A for singlycoordinated atoms. A small residual peak above 1.6 ˚A is also found at very high energy.This peak results from the few carbon atoms involved in three-carbon rings.The fullerene cage C and its cation were both recently observed in the ISM. Identifica-tion of the neutral buckminsterfullerene could be achieved owing to its few and very specificactive bands [84]. The observation of other clusters, even of the fullerene form, shouldbe far more challenging. One important issue is to determine how such highly organized20 R a d i a l d i s t r i b u t i o n f un c t i o n r (Å) 0.1 eV0.3 eV0.5 eV0.7 eV FIG. 9. (Color online) Pair distribution functions as a function of the carbon-carbon distance forC , C , and C at various energies. molecules could be formed in extremely dilute media, notably in the top-down assumptionwhere C would originate from dehydrogenation and subsequent rearrangement of largerpolycyclic aromatic hydrocarbons. Under such a model, our results indicate that the routeto buckminsterfullerene is far from straightforward and encompasses multiple branched iso-mers along the way. It is noteworthy that the high energy configurations found in this workresemble the pretzel phase previously identified by Kim and Tom´anek [83] in their simula-tion of fullerene melting. With the threshold energy reported here, the pretzel phase of C might be present as well in the ISM. IV. STATISTICAL DATA ANALYSIS In the previous section the relation between isomer energy and various structural pa-rameters was highlighted, emphasizing the large configurational diversity arising as energy21 C C R PCC R PCC R PCC R g A S -0.285 -0.114 -0.906 -0.852 -0.965 -0.933TABLE III. Linear ( R ) and partial correlation coefficients (PCC) between various parameters andisomer energy obtained from the databases of minima for C , C and C . increases. Here we adopt a different but complementary point of view by questioning towhich extent energetic stability is statistically related to these geometric parameters. Moreprecisely, we have performed a systematic correlation analysis between isomer energy andthe square gyration radius R g , the asphericity A , the prolateness S and the sp fraction,for the three databases of minima obtained for C , C , and C .The most straightforward way to quantify linear correlation between two sets of data isbased 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 , A , S and sp fraction. We aremostly interested in correlations between energy (acting as X i , our output parameter) andany of the four other quantities (acting as X j for j (cid:54) = i and treated as input parameters).One known issue with Pearson correlation coefficients is that they incorporate possiblystrong correlations between the various input parameters, and an efficient way of removingthese correlations consists of considering partial correlation coefficients (PCCs) [85] instead.The PCC between variables X i and X j involves the full linear correlation matrix R = R ij (including also those elements R jk between geometrical parameters X j and X k for k (cid:54) = i )and reads PCC( { X i } , { X j } ) = − P ij (cid:112) P ii P jj , (11)where P ij is the cofactor of the element R ij in the determinant | R | of matrix R .The linear ( R ij ) and partial correlation coefficients obtained between the four geometricalparameters and the isomer energy are given in Table III for the three cluster sizes. From thesedata, the smaller cluster C appears to behave somewhat differently from the two largerclusters C and C . For these sizes that support fullerene-type isomers, linear correlations22een from R quantities are highest with the sp fraction, the negative sign obtained for R being consistent with the intuitive observation that isomers lowest in energy have the highestsp fraction. The three shape parameters always correlate positively but with Pearsoncoefficients never exceeding 0.65, such a poor linear correlation appears consistent with thescatter plots reported in Fig. 4.Removing the correlations between the various parameters in the PCC quantities confirmsthis trend for the two larger clusters, and highlights the sp fraction as the most sensitiveparameter causing relative energetic stability among conformers. However, for C thisquantity performs not as well, and none of the quantities shows a strong partial correlationwith energy. This can be understood because for C , within the main peak of the energydistribution (between 0.2 and 0.6 eV/atom) the sp fraction remains relatively constant.The negative sign obtained with the square gyration radius is related with the markedlydifferent distributions displayed in Fig. 4, where for C structures near 0.2 eV/atom arefound to be more compact than the planar global minimum, at variance with the behaviornoted for the two larger clusters.The present statistical analysis thus suggests that, for fullerenes, the fraction of sp atomsis the main factor responsible for the energetic stability of the various conformers, geometricshape parameters appearing of lesser importance. Proceeding further, we have attemptedto represent the energy X as a simple function of the other parameters { X i } , 1 ≤ i ≤ (cid:101) E = a + (cid:88) i a i X i + (cid:88) i The authors gratefully acknowledge financial support by the Agence Nationale de laRecherche (ANR) grant ANR-16-CE29-0025. This work was supported by grants fromR´egion Ile-de-France and by the GDR EMIE 3533. [1] G. von Helden, M. Hsu, P. Kemper, and M. Bowers, J. Chem. Phys. , 3835 (1991).[2] G. von Helden, M. Hsu, N. Gotts, and M. Bowers, J. Phys. Chem. , 8182 (1993).[3] M. J. L´opez, P. A. Marcos, A. Rubio, and J. A. Alonso, Z. Phys. D: At. Mol. Clusters ,385 (1997).[4] A. von Orden and R. Saykally, Chem. Rev. , 2313 (1998). 5] H. Handschuh, G. Gantefor, B. Kessler, P. S. Bechthold, and W. Eberhardt, Phys. Rev. Lett. , 1095 (1995).[6] C. Lifshitz, Int. J. Mass Spectrom. , 423 (2000).[7] P. Kent, M. Towler, R. Needs, and G. Rajagopal, Phys. Rev. B , 15394 (2000).[8] Y. Ueno and S. Saito, Phys. Rev. B , 085403 (2008).[9] D. P. Kosimov, A. A. Dzhurakhalov, and F. M. Peeters, Phys. Rev. B , 195414 (2010).[10] C. J. E. Straatsma, M. I. Fabrikant, G. E. Douberly, and H. J. Lewandowski, J. Chem. Phys. , 124201 (2017).[11] N. Kono, R. Suzuki, T. Furukawa, J. Matsumoto, H. Tanuma, H. Shiromaru, T. Azuma, andK. Hansen, Phys. Rev. A , 063434 (2018).[12] B. G. A. Brito, G.-Q. Hai, and L. Cˆandido, Phys. Rev. A , 062508 (2018).[13] D. Tom´anek and M. A. Schluter, Phys. Rev. Lett. , 2331 (1991).[14] R. Delaunay, M. Gatchell, P. Rousseau, A. Domaracka, S. Maclot, Y. Wang, M. H. Stockett,T. Chen, L. Adoui, M. Alcam´ı, F. Mart´ın, H. Zettergren, H. Cederquist, and B. A. Huber,J. Phys. Chem. Lett. , 1536 (2015).[15] A. Simon, J.-P. Champeaux, M. Rapacioli, P. Moretto Capelle, F. X. Gad´ea, and M. Sence,Theor. Chem. Acc. , 106 (2018).[16] A. Marciniak, V. Despr´e, T. Barillot, A. Rouz´ee, M. C. E. Galbraith, J. Klei, C. H. Yang,C. T. L. Smeenk, V. Loriot, S. N. Reddy, A. G. G. M. Tielens, S. Mahapatra, A. I. Kuleff,M. J. J. Vrakking, and F. L´epine, Nature Comm. , 7909 (2015).[17] M. Ji, L. Chen, R. Bredy, C. Ortega, C. Joblin, A. Cassimi, and S. Martin, J. Chem. Phys. , 044301 (2017).[18] J. F. Zhen, S. R. Castillo, C. Joblin, G. Mulas, H. Sabbah, A. Giuliani, L. Nahon, S. Martin,J.-P. Champeaux, and P. M. Mayer, Astrophys. J. , 113 (2016).[19] P. F. Bernath, K. H. Hinkle, and J. J. Keady, Science , 562 (1989).[20] J. Maier, N. Lakin, G. Walker, and D. Bohlender, Astrophys. J. , 267 (2001).[21] K. Sellgren, M. Werner, J. Ingalls, J. Smith, T. Carleton, and C. Joblin, Astrophys. J. Lett. , L54 (2010).[22] J. Cami, J. Bernard-Salas, E. Peeters, and S. Malek, Science , 1180 (2010).[23] O. Bern´e, G. Mulas, and C. Joblin, Astron. Astrophys. , L4 (2013).[24] A. L´eger and J. L. Puget, Astron. Astrophys. , L5 (1984). 25] L. J. Allamandola, A. G. G. M. Tielens, and J. R. Barker, Astrophys. J. , L25 (1985).[26] G. Sloan, M. Jura, W. Duley, K. Kraemer, J. Bernard-Salas, W. Forrest, B. Sargent, A. Li,D. Barry, C. Bohac, D. Watson, and J. R. Houck, Astrophys. J. , 1144 (2007).[27] T. Pino, E. Dartois, A. T. Cao, Y. Carpentier, T. Chamaille, R. Vasquez, A. P. Jones,L. d’Hendecourt, and P. Br´echignac, Astron. Astrophys. , 665 (2008).[28] B. Acke, J. Bouwman, A. Juhasz, T. Henning, M. van den Ancker, G. Meeus, A. Tielens, andL. Waters, Astrophys. J. , 558 (2010).[29] G. C. Sloan, E. Lagadec, A. A. Zijlstra, K. E. Kraemer, A. P. Weis, M. Matsuura, K. Volk,E. Peeters, W. W. Duley, J. Cami, J. Bernard-Salas, F. Kemper, and R. Sahai, Astrophys.J. , 28 (2014).[30] P. Pilleri, C. Joblin, F. Boulanger, and T. Onaka, Astron. Astrophys. , A16 (2015).[31] Y. Yamaguchi and S. Maruyama, Chem. Phys. Lett. , 336 (1998).[32] S. Irle, G. Zheng, M. Elstner, and K. Morokuma, Nano Lett. , 1657 (2003).[33] G. Zheng, S. Irle, and K. Morokuma, J. Chem. Phys. , 014708 (2005).[34] O. Bern´e, J. Montillaud, and C. Joblin, Astron. Astrophys. , A133 (2015).[35] A. Chuvilin, U. Kaiser, E. Bichoutskaia, N. A. Besley, and A. N. Khlobystov, Nature Chem. , 450 (2010).[36] J. Montillaud, C. Joblin, and D. Toublanc, Astron. Astrophys. , A15 (2013).[37] F. Pietrucci and W. Andreoni, J. Chem. Theory Comput. , 913 (2014).[38] J. Zhang, F. L. Bowles, D. W. Bearden, W. K. Ray, T. Fuhrer, Y. Ye, C. Dixon, K. Harich,R. F. Helm, M. M. Olmstead, O. A. L. Balch, and H. C. A. Dorn, Nature Chem. , 880(2013).[39] R. C. Powles, N. A. Marks, and D. W. M. Lau, Phys. Rev. B , 075430 (2009).[40] A. P. Jones, M. Kohler, N. Ysard, M. Bocchio, and L. Verstraete, Astron. Astrophys. ,A46 (2017).[41] A. S. Sinitsa, I. V. Lebedeva, A. M. Popov, and A. A. Knizhnik, J. Phys. Chem. C ,13396 (2017).[42] C. W. Bauschlicher, Jr., A. Ricca, C. Boersma, and L. J. Allamandola, Astrophys. J. Suppl.S. , 32 (2018).[43] W. W. Duley and A. L. Hu, Astrophys. J. , 115 (2012). 44] A. P. Jones, L. Fanciullo, M. K¨ohler, L. Verstraete, V. Guillet, M. Bocchio, and N. Ysard,Astron. Astrophys , A62 (2013).[45] T. Bout´eraon, E. Habart, N. Ysard, A. P. Jones, E. Dartois, and T. Pino, arXiv e-prints(2019), arXiv:1901.07332.[46] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, J.Phys. Cond. Mat. , 783 (2002).[47] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. , 2607 (1986).[48] E. Marinari and G. Parisi, Europhys. Lett. , 451 (1992).[49] S. Plimpton, J. Comput. Phys. , 1 (1995).[50] Y. Sugita and Y. Okamoto, Chem. Phys. Lett. , 141 (1999).[51] A. Mitsutake, Y. Sugita, and Y. Okamoto, Biopolymers , 96 (2001).[52] J. Yon, A. Bescond, and F.-X. Ouf, J. Aerosol Sci. , 28 (2015).[53] L. Montagnon and F. Spiegelman, J. Chem. Phys. , 084111 (2007).[54] Y. Wang and C. H. Mak, Chem. Phys. Lett. , 37 (1995).[55] N.-T. Van-Oanh, P. Parneix, and P. Br´echignac, J. Phys. Chem. A , 10144 (2002).[56] M. Rapacioli, A. Simon, C. Marshall, J. Cuny, D. Kokkin, F. Spiegelman, and C. Joblin, J.Phys. Chem. A , 12845 (2015).[57] J. Tersoff, Phys. Rev. Lett. , 2879 (1988).[58] D. W. Brenner, Phys. Rev. B , 9458 (1990).[59] S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem. Phys. , 6472 (2000).[60] N. A. Marks, Phys. Rev. B , 035401 (2000).[61] L. M. Ghiringhelli, J. H. Los, A. Fasolino, and E. J. Meijer, Phys. Rev. B , 214103 (2005).[62] X. W. Zhou, D. K. Ward, and M. E. Foster, J. Comput. Chem. , 1719 (2015).[63] B. Ni, S. B. Sinnott, P. T. Mikulski, and J. A. Harrison, Phys. Rev. Lett. , 205505 (2002).[64] S. Makino, T. Oda, and Y. Hiwatari, J. Phys. Chem. Solids , 1845 (1997).[65] N. Patra, P. Kr´al, and H. R. Sadeghpour, Astrophys. J. , 6 (2014).[66] S. K. Lai, I. Setiyawati, T. W. Yen, and Y. H. Tang, Theor. Chem. Acc. , 20 (2016).[67] R. Sure, A. Hansen, P. Schwerdtfeger, and S. Grimme, Phys. Chem. Chem. Phys. , 14296(2017).[68] P. Fowler and D. Manolopoulos, An Atlas of Fullerenes (Clarendon Press, Oxford, 1995).[69] G. Sun, M. C. Nicklaus, and R. H. Xie, J. Phys. Chem. A , 4617 (2005). 70] E. Malolepsza, Y. P. Lee, H. A. Witek, S. Irle, C. F. Lin, and H. M. Hsieh, Int. J. QuantumChem. , 1999 (2009).[71] M. S. Dresselhaus, G. Dresselhaus, and P. C. Eklund, Science of Fullerenes and CarbonNanotubes (Academic Press, California, 1996).[72] D. J. Wales, Energy Landscapes (Cambridge University Press, Cambridge, 2003).[73] T. V. Bogdan, F. Calvo, and D. J. Wales, J. Chem. Phys. , 044102 (2006).[74] K. ˇSolc and W. H. Stockmayer, J. Chem. Phys. , 2756 (1971).[75] D. N. Theodorou and U. W. Suter, Macromolecules , 1206 (1985).[76] V. Blavatska and W. Janke, J. Chem. Phys. , 184903 (2010).[77] F. Calvo, F. Chirot, F. Albrieux, J. Lemoine, Y. O. Tsybin, P. Pernot, and P. Dugourd, J.Am. Soc. Mass Spectrom. , 1279 (2012).[78] G. Galli, R. M. Martin, R. Car, and M. Parrinello, Phys. Rev. Lett. , 555 (1989).[79] G. Galli, R. M. Martin, R. Car, and M. Parrinello, Phys. Rev. B , 7470 (1990).[80] C. Z. Wang, K. M. Ho, and C. T. Chan, Phys. Rev. B , 14835 (1993).[81] B. Bhattarai and D. A. Drabold, Carbon , 532 (2017).[82] V. S. Dozhdikov, A. Y. Basharin, P. R. Levashov, and D. V. Minakov, J. Chem. Phys. ,214302 (2017).[83] S. Kim and D. Tom´anek, Phys. Rev. Lett. , 2418 (1994).[84] A. Omont, Astron. Astrophys. , A52 (2016).[85] D. M. Hamby, Environ. Monit. Assess. , 135 (1994).[86] A. K. Rapp´e and W. A. Goddard, J. Phys. Chem. , 3358 (1991).[87] M. Elstner and G. Seifert, Phil. Trans. R. Soc. A , 20120483 (2014)., 20120483 (2014).