The unimolecular dissociation and infrared radiative cooling rates of cationic 1-hydroxypyrene (OHPyr^{+}, C_{16}H_{10}O^{+}) and 1-bromopyrene (BrPyr^{+}, C_{16}H_{9}Br^{+}) are measured using a cryogenic electrostatic ion beam storage ring. A novel numerical approach is developed to analyze the time dependence of the dissociation rate and to determine the absolute scaling of the radiative cooling rate coefficient. The model results show that radiative cooling competes with dissociation below the critical total vibrational energies *E*_{c} = 5.39(1) eV for OHPyr^{+} and 5.90(1) eV for BrPyr^{+}. These critical energies and implications for radiative cooling dynamics are important for astrochemical models concerned with energy dissipation and molecular lifecycles. The methods presented extend the utility of storage ring experiments on astrophysically relevant ions.

## I. INTRODUCTION

Polycyclic Aromatic Hydrocarbons (PAHs) have long been thought to be ubiquitous in the Interstellar Medium (ISM), as evidenced by the so-called Aromatic Infrared Bands (AIBs) observed in emission at wavelengths coincident with vibrational transition energies of PAHs.^{1,2} Recently, radio astronomy was used to identify the first individual PAHs in space, including two isomers of cyanonaphthalene.^{3,4} A notable conclusion of these reports is that the observed abundances of the cyanonaphthalenes are six orders of magnitude higher than can be explained by astrochemical modeling.^{4} Furthermore, the observed species are much smaller than the PAHs previously assumed to be prevalent in space. These inconsistencies highlight the need for reliable experimental benchmarks of the stability and lifecycle of PAHs under interstellar conditions.

Quantitatively accurate rates of formation and destruction of PAHs are of considerable importance in interpreting astronomical observations and constructing astrochemical models, but laboratory data are scarce, particularly for substituted PAHs.^{5–8} The survival of PAHs in harsh interstellar environments depends on the interplay between photo- and collision-induced dissociation and radiative stabilization.^{9} Electrostatic Storage Devices (ESDs) allow highly excited PAH ions to be isolated in collision-free environments for timescales up to tens of seconds, enabling the determination of radiative cooling rates over several orders of magnitude in time.^{10–12} An important finding of such studies has been the significance of recurrent fluorescence in stabilizing hot PAH cations.^{13,14} One limitation of previous studies is that they generally rely on analytical fitting procedures that are often based on several approximations regarding the ions’ internal energy distributions.

We present a quantitative study of two substituted PAH cations, 1-hydroxypyrene (OHPyr^{+}, C_{16}H_{10}O^{+}) and 1-bromopyrene (BrPyr^{+}, C_{16}H_{9}Br^{+}), in which we have measured the spontaneous dissociation rate of hot ensembles of ions stored in a cryogenic ESD. We introduce a novel numerical approach for analyzing the time dependence of the dissociation rate, in which the rate is numerically transformed into the reciprocal space of dissociation rate coefficients. This framework provides an estimate of the initial temperature of the ion ensemble and sets bounds on the absolute scale of the radiative cooling rate coefficient.

^{7,15}the following dissociation channels are considered:

Hydroxy-substituted PAHs, among other species, have been shown to form in interstellar ice analogs irradiated by ultraviolet light and high-energy particles in experiments intended to simulate the photochemistry of cold molecular clouds.^{16–20} These photolysis products could be released into the ISM from dust grains heated by starlight at the edges of clouds.^{21} Intriguingly, neither hydroxy nor nitrile substituted PAHs are found to contribute significantly to the AIBs,^{1} even though the former are the dominant photoproducts of PAH photolysis in ice^{19,20} and the latter have been positively identified in cold clouds.^{4,22} Evidently, substituted PAHs formed in dark clouds do not survive in the more diffuse interstellar medium from which the AIBs emanate, and a more detailed understanding of their destruction routes is therefore needed.^{19} Bromo-PAHs have been investigated from a more fundamental perspective, as their “flat” transition states lend themselves to accurate rate calculations, providing robust benchmarking data.^{23}

## II. METHODS

### A. Experiments

