# Cascade infrared thermal photon emission

Klavs Hansen  Lanzhou Center for Theoretical Physics, Key Laboratory of
Theoretical Physics of Gansu Province, Lanzhou University, Lanzhou, Gansu
730000, China  Center for Joint Quantum Studies and Department of Physics,
School of Science, Tianjin University, 92 Weijin Road, Tianjin 300072, China
Ori Licht  Department of Physics and Insitute for Nanotechnology and Advanced
materials, Bar Ilan University, Ramat-Gan 5290002, Israel  Adeliya Kurbanov
Department of Physics and Insitute for Nanotechnology and Advanced materials,
Bar Ilan University, Ramat-Gan 5290002, Israel  Yoni Toker  Department of
Physics and Insitute for Nanotechnology and Advanced materials, Bar Ilan
University, Ramat-Gan 5290002, Israel

###### Abstract

The time development of the excitation energy of molecules and clusters
cooling by emission of thermal vibrational infrared radiation has been
studied. The energy distributions and the photon emission rates develop into
near-universal functions that can be characterized with only a few parameters,
irrespective of the precise vibrational spectra and oscillator strengths of
the systems. The photon emission constant and emitted power averaged over all
thermally populated states vary linearly with total excitation energy with a
small offset. The time developments of ensemble internal energy distributions
are calculated with respect to their first two moments. For the derived linear
dependence of the emission rate constant, these results are exact.

##  1 Introduction

Emission of photons is an important cooling channel for molecules and clusters
in the absence of collisions. The decay channel is becoming increasingly
interesting in astrophysical context due to the growing number of molecules
observed in interstellar space 1, 2. The vastly different frequency factors
and activation energies associated with the different thermal processes
(fragmentation, thermionic emission, recurrent fluorescence, and infrared (IR)
emission) create a hierarchy of time scales on which the possible thermal
decay channels of thermally excited molecules proceed. The IR photon emission
is located at the lowest energies of this hierarchy, making it the energy-loss
channel that will dominate at the latest stage of the cooling process. The
development of cryogenic storage rings has opened the possibility to study
cooling processes on very long time scales in the laboratory 3, 4, 5, 6, 7, 8.
The extension of feasible storage times to thousands of seconds in these
devices places a special emphasis on the understanding of this long time
emission of thermal infrared photons from stored ions 9, 4, 6, 10, 7, 8, 11,
12

Before embarking on such a calculation of rate constants and energy
distributions, it is worthwhile to consider two aspects of the thermal
description of the vibrational degrees of freedom of a molecule or cluster.
One is the question of ergodicity. The low internal energies of the emitting
molecules/clusters raise this question with more force than for other decay
channels that need higher excitation energies to be sustained. Harmonic
oscillators do not couple, and when a collection of rigorously uncoupled
harmonic oscillators is used to describe the thermal properties of the
vibrational degrees of freedom, each of the modes will cool down
independently; clearly, for a system of this type the cooling proceeds
independently for all oscillators. The time dependence and average values can
be calculated for such a system without any further ado.

However, as molecular vibrations are not harmonic, energy can flow between
different vibrational modes, a process know as intramolecular vibrational
redistribution (IVR) 13, 14, 15. IVR is indeed an essential component of our
understanding of many phenomena at the heart of physical chemistry, such as
reaction rate theories, fluorescence and coherent control.

The time scales associated with IVR are determined by the degree of
anharmonicity, and scale with the density of states in the simplest large
molecule case (see ref. 16 for an introduction to the subject). These IVR time
scales are short enough to make IVR efficient on the timescales of
milliseconds and longer, which are the characteristic times of IR radiation.
Radiative cooling measurements in storage rings do indeed indicate that IVR is
the correct framework to describe the radiation at long times 6, 7, 12. The
’rapid exchange limit’ denotes the case in which IVR is much faster than IR
emission, and consequently, all states with the same total internal energy
will mix rapidly. We will use the rapid exchange limit Ansatz here.

Another point worth mentioning is the little considered fact that vibrational
motion carries angular momentum. The simplest case of linear molecules has
been analyzed in 17. It is also well-known from the description of the
excitation modes of helium droplets 18, 19, 20. The conservation of this
quantity imposes a constraint on the states that mix in the IVR process. This
clearly has consequences also for the photon emission rates. However, a
general theory is not yet available for molecules, and an in-depth treatment
of the subject is beyond the scope of this article. We will restrict the work
here to consider only the case of unrestricted IVR.

Within this framework we will show that vibrational cooling produces an almost
linear scaling of the cooling rate constants with the internal energy of the
system and, as a consequence, that the photon emission rate constant decays
exponentially with time and with a rate constant that can be computed. We will
use this scaling to calculate the time development of the general case where
energy distributions are not delta functions.

##  2 The harmonic cascade model

