The Effect of Cluster Size on the Intra-Cluster Ionic Polymerization Process

Polyaromatic hydrocarbons (PAHs) are widespread in the interstellar medium (ISM). The abundance and relevance of PAHs call for a clear understanding of their formation mechanisms, which, to date, have not been completely deciphered. Of particular interest is the formation of benzene, the basic building block of PAHs. It has been shown that the ionization of neutral clusters can lead to an intra-cluster ionic polymerization process that results in molecular growth. Ab-initio molecular dynamics (AIMD) studies in clusters consisting of 3–6 units of acetylene modeling ionization events under ISM conditions have shown maximum aggregation of three acetylene molecules forming bonded C6H6+ species; the larger the number of acetylene molecules, the higher the production of C6H6+. These results lead to the question of whether clusters larger than those studied thus far promote aggregation beyond three acetylene units and whether larger clusters can result in higher C6H6+ production. In this study, we report results from AIMD simulations modeling the ionization of 10 and 20 acetylene clusters. The simulations show aggregation of up to four acetylene units producing bonded C8H8+. Interestingly, C8H8+ bicyclic species were identified, setting a precedent for their astrochemical identification. Comparable reactivity rates were shown with 10 and 20 acetylene clusters.


Introduction
Polyaromatic hydrocarbons (PAHs) are widespread in the interstellar medium (ISM) [1][2][3]. They are believed to be ubiquitous and abundant in space, with interstellar IR spectra providing evidence of their presence in different areas by identifying their characteristic features in the spectra [4][5][6][7]. PAHs serve as a bridge between small organic molecules and large carbonaceous materials [8] and thus play an important role in the chemical evolution of the ISM.
Routes for the formation of PAHs and other organic molecules known to be present in the interstellar medium are not thoroughly understood. Various mechanisms for PAH formation and growth have been suggested on the basis of, for example, radical reactions [9][10][11][12][13][14][15]. One such example is the well-known hydrogen abstraction C 2 H 2 addition (HACA) mechanism, which is based on two steps: the first is the formation of the aryl radical, and the second is the addition of acetylene (C 2 H 2 ). The repetition of these two steps results in PAH formation. Other studies suggest the formation of naphthalene via the reaction of the phenyl radical and vinylacetylene or, alternatively, via the formation of naphthalene ions through ion-molecule reactions between benzene cations (C 6 H 6 ) + and ethynyl radicals [16]. The formation of benzene, the basic building block of PAHs, has drawn a great deal of attention. Different pathways for its formation have been proposed, depending on the physical and chemical conditions [17][18][19][20]. Jones et al. suggested a barrierless route for benzene formation via the reaction of the ethynyl radical and 1,3-butadiene, which is relevant to areas such as the Taurus Molecular Cloud (TMC-1) [21]. In several astrochemical models, routes for the formation of benzene are based on ion-molecule reactions [15,[22][23][24][25]. Recently, Zhao et al. showed that molecular growth processes can occur via the reaction between two stabilized propargyl radicals at high temperature and diluted environments. Such a reaction can lead to the formation of benzene molecules, among other species [26].
To understand the intra-cluster ionization process at the molecular level, several studies performed ab-initio molecular dynamics (AIMD) simulations modeling different neutral cluster compositions. In agreement with experimental work, AIMD simulations have demonstrated that the molecular growth process occurs and that covalently bonded structures are formed [38,[41][42][43]. It has been shown that, upon ionization, part of the cluster forms covalently bonded core structures, while the remaining molecules serve as spectators of the process. The role of spectator molecules is important for the growth process, as these molecules change the potential energy surfaces and thus play a catalytic role in the formation of various core structures. Moreover, spectator molecules enable stabilization routes, as they can dissipate the excess energy via evaporation [38,42]. For example, upon ionization of acetylene clusters containing up to six acetylene units, the formation of bonded C 4 H 4 + and C 6 H 6 + including benzene cations was also observed [38]. Ionization of small clusters (up to five units) containing acetylene and hydrogen cyanide also led to molecular growth, in which three and four units bonded to form structures on the potential energy surfaces of C 6 H 6 + , C 5 H 5 N + , C 7 H 7 N + , and C 8 H 8 + . Likewise, when the clusters are built from acetylene and cyanoacetylene units, building blocks that can be found in TMC-1, the ionization process leads to bonding between three and four units [41]. Many of the structures that are formed are aromatic structures and contain a nitrogen atom within the ring or as a side chain; thus, the structures are important from an astrobiological point of view. Among the formed structures we observed was the benzonitrile cation, which was recently identified in its neutral form in TMC-1 [44].
AIMD results from previous studies demonstrate the importance of the composition and size of the cluster in promoting aggregation to yield chemically bonded structures. In reference to pure acetylene clusters, in particular, it has been demonstrated that the larger clusters (five and six acetylene units) enable higher rates of bonded C 6 H 6 + formation in comparison with the smaller ones (three and four acetylene units) [38]. Here, we focus on the study of pure acetylene clusters to understand the influence of a large number (10 and 20, significantly larger than the six clusters studied to date) of acetylene molecules in the aggregation process that produces (C 2 H 2 ) n + species following ionization of the van der Waals clusters. Specifically, our aim was to investigate the effect of cluster size on the percentage of C 6 H 6 + produced in the larger clusters vs. the smaller ones. Additionally, we examined the extent of molecular growth in the larger clusters.

Results
In order to address the aforementioned objectives, we studied large van der Waals clusters containing 10 and 20 units of acetylene. We built 40 structures: 20 decamers and 20 eicosamers (10 and 20 acetylene units, respectively), considering different relative orientations of the acetylene molecules that maximize the CH-π interaction between hydrogen and the π clouds of neighboring acetylenes, and we then optimized them to obtain minima on the neutral potential energy surface (PES) [40]. minima on the neutral potential energy surface (PES) [40]. Figure 1 shows examples of a decamer structure ( Figure 1a) and an eicosamer structure (Figure 1b). The complete set of coordinates corresponding to the optimized AIMD starting structures are available in the Supplementary materials. In order to model the ionization process, the optimized neutral structures were utilized as starting structures in the molecular dynamic simulations on the cationic PES. The neutral cluster geometries are nonoptimal after vertical ionization, and they relax on the cationic PES, yielding bonded structures. Figure 2 presents the structures after optimization on the cationic PES. In both cases (decamer and eicosamer), the clusters relax to core-bonded structures in the C4H4 + , C6H6 + , and C8H8 + PES (Figure 2a-c for the decamer and Figure 2d-f for the eicosamer) without any barriers. The core structures are solvated by additional nonbonded acetylene molecules. We note that the largest core structure that is formed upon optimization is a bonded C8H8 + . The coordinates of all 40 optimized cationic isomers are reported in the SI. (a-c) are decamer clusters optimized on the cationic PES with C4H4 + , C6H6 + , and C8H8 + corebonded structures, respectively; (d-f) correspond to the eicosamer clusters that also show C4H4 + , C6H6 + , and C8H8 + core-bonded structures. The neutral cluster geometries are nonoptimal after vertical ionization, and they relax on the cationic PES, yielding bonded structures. Figure 2 presents the structures after optimization on the cationic PES. In both cases (decamer and eicosamer), the clusters relax to core-bonded structures in the C 4 H 4 + , C 6 H 6 + , and C 8 H 8 + PES (Figure 2a-c for the decamer and Figure 2d-f for the eicosamer) without any barriers. The core structures are solvated by additional nonbonded acetylene molecules. We note that the largest core structure that is formed upon optimization is a bonded C 8 H 8 + . The coordinates of all 40 optimized cationic isomers are reported in the SI. minima on the neutral potential energy surface (PES) [40]. Figure 1 shows examples of a decamer structure ( Figure 1a) and an eicosamer structure (Figure 1b). The complete set of coordinates corresponding to the optimized AIMD starting structures are available in the Supplementary materials. In order to model the ionization process, the optimized neutral structures were utilized as starting structures in the molecular dynamic simulations on the cationic PES. The neutral cluster geometries are nonoptimal after vertical ionization, and they relax on the cationic PES, yielding bonded structures. Figure 2 presents the structures after optimization on the cationic PES. In both cases (decamer and eicosamer), the clusters relax to core-bonded structures in the C4H4 + , C6H6 + , and C8H8 + PES (Figure 2a-c for the decamer and Figure 2d-f for the eicosamer) without any barriers. The core structures are solvated by additional nonbonded acetylene molecules. We note that the largest core structure that is formed upon optimization is a bonded C8H8 + . The coordinates of all 40 optimized cationic isomers are reported in the SI.

AIMD Simulations
Following ionization and the formation of the core structures, the system contains a large amount of excess energy (~5 eV) that can be utilized to cross barriers on the PES. In order to study the evolution of the systems over time and the potential formation of different products upon ionization, we used AIMD. The results of the AIMD simulations shed light on the structures that can be formed at the end of this process and the relative probability of their formation. The structures of the most relevant mono-and bicyclic molecules obtained are presented, and a comparison is made between the decamers and eicosamers as well as smaller acetylene van der Waals clusters previously studied by means of AIMD under similar conditions [38]. We analyze the extent of the molecular growth as a function of the number of acetylene molecules in a cluster and conclude the relevance of the present results in understanding PAH formation routes and astrochemical identification of cyclic molecules.

AIMD of the Acetylene Decamer Clusters
For each of the 20 starting structures, 30 trajectories were run, totaling 600 trajectories. All trajectories were performed for 2.4 ps. For trajectories ending in structures where further potential rearrangements or growth were suspected, the simulation time was extended up to 4.8 ps. The resulting core structures are shown in Figure 3. From the 600 trajectories for this cluster size, we found that 56.5% of the trajectories led to a molecular product of the reaction of 2 units of acetylene to yield bonded C 4 H 4 + . The vast majority of the C 4 H 4 + species produced corresponded to cyclobutadiene (Figure 3a), while a minor percentage was methylenecyclopropane (Figure 3b). In other trajectories, which represent 40.8% of the trajectories, three units of acetylene react to yield bonded C 6 H 6 + . Multiple C 6 H 6 + products were identified. Among them were benzene ( Figure 3h) and its conformational isomers fulvene (Figure 3e), dewar-benzene (Figure 3f), and benzvalene (Figure 3g), all of which are known to easily interconvert into benzene. In a previous study, these products were also identified in the products from the ionization of clusters up to hexamer, although in smaller percentages [38]. Interestingly, here we observed a higher percentage of the trajectories leading to bonded C 6 H 6 + , which reinforces the important role played by the cluster environment. In the remaining~2.7% of the trajectories, we observed the aggregation of four acetylene molecules to yield bonded C 8 H 8 + . While other cluster compositions (namely, mixed acetylene with HCN and cyanoacetylene) enable bonding between four units [41,43], C 8 H 8 + species were not obtained from the acetylene hexamers previously studied, clearly demonstrating that more than six acetylene units are required to observe aggregation beyond C 6 structures from a pure acetylene cluster [38,42]. Again, this emphasizes the role of the cluster environment not only in the stabilization of the products but also in their formation.

AIMD Simulations
Following ionization and the formation of the core structures, the system contains a large amount of excess energy (~5 eV) that can be utilized to cross barriers on the PES. In order to study the evolution of the systems over time and the potential formation of different products upon ionization, we used AIMD. The results of the AIMD simulations shed light on the structures that can be formed at the end of this process and the relative probability of their formation. The structures of the most relevant mono-and bicyclic molecules obtained are presented, and a comparison is made between the decamers and eicosamers as well as smaller acetylene van der Waals clusters previously studied by means of AIMD under similar conditions [38]. We analyze the extent of the molecular growth as a function of the number of acetylene molecules in a cluster and conclude the relevance of the present results in understanding PAH formation routes and astrochemical identification of cyclic molecules.

AIMD of the Acetylene Decamer Clusters
For each of the 20 starting structures, 30 trajectories were run, totaling 600 trajectories. All trajectories were performed for 2.4 ps. For trajectories ending in structures where further potential rearrangements or growth were suspected, the simulation time was extended up to 4.8 ps. The resulting core structures are shown in Figure 3. From the 600 trajectories for this cluster size, we found that 56.5% of the trajectories led to a molecular product of the reaction of 2 units of acetylene to yield bonded C4H4 + . The vast majority of the C4H4 + species produced corresponded to cyclobutadiene (Figure 3a), while a minor percentage was methylenecyclopropane (Figure 3b). In other trajectories, which represent 40.8% of the trajectories, three units of acetylene react to yield bonded C6H6 + . Multiple C6H6 + products were identified. Among them were benzene (Figure 3h) and its conformational isomers fulvene (Figure 3e), dewar-benzene (Figure 3f), and benzvalene (Figure 3g), all of which are known to easily interconvert into benzene. In a previous study, these products were also identified in the products from the ionization of clusters up to hexamer, although in smaller percentages [38]. Interestingly, here we observed a higher percentage of the trajectories leading to bonded C6H6 + , which reinforces the important role played by the cluster environment. In the remaining ~2.7% of the trajectories, we observed the aggregation of four acetylene molecules to yield bonded C8H8 + . While other cluster compositions (namely, mixed acetylene with HCN and cyanoacetylene) enable bonding between four units [41,43], C8H8 + species were not obtained from the acetylene hexamers previously studied, clearly demonstrating that more than six acetylene units are required to observe aggregation beyond C6 structures from a pure acetylene cluster [38,42]. Again, this emphasizes the role of the cluster environment not only in the stabilization of the products but also in their formation.  The largest bonded structures obtained corresponded to the aggregation of four acetylene units to yield C 8 H 8 + isomers, as shown in Figure 4. Among these products, we detected structures that contained a three-member (Figure 4a) or four-member (Figure 4b) ring, as observed in the C 4 H 4 + and C 6 H 6 + structures. Trajectories leading to C 8 H 8 + were extended for an additional 2.4 ps. At these longer times, C 8 H 8 + products further react, yielding bicyclic products with five-and three- (Figure 4c), five-and four- (Figure 4d), five-and five- (Figure 4e), and six-and four- (Figure 4f) member ring structures. In the literature, C 8 H 8 + is attributed to a benzene cation complexed to an acetylene molecule, which we also observed in our simulations [33]. The PES for the formation of the structure in Figure 4f is presented in Figure 5, along with the corresponding time of occurrence of each structure during the simulation time. At the first stages of the simulation, two acetylene molecules were bonded together, and a three-member ring structure ( Figure 5, 252 fs) was formed. This resembled the association product of acetylene and its cation observed by Bera et al. [45]. Afterward, an additional bond formed with a neighboring acetylene molecule, resulting in a bicyclic structure ( Figure 5, 584 fs), which later reorganized into a three-member cyclic structure with a long chain ( Figure 5, 1067 fs). A similar structure was also seen when all four acetylene molecules were bonded ( Figure 5, 1592 fs). We observed a transition to a four-member ring ( Figure 5, 1729 fs), which reorganized to produce the final product. The largest bonded structures obtained corresponded to the aggregation of four acetylene units to yield C8H8 + isomers, as shown in Figure 4. Among these products, we detected structures that contained a three-member (Figure 4a) or four-member (Figure 4b) ring, as observed in the C4H4 + and C6H6 + structures. Trajectories leading to C8H8 + were extended for an additional 2.4 ps. At these longer times, C8H8 + products further react, yielding bicyclic products with five-and three- (Figure 4c), five-and four- (Figure 4d), five-and five- (Figure 4e), and six-and four- (Figure 4f) member ring structures. In the literature, C8H8 + is attributed to a benzene cation complexed to an acetylene molecule, which we also observed in our simulations [33]. The PES for the formation of the structure in Figure 4f is presented in Figure 5, along with the corresponding time of occurrence of each structure during the simulation time. At the first stages of the simulation, two acetylene molecules were bonded together, and a three-member ring structure ( Figure 5, 252 fs) was formed. This resembled the association product of acetylene and its cation observed by Bera et al. [45]. Afterward, an additional bond formed with a neighboring acetylene molecule, resulting in a bicyclic structure ( Figure 5, 584 fs), which later reorganized into a three-member cyclic structure with a long chain ( Figure 5, 1067 fs). A similar structure was also seen when all four acetylene molecules were bonded ( Figure 5, 1592 fs). We observed a transition to a four-member ring ( Figure 5, 1729 fs), which reorganized to produce the final product.
A bicyclic C8H8 + molecule with five-and three-member rings (Figure 4c) was observed in two of the trajectories. One of them dissociated in a later stage into C5H5 and c-C3H3 + , which was experimentally detected [46], and its direct observation in the ISM is only feasible via ro-vibrational transitions because of the absence of pure rotational transitions [46].

Comparison of the Growth Process in Different Cluster Sizes
Starting from 20 eicosamer optimized structures, we ran a total of 600 trajectories at 2.4 ps. The results for this cluster size were similar to those obtained for the decamer clusters; no additional structures were found when the AIMD started from the eicosamer clusters. Figure 6 shows a comparison of the number of core structures obtained at the end of The minima energies are reported in black, and the transition state energies are in red. The corresponding time of each structure during the molecular dynamics trajectory is noted in blue.
A bicyclic C 8 H 8 + molecule with five-and three-member rings (Figure 4c) was observed in two of the trajectories. One of them dissociated in a later stage into C 5 H 5 and c-C 3 H 3 + , which was experimentally detected [46], and its direct observation in the ISM is only feasible via ro-vibrational transitions because of the absence of pure rotational transitions [46].

Comparison of the Growth Process in Different Cluster Sizes
Starting from 20 eicosamer optimized structures, we ran a total of 600 trajectories at 2.4 ps. The results for this cluster size were similar to those obtained for the decamer clusters; no additional structures were found when the AIMD started from the eicosamer clusters. Figure 6 shows a comparison of the number of core structures obtained at the end of the AIMD simulations of different cluster sizes consisting of 6, 10, and 20 acetylene units. In the case of the hexamer clusters, it can be seen that the majority of trajectories formed C 4 H 4 + core structures, while the rest of the trajectories formed C 6 H 6 + core structures. No C 8 H 8 + structure was observed. In the case of the decamer and eicosamer clusters, more C 6 H 6 + was obtained with respect to the hexamer clusters. Additionally, these cluster sizes enabled the formation of C 8 H 8 + structures that were not observed in the smaller clusters. In fact, a previous photoionization mass spectrometry study demonstrated that C 8 H 8 + can be obtained from acetylene and is identified by a peak at m/z = 104 [38]. The percentages for the formation of different structures were similar for the decamers and eicosamers; no significant differences were observed between them. The C6H6 + /C4H4 + ratio is a way to express the capacity of molecular growth in acetylene clusters. A previous study found that the C6H6 + /C4H4 + ratio was ~0.3 for acetylene dimers and trimers, and the C6H6 + /C4H4 + ratio was ~0.5 for tetramers and pentamers [38]. Here, we found that the C6H6 + /C4H4 + ratio was 0.7 for decamers and eicosamers. Figure 7 summarizes the trend of the C6H6 + /C4H4 + ratio for pure acetylene clusters as a function of the number of acetylenes and is an adaptation of a figure from a previous work [38] in order to add the results from the present study. Again, the results demonstrate that a growing number of acetylene spectator molecules in a cluster promotes the molecular growth process by the process described above. As can be seen from the figure, decamers and eicosamers present comparable reactivity; eicosamers do not lead to further aggregation (e.g., C10H10 + and beyond) despite having double the number of acetylenes. This suggests that the addition of acetylene molecules to the large clusters' environment (beyond the decamers studied here) does not affect the environment substantially and does not prompt additional growth. We note, however, that it is also plausible that clusters containing species other than acetylene can be involved in catalyzing the formation of structures larger than C8H8 + . The C 6 H 6 + /C 4 H 4 + ratio is a way to express the capacity of molecular growth in acetylene clusters. A previous study found that the C 6 H 6 + /C 4 H 4 + ratio was~0.3 for acetylene dimers and trimers, and the C 6 H 6 + /C 4 H 4 + ratio was~0.5 for tetramers and pentamers [38]. Here, we found that the C 6 H 6 + /C 4 H 4 + ratio was 0.7 for decamers and eicosamers. Figure 7 summarizes the trend of the C 6 H 6 + /C 4 H 4 + ratio for pure acetylene clusters as a function of the number of acetylenes and is an adaptation of a figure from a previous work [38] in order to add the results from the present study. Again, the results demonstrate that a growing number of acetylene spectator molecules in a cluster promotes the molecular growth process by the process described above. As can be seen from the figure, decamers and eicosamers present comparable reactivity; eicosamers do not lead to further aggregation (e.g., C 10 H 10 + and beyond) despite having double the number of acetylenes. This suggests that the addition of acetylene molecules to the large clusters' environment (beyond the decamers studied here) does not affect the environment substantially and does not prompt additional growth. We note, however, that it is also plausible that clusters containing species other than acetylene can be involved in catalyzing the formation of structures larger than C 8 H 8 + .
summarizes the trend of the C6H6 + /C4H4 + ratio for pure acetylene clusters as a function of the number of acetylenes and is an adaptation of a figure from a previous work [38] in order to add the results from the present study. Again, the results demonstrate that a growing number of acetylene spectator molecules in a cluster promotes the molecular growth process by the process described above. As can be seen from the figure, decamers and eicosamers present comparable reactivity; eicosamers do not lead to further aggregation (e.g., C10H10 + and beyond) despite having double the number of acetylenes. This suggests that the addition of acetylene molecules to the large clusters' environment (beyond the decamers studied here) does not affect the environment substantially and does not prompt additional growth. We note, however, that it is also plausible that clusters containing species other than acetylene can be involved in catalyzing the formation of structures larger than C8H8 + .

Materials and Methods
Every calculation in this manuscript was performed using Q-Chem 5.1 software [47]. Optimizations of neutral and cationic structures were accomplished with the wB97X-V functional [48] and cc-pVTZ basis set [49]. Using the same level of theory, we performed vibrational frequency calculations to categorize the structures as minima or saddle points on the PES.
Starting from the optimized minima structures of isomers of clusters consisting of 10 or 20 units of acetylene, we ran AIMD simulations on the cationic PES utilizing the ωB97 functional [50] and 6-31G* basis set [51]. The initial velocities for the dynamics were the thermal ones corresponding to the randomly sampled temperatures in the range of 30-80 K. Each of the 30 trajectories for each isomer was run for~2.4 ps (time step = 1.21fs) unless otherwise stated. Twenty different starting structures were considered per system size (decamer and eicosamer), giving a total of 1200 trajectories.

Conclusions
We studied ionic polymerization processes on large clusters containing 10 and 20 acetylene molecules in conditions relevant to the ISM by means of AIMD. Our results reinforce the observations of previous studies on clusters consisting of 3-6 acetylene units, i.e., an increased number of spectator acetylenes promotes molecular growth. We observed a larger C 6 H 6 + /C 4 H 4 + ratio for decamers and eicosamers than what was previously found for acetylene clusters up to hexamers. Additionally, we observed C 8 H 8 + structures that were not obtained from smaller clusters (consisting of 3-6 units of acetylene). Eicosamer clusters exhibit results similar to those of the decamer clusters, demonstrating that the additional 10 acetylene molecules do not significantly influence the environment to promote further growth. Moreover, we predict that the formation of mono-and bi-cyclic structures with formulas (C 2 H 2 ) n + n = 2-4 produced by ionic polymerization processes can serve as a guide for astronomers in search of new molecules in the ISM.