Experiments were conducted at the DESIREE (Double ElectroStatic Ion Ring ExpEriment) infrastructure at Stockholm University.^{24} Cryogenic cooling of the DESIREE storage ring, which is schematically shown in Fig. 1, to ≈13 K results in a residual gas density on the order of $\u223c104$ cm^{−3}, consisting mostly of H_{2}.^{25} The excellent vacuum reduces the background rate due to collisions between stored ions and residual gas, increasing the dynamic range of the measurement by two to three orders of magnitude compared to a room-temperature ESD.^{26}

OHPyr [Tokyo Chemical Industry (TCI) $>98%$] and BrPyr (Sigma–Aldrich $>96%$) were sublimed from powder in a resistively heated oven coupled to an electron cyclotron resonance (ECR) ion source (Pantechnik Monogan). Helium was used as a support gas. Cations extracted from the source were accelerated to 56 keV. Mass-selected beams of cationic OHPyr^{+} (*m*/*z* = 218) and ^{81}BrPyr^{+} (*m*/*z* = 282) were stored in the DESIREE ion storage ring illustrated in Fig. 1. Beam oscillations due to ions injected along unstable trajectories were minimized using 12.5 mm diameter apertures placed before and after the lower straight section in Fig. 1.

After ion injection into the DESIREE storage ring, neutral fragments are formed from ions that retain significant internal energy from their formation in the ion source. Neutrals formed in the observation arm (upper straight section in Fig. 1) continue with high velocity toward a microchannel plate (MCP) detector comprised of custom ultra-high dynamic range MCPs (Photonis) that are suitable for high count rates at cryogenic temperatures.^{27} The measured neutral yield rate as a function of time *t* after the ions left the source is *R*(*t*).

### B. Data analysis and modeling

#### 1. Master equation simulations

*g*(

*E*,

*t*) of intact OHPyr

^{+}and BrPyr

^{+}was simulated using the master equation approach.

^{28}Further details of the methodology are presented elsewhere.

^{12}The vibrational level density

*ρ*(

*E*) is computed using the Beyer–Swinehart algorithm.

^{29}The vibrational frequencies

*ν*

_{s}and Einstein coefficients

*A*

_{s}for each vibrational mode

*s*are calculated at the B3LYP/6-31G(d,p) level of Density Functional Theory (DFT) as implemented in Gaussian 16.

^{30}The infrared radiative (vibrational) cooling rate is calculated within the simple harmonic cascade approximation

^{31}

*v*is the vibrational quantum number.

^{9}

^{,}

*A*

_{d}is the pre-exponential factor and

*E*

_{a}is the activation energy. We adopt the measured values of

*A*

_{d}= 2 × 10

^{14}s

^{−1}and

*E*

_{a}= 2.91 eV from the work of Lesniak

*et al.*for the loss of CO from OHPyr

^{+}. For BrPyr

^{+}, we use the values

*A*

_{d}= 1.5 × 10

^{15}s

^{−1}and

*E*

_{a}= 3.3 eV for Br-loss as for smaller bromo-PAHs.

^{15}The contribution to

*E*

_{a}from rotational excitation of the molecules is estimated to be small and comparable to the statistical uncertainties and is therefore assumed to be negligible.

^{32}

Recurrent fluorescence rate coefficients^{9} were calculated using the transition energies and oscillator strengths delivered by time-dependent DFT calculations. They were found to be insignificant compared to dissociation and infrared cooling and were not included in the simulations presented here.

*g*(

*E*,

*t*= 0) normalized such that

*N*(0) = ∫

*g*(

*E*,

*t*= 0)

*dE*= 1, the distribution is propagated according to the master equation

The first term gives the depletion of the population by unimolecular dissociation. The first term in brackets represents *v* + 1 → *v* IR photon emission from levels above *E*, while the second term represents *v* → *v* − 1 emission to levels below *E*. The time step *dt* is chosen to match the experimental data, with 32 extra points prior to the first experimental time bin to allow for the ∼100 *µ*s ion transit time from the source to the storage ring.

#### 2. Dissociation rates

*k*

_{d}varies rapidly with

*E*, the dissociation rate will not follow a simple exponential dependence, rather it is given by

*t*) is the per-particle decay rate, which is related to the experimental neutral yield rate

*R*(

*t*) by

*ϵ*

_{det}is the detector efficiency,

*L*

_{SS}= 0.95 m is the length of the straight section seen by the detector,

*C*= 8.7 m is the circumference of the storage ring, and

*N*