We describe the $s=3N-6$ nuclear vibrations of a molecule with $N$ atoms as
harmonic oscillators with energies $\epsilon_{i}$, $i=1,...,s$. In this
approximation, which is sufficiently accurate at the low energies we consider,
the total excitation energy, $\varepsilon$, of a molecule is therefore a sum
over the mode energies:

$\varepsilon=\sum_{i}\nu_{i}\epsilon_{i}$ (1)

where $\nu_{i}$ are the number of excitation quanta in mode $i$. A list of
these occupation numbers, $\\{\nu_{i}\\}$ specifies the system and will be
denoted its state. Within the harmonic approximation, IR radiative emission
only occurs between neighboring levels of each oscillator, i.e,
$\nu_{i}\rightarrow\nu_{i}-1$. We will denote the rate constant for this
transition,
$k_{\nu_{i}\rightarrow\nu_{i}-1}^{i}=1/\tau_{\nu_{i}\rightarrow\nu_{i}-1}^{i}$,
by the term emission rate constant of the state $\nu_{i}$ of oscillator $i$.
The matrix element that determines this rate constant for a particular
vibrational mode, $i$, is the dipole matrix element:

$\langle\nu_{i}-1|x|\nu_{i}\rangle_{i}\propto\sqrt{\nu_{i}}$ (2)

The rate constant for this particular mode is then given by:

$k_{\nu_{i}\rightarrow\nu_{i}-1}^{i}=\nu_{i}k^{i}_{1\rightarrow 0}.$ (3)

In the case of a single vibrational mode, with quantum energy $\epsilon_{1}$,
the photon emission rate constant is therefore proportional to energy;
$k(\varepsilon)=\left(\varepsilon/\epsilon_{1}\right)k_{1\rightarrow 0}^{1}$.

###  2.1 A Two Mode Example

Before treating the general case it is instructive to examine the simple case
of a system with only two vibrational modes in which
$\epsilon_{2}=2\epsilon_{1}$ and only one of the modes is IR active. We set
$k_{\rightarrow 0}^{2}=0$. In this case, illustrated in Fig. 1, the number of
states of the system where $\varepsilon=2m\epsilon_{1}$ or
$\varepsilon=2m\epsilon_{1}+1$ are both equal equal to $m+1$, where $m$ is a
non-negative integer. For $\varepsilon=2\epsilon_{1}$ the system is half the
time in the state $\\{2,0\\}$ which decays at a rate of $k_{2\rightarrow
1}^{1}=2k_{1\rightarrow 0}^{1}$, and half the time in the state $\\{0,1\\}$
which decays at the rate of $k_{1\rightarrow 0}^{2}=0$, and therefore:
$k(2\epsilon_{1})=k_{1\rightarrow 0}^{1}$. It is not hard to show that
$k(2m\epsilon_{1})$ and $k(2m\epsilon_{1}+1)$ are both equal to
$mk_{1\rightarrow 0}^{1}$. Thus, once again we see a general linear dependence
of the cooling rate on internal energy, albeit with fluctuations around the
general trend.

Figure 1: An illustration of the harmonic cascade model for the case of a
system with two modes for which $\epsilon_{2}=2\epsilon_{1}$, and only the low
energy mode optically active. The allowed quantum states, defined by
$\\{\nu_{1},\nu_{2}\\}$, are represented by points on a two dimensional
lattice. States with the same total energy
$\varepsilon=\nu_{1}\epsilon_{1}+\nu_{2}\epsilon_{2}$ lie on diagonals, as
shown. IVR allows for the system to freely transit between states of the same
energy. The yellow arrows indicate the transitions that can occur via
radiative cooling.

###  2.2 General Formulation

In more realistic situations where the different vibrational energies are not
exact multiples of each other, the rate constant for emission from the
thermally populated states of oscillator $i$ is still calculated by summing
over the products of populations and the energy-specified rate constants. The
populations are calculated with the level densities (densities of states). Let
$\rho(\varepsilon,\nu_{i}=m)$ denote the density of states of total energy
$\varepsilon$ for which the vibrational mode $i$ is occupied $m$ times. The
probability for the vibrational mode to be occupied exactly $m$ times is then
$\rho(\varepsilon,\nu_{i}=m)/\rho(\varepsilon)$. Hence the rate of decay
through mode $i$ is given by:

$k^{i}(\varepsilon)=\frac{1}{\rho(\varepsilon)}\sum_{m=1}^{\varepsilon/\epsilon_{i}}\rho(\varepsilon,\nu_{i}=m)k_{m\rightarrow
m-1}^{i}=\frac{1}{\rho(\varepsilon)}\sum_{m=1}^{\varepsilon/\epsilon_{i}}m\rho(\varepsilon,\nu_{i}=m)k_{1\rightarrow
0}^{i},$ (4)

which can be simplified using:

