This investigation implemented the DFT/B3LYP methodological framework. A tiered strategy guided basis set selection: initial geometry optimizations and transition state calculations were conducted by the cost-effective 6-31G(d) basis set, followed by 6-311G(d,p) for improving precision. The geometric optimizations and vibrational frequency analysis were ultimately performed using the 6-311++G(d,p) basis set. The incorporation of diffuse functions enabled superior treatment of electron correlation effects and weak intermolecular weak interactions. Method validation and computational reliability were verified through systematic comparison with existing literature data, geometry optimizations and frequency analyses were performed for n-C₄H₃, i-C₄H₃, C₂H₂, transition states, and products within the Gaussian^[16] program package using the 6-311++G(d,p) basis set, complemented by CBS-QB3 method for single-point energy calculations. The optimized reactants and products showed no imaginary frequencies, while identified transition states had one and only one imaginary frequency, which was consistent with the single imaginary frequency criterion of TST. To verify the transition states, Intrinsic Reaction Coordinate (IRC) calculations were performed on the transition states using the same basis set and computational method. The results confirmed reactants and products reside at local minima, transition states occupy first-order saddle points, and each properly connects corresponding reactants with products.

Transition State Theory provides the fundamental framework for describing chemical reaction rates in this study. Its central premise maintains that all chemical transformations proceed through a highest-energy transition state configuration during conversion from reactants to products. The primary mathematical description resides in the Eyring equation, which enables theoretical calculation of chemical reaction rate constants.

Based on the form involving the Gibbs free energy of activation ($ \Delta {G^\ddag } $):where k represents the reaction rate constant, $ {k_{\mathrm{b}}} $ is the Boltzmann constant (1.380649×10⁻²³ J/K), h is the Planck constant (6.62607015×10⁻³⁴ J·s) and T is the absolute temperature (K). The terms $ {Q^\ddag } $, $ {Q_{\mathrm{A}}} $ and $ {Q_{\mathrm{B}}} $ correspond to molecular partition functions for the transition state, reactant A and reactant B, respectively. $ \Delta {G^\ddag } $ signifies the standard Gibbs free energy of activation (J/mol) and R is the ideal gas constant (8.314 J/(mol·K)).

This study determined reaction rate constants using KiSThelP^[17] and Gaussian software based on Transition State Theory. All calculated rate constants refer to the high-pressure limit values. The pre-exponential factor A, temperature exponent n, and activation energy E_a were obtained by fitting the Arrhenius equation. Following frequency calculations at the B3LYP/6-311++G(d,p) level, energy values, vibrational frequencies, and molecular structure data were extracted from the Gaussian output files. Based on this information, KiSThelP was applied to obtain theoretical predictions of the benchmark gas-phase reaction rate constants. Generally, the empirical equation derived from the conventional TST method is adopted for this purpose:

In the equation above, $ \sigma $ represents the reaction path degeneracy, and N_A signifies Avogadro’s constant (which cancels out for unimolecular reaction rate constants with units of s⁻¹). $ {V^\ddag } $is the potential energy difference excluding zero-point vibrational energy. The terms ($ {Q^{{\mathrm{TS}}}} $) and ($ {Q^{\mathrm{R}}} $) denote the total partition functions for the transition state and the reactants, respectively. In rate constant calculations, the partition function incorporated the thermodynamic contributions from translational, rotational, and vibrational motions for all species (reactants and transition state). Consequently, the rate constants computed using the KiSThelP software inherently accounted for these effects, thus ensuring physical accuracy under high-temperature conditions. Comprehensive details concerning partition function treatment and rate constant derivation are available in reference^[18], and the complete computational methodology can be foundin the KiSThelP software manual.

In KiSThelP, one-dimensional quantum mechanical effects were characterized using the transmission coefficient, expressed as k^TST/T(T) = χ(T)×k^TST(T). This TST/T method introduced quantum tunneling corrections based on conventional TST. As an initial approximation, the Wigner correction method^[19] was applied, which relies solely on the imaginary frequency of the reaction mode to characterize the vibrational mode and quantifies its contribution as a multiplicative factor.where in $ {Im} {v^\ddag } $ denotes the magnitude (absolute value) of the imaginary frequency at the transition state.