_{stored}is the number of ions in the stored ensemble measured using the Faraday cup (Fig. 1), averaged over the number of injection cycles.

To allow a quantitative comparison between the experimental results and those of our master equation simulations, we must put the measured neutral yield on an absolute scale, i.e., we must find *α* in Eq. (7). In a previous report,^{12} *ϵ*_{det} was determined by comparing the neutral yield at long times, after the spontaneous decay rate reaches the floor set by the residual gas collision rate, with the beam storage lifetime determined by measuring *N*_{stored} as a function of storage time. In the present experiments, however, the ion current at the end of the storage cycle was too low to be measured with the Faraday cup due to the need to limit saturation of the detector by high count rates. Here, we take a different approach, determining *α* directly by comparing the experimental neutral yield rate *R*(*t*) with the simulated Γ(*t*). This requires that we first fix the initial energy distribution for the simulation using the method detailed below.

^{33,34}decay rate curves measured in ESDs are fit with an analytical expression such as

*t*

^{−p}, where

*p*≈ 1, results from the broad energy distribution.

^{35}Measured deviations of

*p*from unity have been variously ascribed to the ion’s finite heat capacity,

^{35}competition between dissociation channels,

^{36}the shape of the internal energy distribution,

^{14,37}and blackbody infrared radiative dissociation.

^{26}The exponential factor is attributed to radiative cooling, which “quenches” the decay rate after a critical time $kc\u22121$, where the dissociation and cooling rate coefficients are comparable. The value of

*k*

_{c}is, unlike

*p*, an intrinsic property of the molecule or cluster ion under study. Note that the exponential quenching of the decay rate is

*not*due to loss of stored ion current due to residual gas collisions. The collision-limited lifetimes of molecular ion beams in DESIREE are measured in hundreds of seconds,

^{10,38–40}five orders of magnitude longer than typical values of $kc\u22121$. On the other hand, $kc\u22121$ is typically orders of magnitude longer than the revolution period of the ions around the ring. Thus, quenching of the decay rate is not likely to be influenced by the loss of ion current due to ions injected along unstable trajectories, switching transients, etc., or due to detector saturation effects, both of which are manifest mainly in the first few revolutions. However, when fitting experimental data to an equation such as Eq. (8), such experimental artifacts at early times can spuriously influence the determination of

*k*

_{c}. This motivates the development of an alternative analytical approach.

*p*= 1, Eq. (8) may be recognized as the Laplace transform of the Heaviside step function

*u*,

*t*) as the Laplace transform of some weighting function $F(k\u2032)$. If we neglect the redistributive effect of radiative cooling on the energy distribution, we can approximate $g(E,t)\u2248g\u2020(E)e\u2212ktott$, where

*k*

_{tot}=

*k*

_{d}+

*k*

_{IR}and

*g*

^{†}(

*E*) is the energy distribution of ions that contribute to Γ(

*t*), i.e., ions that are

*not*stabilized by radiative cooling. Changing variables with the substitution $dk\u2032dE=k\u2032ddElog(k\u2032)$, we have

*g*

^{†}(

*E*) by equating the integrands of Eq. (10) and identifying the dummy variable

*k*′ with the total decay rate

*k*

_{tot}.

*t*) for which inverse Laplace transforms exist, we numerically compute the Reimann sum,

*t*

_{i}are the experimental time bins and

*k*

_{j}are chosen to span the range of rate coefficients that can contribute to the experimental signal

*R*(

*t*), spaced evenly on a logarithmic scale, and Δ

*k*

_{j}=

*k*

_{j+1}−

*k*

_{j}. The number of

*k*

_{j}is chosen to be one less than the number of time points.

We estimate the solution to Eq. (11) by finding the pseudo-inverse $L+$ of the transform matrix $L$ using the singular value decomposition routine in the NumPy library.^{41} We truncate the decomposition at the largest number of singular values of $L$, which gives a physical (non-negative) resultant $F(kj)$.

*g*

^{†}(

*E*), we must find the energies

*E*

_{j}for which

*k*

_{tot}(

*E*

_{j}) =

*k*

_{j}. Recognizing that most ions which dissociate have vibrational energies for which

*k*

_{d}≫

*k*

_{IR}, we take

*k*

_{tot}≈

*k*

_{d}and find, from Eq. (10), the simple expression

*g*