$\sum_{m=1}^{\varepsilon/\epsilon_{i}}m\rho(\varepsilon,\nu_{i}=m)=\sum_{m=1}^{\varepsilon/\epsilon_{i}}\rho(\varepsilon,\nu_{i}\geq
m)=\sum_{m=1}^{\varepsilon/\epsilon_{i}}\rho(\varepsilon-m\epsilon_{i}).$ (5)

Here $\rho(\nu_{i}\geq m,\varepsilon)$ denotes the number of states of total
energy $\varepsilon$ where mode $i$ is occupied at least $m$ times. To prove
the first equality we note that every state appears once in the sum on the
left multiplied by $m$ while it is counted $m$ times in the sum on the right
of this equality sign. For the equation on the right we note that once the
mode $i$ is occupied at least $m$ times there remains an amount of
$\varepsilon-m\epsilon_{i}$ energy to distribute within the system. Using Eq.
4 with Eq. 5 leads to:

$k^{i}(\varepsilon)=\frac{1}{\rho(\varepsilon)}\sum_{m=1}^{[\varepsilon/\epsilon_{i}]}\rho(\varepsilon-m\epsilon_{i})k_{1\rightarrow
0}^{i}=\frac{1}{\rho(\varepsilon)}\sum_{m=0}^{[\varepsilon/\epsilon_{i}]-1}\rho((m+\Delta
m)\epsilon_{i})k_{1\rightarrow 0}^{i}.$ (6)

The precise upper limit of the summation is indicated here with by taking the
integer part of the energy in units of $\varepsilon_{i}$,
$[\varepsilon/\varepsilon_{i}]$. The number $\Delta m$ is the fractional part
of $\varepsilon/\epsilon_{i}$, $\Delta
m\equiv\varepsilon/\varepsilon_{i}-[\varepsilon/\varepsilon_{i}]$. The total
cooling rate constant of the molecule is the sum of such contributions from
each vibration:

$\displaystyle k(\varepsilon)=\sum_{i}k^{i}(\varepsilon)$ (7)

The use of Eqs. 6 and Eq. 7 along with the density of states is known as the
’harmonic cascade model. 21, 6

An example of the results of such a calculation is shown in Fig. 2. As with
the single oscillator and the two mode examples discussed above, we observe
that a general linear dependence of $k$ on $\varepsilon$, as well as
fluctuations on top of the general trend. The figure also shows the effect of
a variation of the width of the energy bin, $\delta\varepsilon$, used for the
calculation of the density of states. The curves for different values of
$\delta\varepsilon$ agree, up to the fluctuations on top of the general trend.
We need therefore not consider this issue any further.

Figure 2:  An example of the calculation of $k(\varepsilon)$ for the case of
Al${}_{4}^{-}$, based on the vibrational frequencies and Einstein coefficients
calculated in Ref.12.

##  3 Rate of Decay

In the special case where all vibrational frequencies in a molecule are
identical (the Einstein molecule), simple combinatorics can be used to show
that $k(\varepsilon)=\sum_{i}k_{1\rightarrow 0}^{i}\varepsilon/s\epsilon_{\rm
E}=\langle k_{1\rightarrow 0}\rangle\varepsilon/\epsilon_{\rm E}$, with
$\epsilon_{\rm E}$ the common vibrational quantum energy. Trianglar brackets
indicate averaging over the different vibrational modes, i.e.:

$\langle x\rangle\equiv\frac{1}{s}\sum_{i}x_{i}.$ (8)

To treat the general case in which the vibrational levels are not degenerate,
we apply the criterion of agreement of the level density with the known high
energy limit. Inserting the asymptotic limit for the density of states,

$\displaystyle\rho(\varepsilon)=\frac{\left(\varepsilon+\frac{1}{2}\sum_{j=1}^{s}\epsilon_{j}\right)^{s-1}}{(s-1)!\prod_{j=1}^{s}\epsilon_{j}}$
(9)

into Eq. 6 gives

$\displaystyle\frac{k^{i}(\varepsilon)}{k^{i}_{1\rightarrow 0}}$
$\displaystyle=$
$\displaystyle\frac{\sum_{m=0}^{\varepsilon/\epsilon_{i}-1}\left(m\epsilon_{i}+\frac{1}{2}\sum_{i}\epsilon_{i}\right)^{s-1}}{\left(\varepsilon+\frac{1}{2}\sum_{j}\epsilon_{j}\right)^{s-1}}=\frac{\sum_{m=0}^{N_{i}-1}(m+V_{i})^{s-1}}{(N_{i}+V_{i})^{s-1}}$
$\displaystyle=$
$\displaystyle\frac{1}{(N_{i}+V_{i})^{s-1}}\left(\sum_{m=0}^{N_{i}+V-1}m^{s-1}-\sum_{m=0}^{V-1}m^{s-1}\right),$

where $N_{i}\equiv\varepsilon/\epsilon_{i}$ and
$V_{i}\equiv\sum_{j}\epsilon_{j}/2\epsilon_{i}$. The $\sum_{m=0}^{V-1}m^{s-1}$
term is asymptotically negligible compared with the first and can therefore be
ignored. The sum can be approximated using Faulhaber’s formula (which are also
the first two terms the Euler-Maclaurin formula):