To more accurately describe the relationship between the reaction rate constant k and temperature T and enhance predictive accuracy for reaction rates, the three-parameter Arrhenius equation^[20] was applied to calculate the rate constants over the temperature range of 300 to 4000 K. The equation is given as follows:where A is the pre-exponential factor, n refers to the temperature exponent, and E_a denotes the activation energy.

For simplified kinetic expressions, the two-parameter Arrhenius formula was also used for rate constant fitting, which is expressed as follows:

The addition reaction mechanism of n-C₄H₃ with C₂H₂ was explored. Based on Gaussian calculations, the reactive site of the n-C₄H₃ radical was determined through analysis of its vibrational modes. The radical active center at the carbon 1 (C1) atom was selected for calculating the addition reaction pathways with the acetylene molecule. Molecular structures of n-C₄H₃ and i-C₄H₃ are presented in Figure 1, where TS denotes the transition state, and IM represents the intermediate. For potential energy surface calculations, the total energy was determined by combining single-point energy (SPE) and the zero-point energy (ZPE) correction, which is expressed as: E(total) = E(SPE) + E(ZPE).

As shown in Figure 2, DFT calculations reveal the reaction mechanism for acetylene addition to the radical site of n-C₄H₃. The acetylene molecule initially overcame an energy barrier of 2.99 kcal/mol to attack the C1 radical site of n-C₄H₃, forming intermediate IM1 with an exothermicity of 36.73 kcal/mol. This calculated barrier compares favorably with the value of 4.85 kcal/mol reported by Lorries et al.^[21] using B3LYP/6-31G(d)//G3B3 methodology, with the 1.86 kcal/mol difference. This is attributed to our implementation of basis sets incorporating diffuse and polarization functions that more accurately described electron distribution, particularly in high electron-density regions like multiple bonds, and thus the precision of both geometric optimization and energy calculations was improved. Subsequently, the newly formed C−C single bond in IM1 underwent a conformational transformation through a rotational transition state, absorbing 2.29 kcal/mol of energy to form intermediate IM2. IM2 is then converted to IM3 via C−H bond torsional motion with a rotation angle of 81.8°. Finally, IM3 proceeded through three competing cyclization pathways: forming a cyclobutadiene derivative via the four-membered ring transition state TS6 (energy barrier 33.4 kcal/mol), generating a fulvene structure through the five-membered ring transition state TS5 (energy barrier 4.2 kcal/mol), and constructing a benzene ring skeleton via the six-membered ring transition state TS4 (energy barrier 4.38 kcal/mol). The IRC analysis confirmed that the six-membered ring formation proceeded through multiple intermediates (IM1→IM2→IM3→TS4→IM4) without evidence of a concerted transition state for direct six-membered ring formation.

Although a theoretical possibility exists for concerted [4+2] cycloaddition of acetylene to the n-C₄H₃ conjugated system directly generating a phenyl radical, the initial reaction of the n-C₄H₃ radical preferentially followed the stepwise radical addition mechanism, which represents the commonly accepted pathway in the literature. The energy barriers for IM3 to IM4, IM5 and IM6 calculated in this work show differences of 1.15, 1.48 and 0.06 kcal/mol, respectively, compared to the corresponding values of 5.53, 2.72 and 33.34 kcal/mol reported in the reference [21].

Geometrical optimization and frequency analysis of the i-C₄H₃ radical were performed using DFT with the B3LYP/6-311++G(d,p) functional and basis set. It is revealed that the C₂H₂ molecule preferentially underwent addition reactions with the carbon atoms at the C2 and C4 positions of the i-C₄H₃ radical.