^{†}(

*E*) distribution obtained from the experimental decay rate

*R*(

*t*) to the results of our master equation simulations for different initial temperatures. The simulations directly record

*g*

^{†}(

*E*) by comparing the full energy distribution

*g*(

*E*,

*t*) before and after each dissociation step. We choose the initial temperature that best reproduces the experimental result according to the overlap integral,

*t*) to the experimental

*R*(

*t*) determines the normalization factor

*α*.

## III. RESULTS AND DISCUSSION

From the weighting function $F(k\u2032)$ determined from Eq. (11) and the *k*_{d}(*E*) found in Eq. (4), we evaluated the effective energy distribution of dissociating ions contributing to the measured neutral yield *g*^{†}(*E*) from Eq. (12). The result for OHPyr^{+} is shown in Fig. 2. Also shown are results from our master equation simulations, where *g*^{†}(*E*) is tracked directly, for several initial temperatures. The distribution is bell-shaped, where the low-energy edge is set by the competition between dissociation and radiative cooling. The high energy edge reflects the warmest ions that survive transport from the source to the storage ring. The spectral weight near these edges is influenced by the initial energy distribution *g*(*E*, *t* = 0) and, thus, the initial temperature. The optimal temperature for reproducing the *g*^{†}(*E*) distribution derived from experiment is 1725 K for OHPyr^{+}.

Figure 3 shows the measured neutral particle detection rates *R*(*t*) for OHPyr^{+} and BrPyr^{+} on double-logarithmic axes. Also shown are the reconstructions of the experimental data according to Eq. (11). By comparing the simulated Γ(*t*) for the optimal initial temperatures to the experimental *R*(*t*), we find the constant of proportionality relating the two [Eq. (7)]. This normalization has been applied in Fig. 3 to align the simulated curves with those of the experiment. For both OHPyr^{+} and BrPyr^{+}, the agreement of the reconstruction with the experimental data is quite good at early times but diverges somewhat after the effect of radiative cooling becomes more pronounced. This is due to the truncation of the singular value decomposition of the transform matrix at rank 6 and 3 for OHPyr^{+} and BrPyr^{+}, respectively.

For BrPyr^{+}, our simulations did not reproduce the experimental results without the inclusion of a scaling factor *f*_{IR} modifying the radiative cooling rate coefficient *k*_{IR}(*E*) [Eq. (3)]. The procedure for determining *f*_{IR} is illustrated in Fig. 4. Recalling that our first approximation to the weighting function $F(k\u2032)$ was the step function *u*(*k*′ − *k*_{c}), which conventionally takes the value 0.5 at *k*′ = *k*_{c}, we find the value of *k*′ where our experimentally determined $F(k\u2032)$ first reaches half its maximum value (blue dot in Fig. 4). Next, we require that the scaled cooling rate coefficient crosses the dissociation rate coefficient at the critical energy *E*_{c}, where *f*_{IR}*k*_{IR}(*E*_{c}) = *k*_{d}(*E*_{c}) = *k*_{c} (orange dot in Fig. 4). For BrPyr^{+}, this results in values of *k*_{c} = 168(3) s^{−1} and *f*_{IR} = 1.74(3). Empirical scaling factors of this magnitude have been reported for other systems and result from inaccuracies in the calculated vibrational frequencies and the approximations inherent to the simple harmonic cascade model.^{38,39} The same procedure applied to OHPyr^{+} gives *k*_{c} = 140(4) s^{−1} and *f*_{IR} = 0.98(3). With the scaling factor *f*_{IR} = 1.74, the initial temperature that best reproduces our experimental results for BrPyr^{+} is 1850 K. This simulation is shown in Fig. 3. The simulated decay rates for both OHPyr^{+} and BrPyr^{+} agree well with the experimental data across three orders of magnitude in time. The critical vibrational energy at which dissociation and radiative cooling are competitive is *E*_{c} = 5.39(1) eV for OHPyr^{+} and 5.90(1) eV for BrPyr^{+}, which in each case is about twice the activation energy.