$\sum_{m=0}^{N_{i}+V_{i}-1}m^{s-1}\simeq\frac{(N_{i}+V_{i}-1)^{s}}{s}+\frac{(N_{i}+V_{i}-1)^{s-1}}{2}\simeq\frac{(N_{i}+V_{i})^{s}}{s}-\frac{(N_{i}+V_{i})^{s-1}}{2}+O\left((N_{i}+V_{i})^{s-2}\right)$
(11)

Here the rightmost expression was derived using a Taylor expansion of the two
powers. Inserting this approximation into Eq. 3 leads to:

$\frac{k^{i}(\varepsilon)}{k^{i}_{1\rightarrow
0}}\simeq\frac{N_{i}+V_{i}}{s}-\frac{1}{2}=\frac{1}{s\epsilon_{i}}\varepsilon-\frac{1}{2s\epsilon_{i}}\sum_{j}\epsilon_{j}-\frac{1}{2}$
(12)

which can be written as:

$k^{i}(\varepsilon)=\frac{k^{i}_{1\rightarrow
0}}{s\epsilon_{i}}\varepsilon-\frac{k^{i}_{1\rightarrow
0}}{2}\left(\frac{\langle\varepsilon\rangle}{\epsilon_{i}}-1\right)\equiv
a_{i}\varepsilon+b_{i}$ (13) $a_{i}=\frac{k_{1\rightarrow
0}^{i}}{\epsilon_{i}s}~{}~{}~{}~{}b_{i}=\frac{k^{i}_{1\rightarrow
0}}{2}\left(1-\frac{\langle\epsilon\rangle}{\epsilon_{i}}\right)$ (14)

An interesting consequence is that:

$\sum_{i}\epsilon_{i}k^{i}(\varepsilon)=\frac{\sum k^{i}_{1\rightarrow
0}}{s}\varepsilon=\langle k_{1\rightarrow 0}\rangle\varepsilon$ (15)

The total decay rate $k(\varepsilon)$ is therefore given by:

$k(\varepsilon)=\sum_{i}k^{i}(\varepsilon)\equiv a\varepsilon+b$ (16)

with:

$a=\sum_{i}a_{i}=\left\langle\frac{k_{1\rightarrow
0}}{\epsilon}\right\rangle,~{}~{}~{}b=\sum_{i}b_{i}=\frac{s}{2}\left(\langle
k_{1\rightarrow 0}\rangle-\langle\frac{k_{1\rightarrow
0}}{\epsilon}\rangle\langle\epsilon\rangle\right)$ (17)

In particular we find that, in the asymptotic limit, $k$ has a linear
dependence on internal energy $\varepsilon$. A similar scaling has previously
been given in Ref. 22, albeit without the offset in the linear dependence of
the rate constant that will be derived here. Figure 3 shows a comparison of
the calculation of the decay rate $k(\varepsilon)$ using the harmonic cascade
model calculated with the Beyer-Swinehart level density, compared with the
results of Eq. 16. Notably, while Eq. 16 does not account for the fluctuations
on top of the general trend, it does seem to capture the general trend
extremely well.

Figure 3:  Calculation of the IR cooling rate $k(\varepsilon)$ using the
harmonic cascade model (Blue), for a few representative cases. For the two-
level model (Upper left panel) we calculated two cases, where either the first
or the second vibration is IR active and the other is not. For the cases of
Perylene 8 (upper right frame), carbon clusters 6(lower left frame) and
Al${}_{4}^{-}$ (lower right frame), we use the vibrational frequencies and
Einstein coefficients calculated in the respective works. The red line in each
panel is the theoretical $k(\varepsilon)$ calculated according to Eq. 7.

##  4 Population Evolution

Figure 4:  A: calculation of the evolution $g(\varepsilon)$ according to the
harmonic cascade model, for the case of Al${}_{4}^{-}$, starting with an
initial Gaussian distribution with a width of $0.2~{}$eV. B: mean energy,
$\,\overline{\\!{\varepsilon}}$ as a function of time, the red line
corresponds to Eq. 22, C: the variance as a function of time, compared with
Eq. 4.

Let $g(\varepsilon)$ denote the normalized energy distribution of a system,

$\displaystyle\int_{0}^{\infty}g(\varepsilon)d\varepsilon=1.$ (18)

The harmonic cascade model allows a calculation of the evolution of
$g(\varepsilon)$ with time according to:

$\frac{{\rm d}g(\varepsilon)}{{\rm
d}t}=\sum_{i}k^{i}(\varepsilon+\epsilon_{i})g(\varepsilon+\epsilon_{i})-k(\varepsilon)g(\varepsilon),$
(19)