As illustrated in Figure 3, different from the study by Lorries et al.^[21], this work computed the addition-cyclization reaction pathway involving the C4 site of i-C₄H₃ and acetylene. The C₂H₂ molecule and the i-C₄H₃ radical first overcame an energy barrier of 11.02 kcal/mol to form intermediate IM7 via addition. However, computational results from Wang et al.^[22] at the B2PLYP-D3/6-311++G(d,p), B2PLYP-D3/cc-pVTZ, B3LYP/6-311++G(d,p), B3LYP/cc-pVTZ, and B3LYP/cc-pV5Z levels reported a barrier of 14.59 kcal/mol for this step, which exceeded our value by 3.57 kcal/mol. Subsequently, the newly formed C–C single bond in IM7 underwent conformational rotation, generating intermediate IM8 via an endothermic process (ΔH = +2.41 kcal/mol). IM8 then underwent intramolecular C–H bond rotation to adjust its spatial configuration, forming intermediate IM9 with a suitable orientation for subsequent cyclization. IM9 proceeded through competitive pathways to form four-membered, five-membered, and six-membered ring products, with corresponding transition state barriers of 20.65 (four-membered ring, TS12), 16.70 (five-membered ring, TS10) and 33.14 kcal/mol (six-membered ring, TS13), respectively. A 1,3-hydrogen migration occurred within the five-membered ring intermediate IM10 (ΔG^‡ = 54.14 kcal/mol, TS11), which transferred a hydrogen atom to the meta-unsaturated carbon atom and yielded the fulvenyl intermediate IM11 with an extended conjugated system. Although this pathway produced a thermodynamically more stable product, its significantly elevated activation energy severely restricted the reaction. In contrast, the six-membered ring intermediate IM13 underwent a 1,2-hydrogen migration (ΔG^‡ = 23.32 kcal/mol, TS14) to form the aromatic phenyl-structured IM14 via ortho-hydrogen transfer.

As shown in Figure 4, an alternative reaction pathway involved the addition of the C₂H₂ molecule to the radical site at the C2 position of i-C₄H₃. This step needs to surmount a 12.55 kcal/mol barrier to form intermediate IM15. The calculated barrier reported in reference [21] is 8.49 kcal/mol, which results in a difference of 4.06 kcal/mol between the computational results. In the subsequent step, IM15 underwent conformational isomerization to form IM16 by surmounting a barrier of 2.14 kcal/mol, which was consistent with typical rotational barriers for C–C single bonds. IM16 then experienced an intramolecular conformational adjustment induced by C–H bond rotation, and yielded IM17 through a low-barrier structural rearrangement driven by non-bonding interactions. The final step involved a concerted ring-closure reaction of IM17, which overcame an activation energy of 14.92 kcal/mol to form IM18, with the cyclization barrier differing from reference [21] by 1.29 kcal/mol. Subsequently, the transformation from IM18 to IM19 followed an identical pathway to that from IM10 to IM11.

This study combined DFT with TST and the Wigner tunneling correction (TST/W) method to calculate reaction rate constants. All calculations were performed at the B3LYP/6-311++G(d,p) level of theory, and applied a frequency scaling factor of 0.998 (based on Kanchanakungwanhui et al., 2021, Version 5)^[23]. Using KiSThelP software, geometric structures and vibrational frequency data of reactants and transition states were extracted from Gaussian output files to determine rate constants for addition reactions, ring-closure steps, and rate-determining elementary reactions. Compared to the temperature range of 300−3000 K reported in reference [20], this work extended the computational range for reaction rate constants to 300−4000 K. This expansion covered kinetic processes under extreme combustion conditions like those in high-temperature gas turbine combustion chambers, and provided the first reported kinetic predictions for the ultra-high temperature regime (>3000 K), as illustrated in Figures 5−7.

As shown in Figure 5, kinetic calculations revealed strong temperature dependence for the n-C₄H₃ (C1) addition reaction with acetylene, and thus exhibited exponential rate enhancement with increasing temperature. The low energy barrier of this reaction, consistent with prior theoretical analysis, indicated facile reaction proceeding. The results showed a reaction rate constant of 3.08 × 10⁹ cm³·molecule^–1·s^–1 at 300 K.

When the temperature increased to 1000, 2000 and 3000 K, the reaction rate constants reached 3.65 × 10¹², 3.98 × 10¹³ and 1.26 × 10¹⁴ cm³·molecule^–1·s^–1, respectively, which demonstrated a significant temperature effect of spanning five orders of magnitude. This confirmed the substantial promoting effect of high-temperature combustion environments on the formation of soot precursors. This temperature sensitivity also aligned with the Arrhenius-type dependence predicted by transition state theory.