Similar values of *k*_{c} [1.31(2) × 10^{2} and 1.70(7) × 10^{2} s^{−1} for OHPyr^{+} and BrPyr^{+}, respectively] can be found by fitting the experimental decay rate to Eq. (8), as shown in Fig. 3. The agreement with the experimental data is comparable for the analytic fit [Eq. (8)] and the reconstruction [Eq. (11)]. The fit underestimates the decay rate at long times by about the same degree as the reconstruction overestimates it. We conclude that our numerical method based on the weighting function provides similar insight into the cooling process as the analytical approach. However, we have traded out the inscrutable power-law exponent *p* for the directly interpretable *g*^{†}(*E*).

^{+}and two modified rates accounting for different hypothetical experimental errors in Fig. 5. In scenario 1, the true rate

*R*

_{Expt.}has been modified to include the effect of machine-related beam losses due to ions injected into the storage ring along unstable orbits, switching transients, etc. This is modeled assuming half the injected beam is lost with 1/

*e*lifetimes equal to the revolution frequency

*f*

_{rev},

*τ*

_{d}= 1

*µ*s. The values of

*k*

_{c}determined from the weighting function and fitting approaches are given in Table I.

k_{c} (s^{−1})
. | $F(k\u2032)$ . | Fit . |
---|---|---|

Expt. | 140(4) | 131(5) |

Scen. 1 | 146(6) | 119(4) |

Scen. 2 | 130(4) | 146(5) |

k_{c} (s^{−1})
. | $F(k\u2032)$ . | Fit . |
---|---|---|

Expt. | 140(4) | 131(5) |

Scen. 1 | 146(6) | 119(4) |

Scen. 2 | 130(4) | 146(5) |

In both modified scenarios, the modified count rate differs significantly from the measured rate during the first few revolutions. These data points, however, have the highest statistical weight and leverage in a standard least-squares fit with variance weights. In Table I, it can be seen that the critical rate coefficient for radiative cooling is impacted by changing these early points in the modified scenarios. Compared to the experiment, the values of *k*_{c} determined from fitting differ by 1.9 and 2.1 standard deviations in scenarios 1 and 2, respectively. The differences in values determined from the weighting function are slightly smaller (at 0.8 and 1.8 standard deviations). Furthermore, in the expanded inset of Fig. 5, it can be seen that reconstruction of the rate from the weighting function (solid lines) follows all three sets of data points more closely than the fits (dotted lines) at early times. While far from exhaustive, analysis of these scenarios shows the weighting function method to be more robust against experimental artifacts than standard fitting methods.

## IV. CONCLUSIONS

We have measured the absolute unimolecular dissociation rates for two substituted PAH cations. We have applied a novel numerical approach to analyze the decay rates from which we determine the absolute scale of the radiative cooling rate coefficients. For both cations, we found the cooling rate to be consistent with a simple harmonic cascade model of infrared vibrational cooling, with no significant contribution from recurrent fluorescence.

The absence of significant cooling by recurrent fluorescence implies that OHPyr^{+} is less likely to survive in the interstellar medium than the fully aromatic PAH species that have previously been shown to be efficiently stabilized by this mechanism.

Compared to analytical modeling, our numerical approach does not require us to make assumptions about the cooling mechanism or approximations of the vibrational energy distribution. This method can be further improved by implementing a more rigorous solution to the inverse problem in Eq. (11) than the truncated singular value decomposition used here. The decay curve analysis procedure could be applied to laser-induced decay measurements to track the evolution of the vibrational energy distribution on longer timescales.^{42,43} It could also be applied in more complex situations where multiple species are present in the stored ion beam^{12} or where competing channels yield decay curves that are challenging to model analytically.^{33,37}

## ACKNOWLEDGMENTS

This work was supported by the Swedish Research Council (Grant Nos. 2016-03675 and 2020-03437), the Knut and Alice Wallenberg Foundation (Grant No. 2018.0028), the Olle Engkvist Foundation (Grant No. 200-575), and the Swedish Foundation for International Collaboration in Research and Higher Education (STINT, Grant No. PT2017-7328 awarded to J.N.B. and M.H.S.). We acknowledge the DESIREE infrastructure for provisioning of facilities and experimental support and thank the operators and technical staff for their invaluable assistance. The DESIREE infrastructure receives funding from the Swedish Research Council under Grant Nos. 2017-00621 and 2021-00155. This article is based upon the work from COST Action CA18212—Molecular Dynamics in the GAS phase (MD-GAS), supported by COST (European Cooperation in Science and Technology).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are openly available in Zenodo at http://doi.org/10.5281/zenodo.6512271.