where the sum runs over the vibrational modes. The first term comes from
states with higher energy radiatively decaying into the states of interest,
and $k(\varepsilon)\equiv\sum_{i}k^{i}(\varepsilon)$ is the total rate
constant for decay out of the system with energy $\varepsilon$. In this
section we will find the time evolution of $g(\varepsilon)$ in terms of its
mean,
$\,\overline{\\!{\varepsilon}}=\int_{0}^{\infty}g(\varepsilon)\varepsilon
d\varepsilon$, and its variance:
$\sigma^{2}=\int_{0}^{\infty}g(\varepsilon)(\varepsilon-\,\overline{\\!{\varepsilon}})^{2}d\varepsilon$.
We will analyze the high energy limit, with the understanding that this limit
is a fairly good approximation to reasonably low energies, up to the
fluctuation induced by the quantum discreteness (see fig. 3). In this limit
the total emission rate constant is given by Eq. 16. With that expression, the
energy develops as:

$\displaystyle\frac{{\rm d}\varepsilon}{{\rm
d}t}=-\sum_{i}k^{i}(\varepsilon)\epsilon_{i}=-\sum_{i}\epsilon_{i}k^{i}_{1\rightarrow
0}\left(\frac{\varepsilon}{s\epsilon_{i}}+\frac{\sum_{j}\epsilon_{j}}{2s\epsilon_{i}}-\frac{1}{2}\right)=-\frac{\varepsilon}{s}\sum_{i}k^{i}_{1\rightarrow
0}.$ (20)

Note that upon summation the offset terms cancel out
($\sum\epsilon_{i}\left(\frac{\sum_{j}\epsilon_{j}}{2s\epsilon_{i}}-\frac{1}{2}\right)=0$).
The energy therefore decreases exponentially:

$\displaystyle\varepsilon(t)=\varepsilon(0){{\rm
e}}^{-\sum_{i}k^{i}_{1\rightarrow 0}t/s}=\varepsilon(0){{\rm e}}^{-\langle
k_{1\rightarrow 0}\rangle t}.$ (21)

The exponent in this expression is the same for all energies. Averaging over a
distribution will therefore give the same decrease for the mean energy:

$\displaystyle\,\overline{\\!{\varepsilon}}=\,\overline{\\!{\varepsilon}}(t=0){{\rm
e}}^{-\langle k_{1\rightarrow 0}\rangle t}.$ (22)

The time evolution of the variance of $g(\varepsilon)$, $\sigma^{2}$, is
governed by the equation

$\displaystyle\frac{{\rm d}\sigma^{2}}{{\rm d}t}=\frac{{\rm d}}{{\rm
d}t}\int_{0}^{\infty}(\varepsilon-\,\overline{\\!{\varepsilon}})^{2}g(\varepsilon){\rm
d}\varepsilon=\int_{0}^{\infty}\left[\frac{{\rm d}g}{{\rm
d}t}\times\left(\varepsilon-\,\overline{\\!{\varepsilon}}\right)^{2}-2g(\varepsilon)\times(\varepsilon-\,\overline{\\!{\varepsilon}})\frac{{\rm
d}\,\overline{\\!{\varepsilon}}}{{\rm d}t}\right]{\rm d}\varepsilon.$ (23)

The last term is zero by definition of the mean value. The first term requires
the time evolution of $g$, according to Eq.19. This gives:

$\displaystyle\frac{{\rm d}\sigma^{2}}{{\rm
d}t}=\int_{0}^{\infty}\left[\left(\varepsilon-\,\overline{\\!{\varepsilon}}\right)^{2}\left(\sum_{i}k^{i}(\varepsilon+\epsilon_{i})g(\varepsilon+\epsilon_{i})-k^{i}(\varepsilon)g(\varepsilon)\right)\right]{\rm
d}\varepsilon,$ (24)

or, exchanging integration and summation,

$\displaystyle\frac{{\rm d}\sigma^{2}}{{\rm
d}t}=\sum_{i}\int_{0}^{\infty}\left(\varepsilon-\,\overline{\\!{\varepsilon}}\right)^{2}\left[k^{i}(\varepsilon+\epsilon_{i})g(\varepsilon+\epsilon_{i})-k^{i}(\varepsilon)g(\varepsilon)\right]{\rm
d}\varepsilon.$ (25)

The rate constants are represented with the two (known) parameters
$a_{i},b_{i}$ as

$\displaystyle k_{i}=a_{i}\varepsilon+b_{i}.$ (26)

Each term in the sum can then be written as