A notable site sensitivity was observed for the addition reaction of acetylene to i-C₄H₃, with the carbon 4 site having 1–2 orders of magnitude faster than that at the carbon 2 site. The results show rates of 5.69 × 10⁹ cm³·molecule^–1·s^–1 (C4) and 3.23 × 10⁸ cm³·molecule^–1·s^–1 (C2) at 1000 K, which increased to 1.64 × 10¹¹ and 6.62 × 10⁹ cm³·molecule^–1·s^–1, respectively, at 2500 K. The molecular structure of i-C₄H₃ illustrated its conformation: the carbon 4 site experienced less steric hindrance than the carbon 2 site, resulting in a more favorable stereochemical orientation for the addition transition state and a correspondingly lower activation energy barrier (approximately 1.53 kcal/mol), which preferentially directed reaction flux towards the more accessible C4 position. This finding was consistent with the Hammond Postulate prediction regarding the influence of spatial effects on reaction path selectivity. Molecular dynamics simulations further indicated that differences in configurational entropy could also influence the competition between reaction channels.

As shown in Figure 6, investigation of n-C₄H₃ (C1) addition to acetylene revealed distinct ring-size selectivity in subsequent cyclization: the formation rate of the five-membered ring (fulvenyl) was the fastest, the four-membered ring formation was the slowest, and the six-membered ring (phenyl) formation rate was intermediate, indicating a clear preference for five-membered ring formation. These results were consistent with the experimental findings from reference [15], where doping flames with diacetylene (C₄H₂, a precursor to n-C₄H₃) led to a significant increase in aromatic hydrocarbon yields. However, the primary products were not benzene-based compounds but five-membered ring fulvene derivatives, whose mass spectrometric signal intensity showed a linear positive correlation with the doping concentration. Kinetic analysis of the temperature dependence suggested that at 300 K, the reaction rates for the five- and six-membered ring pathways exceeded that of the four-membered ring pathway by 21 orders of magnitude. When the temperature increased to 2500 K, the rate ratio between the five- and six-membered ring pathways and the four-membered ring pathway decreased to 2 orders of magnitude. The calculated transition state energy barriers indicated higher reactivity for the five- and six-membered ring pathways than for the four-membered ring pathway, with an advantage of over 29 kcal/mol compared to the four-membered ring path, consistent with observed kinetic trends in ring-size dependent temperature effects. Furthermore, this kinetic preference for five- and six-membered rings was in agreement with the ring strain effects described by Baldwin’s Rules for ring closure^[24], where reduced ring strain corresponded to enhanced cyclization probability and lower-energy reaction pathways.

The selectivity for cyclic products in the reaction pathway involving the addition of acetylene to the C4 site of i-C₄H₃ exhibits substantial temperature dependence. Within the temperature range of 300–1800 K, the formation rate of the five-membered ring is higher than that of the four-membered ring, with the rate difference between them being approximately one order of magnitude.

However, at temperatures above 1800 K, the four-membered ring pathway becomes faster, owing to its higher temperature sensitivity (as reflected by a larger temperature exponent in the Arrhenius pre-exponential factor). Compared to the four- and five-membered rings, the six-membered ring transition state experiences significant steric hindrance, resulting in a formation rate nearly 12 orders of-magnitude slower at 300 K. As the temperature increases (1000–3000 K), the difference in rate constants gradually narrows to 1–2 orders of magnitude. This indicates that while the influence of energy barrier differences diminishes at high temperatures, it still dominates the pathway selectivity in the low-to-medium temperature combustion range. The rate differences among the various cyclization pathways could span several orders of magnitude, and temperature could even lead to a crossover in pathway dominance (for example, the intersection observed around 1800 K for the i-C₄H₃ pathways), all of which highlight the critical importance of kinetic analysis for predicting the actual product distribution in combustion systems. The influence of ring strain (kinetic control) and entropy effects on ring-size selectivity changed with temperature: lower temperatures suppressed the higher-barrier six-membered ring pathway, while higher temperatures reduced activation energy disparities and moderated selectivity. These findings were consistent with the trend predicted by Baldwin’s rules for ring closure regarding the effect of ring strain on transition state stability. For the reaction pathway involving acetylene addition to the C2 site, the rate constants at 300, 1000 and 3000 K were 3.40 × 10¹, 1.12 × 10⁹ and 1.58 × 10¹¹ s⁻¹, respectively. This represented a span of 10 orders of magnitude (10¹−10¹¹ s⁻¹), with the reaction rate increasing exponentially with temperature.