$\displaystyle\frac{{\rm d}\sigma_{i}^{2}}{{\rm d}t}$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}\left(\varepsilon-\,\overline{\\!{\varepsilon}}\right)^{2}\left[\left(a_{i}(\varepsilon+\epsilon_{i})+b_{i}\right)g(\varepsilon+\epsilon_{i})-\left(a_{i}(\varepsilon)+b_{i}\right)g(\varepsilon)\right]{\rm
d}\varepsilon=$ $\displaystyle-$
$\displaystyle\int_{0}^{\infty}\left(\varepsilon-\,\overline{\\!{\varepsilon}}\right)^{2}\left(a_{i}\varepsilon+b_{i}\right)g(\varepsilon){\rm
d}\varepsilon$ $\displaystyle+$
$\displaystyle\int_{0}^{\infty}\left(\varepsilon-\,\overline{\\!{\varepsilon}}+\epsilon_{i}\right)^{2}\left(a_{i}(\varepsilon+\epsilon_{i})+b_{i}\right)g(\varepsilon+\epsilon_{i}){\rm
d}\varepsilon$ $\displaystyle+$
$\displaystyle\int_{0}^{\infty}\epsilon_{i}^{2}\left(a_{i}\varepsilon+b_{i}\right)g(\varepsilon){\rm
d}\varepsilon$ $\displaystyle-$ $\displaystyle
2\int_{0}^{\infty}\epsilon_{i}\left(\varepsilon+\epsilon_{i}-\,\overline{\\!{\varepsilon}}\right)\left(a_{i}(\varepsilon+\epsilon_{i})+b_{i}\right)g(\varepsilon+\epsilon_{i}){\rm
d}\varepsilon.$

The first and the second terms on the right hand side cancel. Collecting terms
from the remaining two integrals gives