Based on DFT calculations, the cyclization reaction pathways of n-C₄H₃ and i-C₄H₃ with C₂H₂ were studied. Each pathway comprised multiple elementary reaction steps. According to the definition of the rate-determining step (RDS), the overall rate of a multi-step reaction is governed by the elementary step with the highest activation energy (i.e., the slowest step)^[25]. Theoretical calculations identified the respective RDS for each pathway as follows: Path A: IM1 → IM4, IM1 → IM5, IM1 → IM6; Path B: IM7 → IM11, IM7 → IM12, IM7 → IM13; Path C: IM18 → IM19. The calculation results are shown in Figure 7. In this study, due to the broad reaction temperature range (300 – 4000 K) and significant differences in energy barriers among the different reactions, the variation of reaction rates with temperature exhibited a complex trend. Furthermore, the rate constants showed a gradual plateauing with increasing temperature in the high-temperature regime. This phenomenon occurred because the influence of activation energy on the rate diminishes at elevated temperatures, while the influences of factors such as entropy change and tunneling effects become progressively more significant.

In Path A, the calculated kinetic parameters for the initial n-C₄H₃ + C₂H₂ addition revealed a substantial temperature exponent (n = 2.03). It was indicated that the energy of n-C₄H₃ was approximately 13.66 kcal/mol higher than that of i-C₄H₃, suggesting that n-C₄H₃ is a structurally more rigid and thermodynamically less stable radical isomer. Its initial addition reaction with acetylene (via TS1) was a bimolecular process requiring highly oriented collisions, resulting in a highly ordered transition state and consequently a significant entropy penalty (i.e., a large negative activation entropy). In the Arrhenius equation, the large positive n value compensates for this entropic effect. The total rate for forming five- and six-membered rings (fulvenyl and phenyl) from the reaction of n-C₄H₃ (C1) with acetylene exceeded that for four-membered ring formation by approximately 14–15 orders of magnitude or more. This difference exhibits a distinct temperature dependence: at low temperature (300 K), the rate difference is 14–15 orders of magnitude, narrowing to 1–2 orders of magnitude in the high-temperature range (3000–4000 K).

On reaction Path B, the reaction between i-C₄H₃ (C4) and acetylene exhibited distinct ring-size selective kinetic behavior. DFT calculations revealed that the formation of the six-membered ring product requires overcoming a high activation energy due to significant steric hindrance in the transition state, arising from twisting of the conjugated carbon chain.

Consequently, at a low temperature of 300 K, the six-membered ring channel proceeded at a negligible rate (k = 3.12 × 10⁻⁴⁵ s⁻¹), which lags behind the four-membered ring (k = 1.02 × 10⁻⁴ s⁻¹) and the five-membered ring (k = 2.26 × 10⁻¹⁰ s⁻¹) by 41 and 35 orders of magnitude, respectively. However, at a high temperature of 1500 K, the rate difference between the six-membered ring and the four-/five-membered rings narrowed to 6–8 orders of magnitude, which indicates a gradual reduction in disparity as temperature increases. In the reaction of i-C₄H₃ (C2) with acetylene, the hydrogen migration step leading to five-membered ring (fulvenyl) formation was highly sensitive to the energy barrier. Transition state analysis identified a substantial barrier of 54.14 kcal/mol for this process, resulting in a minimal rate constant (k = 1.03 × 10⁻²⁸ s⁻¹). The rate increased exponentially through the medium-temperature range (1000–2000 K), reaching 1.41 × 10¹ s⁻¹ at 1000 K and 3.32 × 10⁷ s⁻¹ at 2000 K. In the high-temperature dominant regime (3000−4000 K), the rate constant surpassed the 10⁹–10¹⁰ s⁻¹ magnitude.

Regarding the experimental evidence supporting the computational results, this study primarily focused on elucidating the reaction mechanisms at the atomic and molecular level through quantum chemistry calculations and transition state theory. The fitted kinetic parameters are listed in Table 1 (three-parameter Arrhenius equation) and Table 2 (two-parameter Arrhenius equation).

Although the present work is a theoretical investigation without direct experimental data, the main conclusions are consistent with experimental observations reported in the literature. Specifically, the finding that the n-C₄H₃ pathway predominantly led to five-membered ring (fulvene) formation aligns with experimental results from flame-doping studies by McEnally et al.^[15]. Additionally, the premise that both i-C₄H₃ and n-C₄H₃ isomers coexisted is directly supported by the experimental discovery of Yang et al.^[14] using synchrotron radiation vacuum ultraviolet photoionization technology.