$\displaystyle\frac{{\rm d}\sigma_{i}^{2}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\epsilon_{i}^{2}\left(a_{i}\varepsilon+b_{i}\right)g(\varepsilon){\rm
d}\varepsilon-2\int_{0}^{\infty}\epsilon_{i}\left(\varepsilon+\epsilon_{i}-\,\overline{\\!{\varepsilon}}\right)\left(a_{i}(\varepsilon+\epsilon_{i})+b_{i}\right)g(\varepsilon+\epsilon_{i}){\rm
d}\varepsilon$ $\displaystyle=$ $\displaystyle
a_{i}\epsilon_{i}^{2}\,\overline{\\!{\varepsilon}}+b_{i}\epsilon_{i}^{2}-2a_{i}\epsilon_{i}\sigma^{2}.$

Inserting the values of $a_{i}$ and $b_{i}$ gives

$\displaystyle\frac{{\rm d}\sigma_{i}^{2}}{{\rm
d}t}=a_{i}\epsilon_{i}\left(\epsilon_{i}\,\overline{\\!{\varepsilon}}-2\sigma^{2}\right)+b_{i}\epsilon_{i}^{2}=\frac{k^{i}_{1\rightarrow
0}}{s}\left(\epsilon_{i}\,\overline{\\!{\varepsilon}}-2\sigma^{2}\right)+\frac{1}{2}k^{i}_{1\rightarrow
0}\varepsilon_{i}^{2}-\frac{k^{i}_{1\rightarrow
0}}{2}\epsilon_{i}\langle\epsilon\rangle.$ (29)

Summation over all modes gives

$\displaystyle\frac{{\rm d}\sigma^{2}}{{\rm
d}t}=\frac{\,\overline{\\!{\varepsilon}}}{s}\sum_{i}k^{i}_{1\rightarrow
0}\epsilon_{i}-\sigma^{2}\frac{2}{s}\sum_{i}k^{i}_{1\rightarrow
0}+\frac{1}{2}\sum_{i}k^{i}_{1\rightarrow
0}\epsilon_{i}^{2}-\sum_{i}\frac{k^{i}_{1\rightarrow
0}}{2}\epsilon_{i}\langle\epsilon\rangle,$ (30)

or

$\displaystyle\frac{{\rm d}\sigma^{2}}{{\rm
d}t}=\,\overline{\\!{\varepsilon}}\langle k\varepsilon\rangle-2\langle
k\rangle\sigma^{2}+\frac{s}{2}\langle
k\varepsilon^{2}\rangle-\frac{s}{2}\langle
k\varepsilon\rangle\langle\epsilon\rangle.$ (31)

This has the solution

$\displaystyle\sigma^{2}$ $\displaystyle=$ $\displaystyle\frac{s}{4\langle
k\rangle}\left(\langle k\varepsilon^{2}\rangle-\langle
k\varepsilon\rangle\langle\varepsilon\rangle\right)+\,\overline{\\!{\varepsilon}}(0)\frac{\langle
k\varepsilon\rangle}{\langle k\rangle}{{\rm e}}^{-\langle k\rangle t}$
$\displaystyle+$ $\displaystyle\left(\sigma^{2}(0)-\frac{s}{4\langle
k\rangle}\left(\langle k\varepsilon^{2}\rangle-\langle
k\varepsilon\rangle\langle\varepsilon\rangle\right)-\,\overline{\\!{\epsilon}}(0)\frac{\langle
k\epsilon\rangle}{\langle k\rangle}\right){{\rm e}}^{-2\langle k\rangle t}$

For the special case where all $\epsilon_{i}=\epsilon_{0}$ are identical, the
equation reduces to

$\displaystyle\sigma^{2}=\,\overline{\\!{\epsilon}}(0)\epsilon_{0}{\rm
e}^{-\langle k\rangle t}\left(1-{\rm e}^{-\langle k\rangle
t}\right)+\sigma^{2}(0){\rm e}^{-2\langle k\rangle t}.$ (33)

At long times the ratio of the standard deviation to the mean value for this
case approaches the value

$\displaystyle\frac{\sigma}{\,\overline{\\!{\epsilon}}}\approx\sqrt{\frac{\epsilon_{0}}{\,\overline{\\!{\epsilon}}(0)}}{\rm
e}^{\langle k\rangle t/2}.$ (34)

Hence the relative widths of these distributions increase as the particle
cools.

Another point worth considering is if the width can increase in absolute
terms. This will happen if the initial distribution is sufficiently narrow.
Taking the derivative at zero time shows that the condition for this is

$\displaystyle\sigma(0)\leq\sqrt{\frac{\,\overline{\\!{\epsilon}}(0)\epsilon_{0}}{2}}.$
(35)

For an ensemble with an initially canonical distribution and $s$ thermally
activated vibrational degrees of freedom, the value of $\sigma$ is
$\sqrt{s}\epsilon$. This translates into the condition

$\displaystyle\,\overline{\\!{\epsilon}}(0)\leq\frac{\epsilon_{0}}{2s}.$ (36)

This is a very small width and is rarely of interest in realistic situations.

It is of interest to compare the results for the mean and variances with those
that describe an emission with a continuum spectrum. These cases were treated
in 23 with a photon absorption cross section proportional to the photon energy
to the power $n$. Both the mean energy and the width varies asymptotically as
a power law with the power $-1/(n+3)$, i.e. in parallel, in contrast to the
exponential and non-parallel variations found here.

##  5 Discussion and Summary

The dynamics of vibrational cooling in the rapid exchange model has been
analyzed, showing that the cooling rate has a linear dependence on internal
energy. This implies that if one observes experimentally a deviation from the
linear behavior, this means that the assumption of rapid IVR would be
incorrect. 12 Furthermore, the consequences regarding the temporal evolution
of the means and the variances of energy distributions that are cooled by
vibrational transitions have been studied.

The main result is an exponential decrease of the mean energy with time. The
time dependence of the width of the energy distribution involves several terms
but will asymptotically approach an exponentially decreasing value, albeit
with a time constant of twice of the one for the mean value. Distributions
will therefore tend to broaden relative to the mean with time.

The high energy results fail at excitation energies around one quantum per
mode and less. This is not surprising and is also clear from inspection of the
literature results on C${}_{5}^{-}$ 9. The energy distributions for this
cluster decrease with a power of time, with the same power for width and mean.
Simulated data reproduce these features with measured/calculated low energy
and ground state properties, and the differences to the present results are
therefore not due to new physics, but rather a consequence of the low energy
limit. In lieu of a closed form expression for the photon emission rate
constants rate in this limit, the equations 6, and 19 can be used to calculate
the temporal development of energy distributions.

## References

  * Herbst and Dishoeck 2009 Herbst, E.; Dishoeck, E. F. V. Complex organic interstellar molecules. _Ann. Rev. Astron. Astr._ 2009, _47_ , 427–480
  * McGuire 2018 McGuire, B. A. 2018 Census of Interstellar, Circumstellar, Extragalactic, Protoplanetary Disk, and Exoplanetary Molecules. _Astrophys. J. Suppl. S._ 2018, _239_ , 17
  * Lange et al. 2012 Lange, M.; Froese, M. W.; Menk, S.; Bing, D.; Fellenberger, F.; Grieser, M.; Laux, F.; Orlov, D. A.; Repnow, R.; Sieber, T.; Toker, Y.; von Hahn, R.; Wolf, A.; Blaum, K. Radiative cooling of $Al_{4}^{-}$ and $Al_{5}^{-}$ in a cryogenic environment. _New J. Phys._ 2012, _14_ , 065007
  * Najafian et al. 2014 Najafian, K.; Pettersson, M. S.; Dynefors, B.; Shiromaru, H.; Matsumoto, J.; Tanuma, H.; Furukawa, T.; Azuma, T.; Hansen, K. Radiative cooling of $C_{7}^{-}$. _J. Chem. Phys._ 2014, _140_ , 104311
  * Breitenfeldt et al. 2018 Breitenfeldt, C.; Blaum, K.; George, S.; Göck, J.; Guzmán-Ramírez, G.; Karthein, J.; Kolling, T.; Lange, M.; Menk, S.; Meyer, C.; Mohrbach, J.; Niedner-Schatteburg, G.; Schwalm, D.; Schweikhard, L.; Wolf, A. Long-Term Monitoring of the Internal Energy Distribution of Isolated Cluster Systems. _Phys. Rev. Lett._ 2018, _120_ , 253001
  * Bull et al. 2019 Bull, J. N.; Scholz, M. S.; Carrascosa, E.; Kristiansson, M. K.; Eklund, G.; Punnakayathil, N.; de Ruette, N.; Zettergren, H.; Schmidt, H. T.; Cederquist, H.; Stockett, M. H. Ultraslow radiative cooling of Cn- (n=3-5). _J. Chem. Phys._ 2019, _151_ , 114304
  * Stockett et al. 2020 Stockett, M.; Bull, J. N.; Buntine, J. T.; Carrascosa, E.; Anderson, E. K.; Gatchell, M.; Kaminska, M.; Nascimento, R. F.; Cederquist, H.; T., S. H.; Zettergren, H. Radiative cooling of carbon cluster anions C${}_{2n+1}^{-}$ ($n=3-5$). _Eur. Phys. J. D_ 2020, _74_ , 150
  * Stockett et al. 2020 Stockett, M. H.; Bull, J. N.; Buntine, J. T.; Carrascosa, E.; Ji, M.; Kono, N.; Schmidt, H. T.; Zettergren, H. Unimolecular fragmentation and radiative cooling of isolated PAH ions: A quantitative study. _J. Chem. Phys._ 2020, _153_ , 154303
  * Goto et al. 2013 Goto, M.; Sundén, A.; Shiromaru, H.; Matsumoto, J.; Tanuma, H.; Azuma, T.; Hansen, K. Direct observation of internal energy distributions of C${}_{5}^{-}$. _J. Chem. Phys._ 2013, _139_ , 054306
  * Ferrari et al. 2019 Ferrari, P.; Janssens, E.; Lievens, P.; Hansen, K. Radiative cooling of size-selected gas phase clusters. _Int. Rev. Phys. Chem._ 2019, _38_ , 405
  * Iida et al. 2021 Iida, S.; Kuma, S.; Kuriyama, M.; Furukawa, T.; Kusunoki, M.; Tanuma, H.; Hansen, K.; Shiromaru, H.; Azuma, T. IR-photon quenching of delayed electron detachment from hot pentacene anions. _Phys. Rev. A_ 2021, _104_ , 043114
  * 12 D. Muell, Y. T. e. a., H. Kreckel In preparation: Radiative cooling of Al${}_{4}^{-}$ and Al${}_{5}^{-}$ in a cryogenic environment.
  * Bixon and Jortner 1968 Bixon, M.; Jortner, J. Intramolecular Radiationless Transitions. _J. Chem. Phys._ 1968, _48_ , 715
  * Gruebele and Wolynes 2004 Gruebele, M.; Wolynes, P. G. Vibrational Energy Flow and Chemical Reactions. _Acc. Chem. Res._ 2004, _37_ , 261–267
  * Makarov et al. 2012 Makarov, A. A.; Malinovsky, A. L.; Ryabov, E. A. Intramolecular vibrational redistribution: from high-resolution spectra to real-time dynamics. _Uspekhi Fizicheskih Nauk_ 2012, _182_ , 977
  * Nesbitt and Field 1996 Nesbitt, D. J.; Field, R. W. Vibrational Energy Flow in Highly Excited Molecules: Role of Intramolecular Vibrational Redistribution. _J. Phys. Chem._ 1996, _100_ , 12735–12756
  * Hansen and Ferrari 2021 Hansen, K.; Ferrari, P. Vibrational angular momentum level densities of linear molecules. _Chem. Phys. Lett._ 2021, _768_ , 138385
  * Brink and Stringari 1990 Brink, D. M.; Stringari, S. Density of states and evaporation rate of helium clusters. _Z Phys D_ 1990, _15_ , 257–263
  * Lehmann 2004 Lehmann, K. K. Bias in the temperature of helium nanodroplets measured by an embedded rotor. _J. Chem. Phys._ 2004, _120_ , 513
  * Hansen et al. 2007 Hansen, K.; Johnson, M. D.; Kresin, V. V. Density of states of helium droplets. _Phys. Rev. B_ 2007, _76_ , 235424
  * Chandrasekaran et al. 2014 Chandrasekaran, V.; Kafle, B.; Prabhakaran, A.; Heber, O.; Rappaport, M.; Rubinstein, H.; Schwalm, D.; Toker, Y.; Zajfman, D. Determination of Absolute Recurrent Fluorescence Rate Coefficients for C6–. _J. Phys. Chem. Lett._ 2014, _5_ , 4078–4082
  * Terzieva and Herbst 2000 Terzieva, R.; Herbst, E. Radiative electron attachment to small linear carbon clusters and its significance for the chemistry of diffuse interstellar clouds. _Int. J. Mass Spetrom._ 2000, _201_ , 135–142
  * Stenfalk and Hansen 2007 Stenfalk, J.; Hansen, K. Energy distributions of clusters cooled by thermal radiation. _Eur. Phys. J. D_ 2007, _43_ , 101