text:
# submitted to The Astrophysical Journal Infrared Emission from Interstellar Dust. I. Stochastic Heating of Small Grains

## 1 Introduction

Very small grains – consisting of tens to hundreds of atoms – play a major role in the infrared emission from the interstellar medium. These grains are small enough that the time-averaged vibrational energy $`E`$ is smaller than or comparable to the energy of the starlight photons which heat the grains. Stochastic heating by absorption of starlight therefore results in transient “temperature spikes”, during which much of the energy deposited by the starlight photon is reradiated in the infrared.

The idea of transient heating of very small grains was first introduced by Greenberg (1968). Its importance was not apparent until the detection of the near infrared emission of reflection nebulae (Sellgren et al. 1983) and detection by the Infrared Astronomical Satellite (IRAS) of 12 and 25$`\mu `$m Galactic emission which was far in excess of the emission expected for interstellar dust at $`T20\mathrm{K}`$ (Boulanger & Pérault 1988). Subsequent measurements by the Diffuse Infrared Background Experiment (DIRBE) instrument on the Cosmic Background Explorer (COBE) satellite confirmed this and detected additional broadband emission at 3.5 and 4.9$`\mathrm{\mu m}`$ (Arendt et al. 1998). More recently, spectrometers aboard the Infrared Telescope in Space (IRTS) (Onaka et al. 1996; Tanaka et al. 1996) and the Infrared Space Observatory (ISO) (Mattila et al. 1996) have shown that the diffuse interstellar medium radiates strongly in emission features at 3.3, 6.2, 7.7, 8.6, and 11.3$`\mathrm{\mu m}`$. These emission features, previously observed from a small number of bright reflection nebulae, planetary nebulae, and HII regions, have been tentatively identified as emission from polycyclic aromatic hydrocarbons, or PAHs (Léger & Puget 1984; Allamandola, Tielens, & Barker 1985).

To understand these observations we need to be able to calculate the infrared emission expected for a given model of interstellar dust. The question we address here is: what approximations can be used to carry out accurate calculations of thermal emission from a population of very small grains exposed to stochastic heating by starlight?

There have been a number of previous studies of the discrete heating of very small grains (Duley 1973; Greenberg & Hong 1974; Gail & Seldmayr 1975; Greenberg 1976; Purcell 1976; Drapatz & Michel 1977; Andriesse 1978; Aannestad & Kenyon 1979; Léger & Puget 1984; Draine & Anderson 1985; Puget et al. 1985; Désert et al. 1986; Dwek 1986; Guhathakurta & Draine 1989; Aannestad 1989; Siebenmorgen et al. 1992; Manske & Henning 1998; Duley & Poole 1998). Most studies have made three approximations: (1) the “thermal” approximation for excitation; (2) use of “bulk” specific heats; and (3) the “continuous cooling” approximation. While accurate for large grains, these approximations could potentially introduce significant inaccuracies for very small grains.

Barker & Cherchneff (1989) and d’Hendecourt et al. (1989) found that the thermal approximation was valid for computing thermal emission spectra. Allamandola, Tielens, & Barker (1989), Schutte, Tielens, & Allamandola (1993), and Cook & Saykally (1998) compared calculations of infrared emission from PAHs using the thermal approximation with the results of a more detailed statistical treatment; they found that the thermal approximation provided a good approximation for the emission from low-frequency modes, although it overestimated the emission from high frequency modes such as the $`3.3\mathrm{\mu m}`$ C-H stretching mode. All of these studies used the “continuous cooling” approximation to discuss the overall emission from a PAH following single-photon heating.

The validity of the “continuous cooling” approximation has been questioned by Siebenmorgen et al. (1992), Manske & Henning (1998), and Duley & Poole (1998). Siebenmorgen et al. (1992) and Manske & Henning (1998) attempted to treat the cooling process as discrete events, but did not correctly evaluate the transition probabilities, as we show below in §6.1.

The present paper addresses the vibrational excitation of interstellar grains, both small and large, and the resulting infrared emission. If energy levels and transition probabilities were known, we could (at least in principle) solve for the statistical steady-state populations of the different energy levels of grains illuminated by a known radiation field. However, this level of detailed information is generally unavailable, for even the smallest and simplest PAH molecules. We therefore construct a “model PAH” grain, CnHm, with a realistic spectrum of energy levels. Following previous workers (e.g., Schutte et al. 1993, Boulanger et al. 1998) we adopt realistic absorption cross section $`C_{\mathrm{abs}}(\lambda )`$, depending on the number of C atoms, the H/C ratio, and whether the PAH is neutral or ionized. We construct a similar model for silicate grains, with appropriate mode spectrum and absorption cross section. We divide the energy level spectrum into a manageable number of “bins”, estimate the “bin-to-bin” transition probabilities, and can then solve for the steady-state probability $`P_j`$ of occupying energy bin $`j`$.

Our analysis allows for discrete quantum states at low vibrational energies, and what is effectively a continuum of vibrational levels at large energies. As we show, the only simplifying assumptions that need to be made are: (1) the absorption cross section $`C_{\mathrm{abs}}(\lambda )`$ is independent of the degree of vibrational excitation, and (2) internal vibrational relaxation rapidly redistributes internal energy among the vibrational degrees of freedom of the grain. While not exact, these are both expected to be good approximations. Once $`C_{\mathrm{abs}}(\lambda )`$ is adopted and the spectrum of vibrational “normal modes” (and from this the vibrational density of states) is specified, the physics of absorption and emission of radiation by the grain is essentially fully determined. With these assumptions, we can solve for the energy distribution function $`P(E)`$ for a given grain in an arbitrary radiation field, with full inclusion of “multiphoton heating” effects which become important for large grains or intense radiation fields. Because we solve for $`P(E)`$ including all radiative transitions, we obtain the “exact” solution against which we can then test solutions obtained using the thermal approximation and the continous cooling approximation.

In §2 we estimate the spectrum of vibrational modes for small PAH and silicate grains, and obtain the vibrational density of states for such grains. The transition probabilities depend on the photon absorption cross sections, and in §3 we describe the cross sections adopted for PAH grains.

We formulate the “exact” statistical treatment in §4. In §5 we discuss how the “temperature” of a grain is related to its vibrational energy content, and in §6 we show how the “thermal” approximation is related to the exact treatment. The radiative cooling time for vibrationally-excited grains is estimated as a function of $`E`$ and grain size in §7.

In §8 we discuss the method used to solve for the steady-state energy distribution. In §9 we show how the infrared emission spectrum can be calculated from the energy distribution, both for the exact statistical treatment and using the thermal approximation. We show some sample solutions in §10, where we assess the accuracy of the “thermal” approximation. We find the thermal approximation to be quite accurate for calculation of the overall emission spectrum from a population of very small grains. Furthermore, we find that the “continuous cooling” approximation is itself quite accurate, and can therefore be used to greatly reduce the amount of computation required to find the temperature distribution functions.

The above studies therefore provide a methodology for computing the infrared emission from dust grains irradiated by starlight. The method is accurate, and can be used for any specified mixture of grain sizes and compositions; we present results for both carbonaceous and silicate materials. We discuss our results in §11 and summarize the main points in §12.

The paper is largely concerned with the physics of infrared emission from stochastically-heated grains, and practical methods for calculating the emission spectrum. Readers who are interested only in the resulting emission spectrum may wish to proceed directly to §10. In a following paper (Li & Draine 2001a) we apply this methodology to construct models to account for the observed infrared emission from the diffuse interstellar medium.

## 2 Vibrational Mode Spectrum and Level Degeneracies

A grain containing $`N_a`$ atoms will have 3 translational degrees of freedom, 3 rotational degrees of freedom, and $`N_m=3N_a6`$ vibrational modes. Each vibrational mode $`j=1,\mathrm{},N_m`$ has a fundamental frequency $`\omega _j`$. If we approximate the vibrational modes as harmonic oscillators, then the vibrational energy of the grain (not including zero-point energy) is

$$E=\underset{j=1}{\overset{N_m}{}}v_j\mathrm{}\omega _j,$$
(1)
where the integer $`v_j0`$ is the vibrational quantum number for mode $`j`$. If the vibrational system is in thermal equilibrium at temperature $`T`$, the expectation value of the energy $`E`$ is

$$\overline{E}(T)=\underset{j=1}{\overset{N_m}{}}\frac{\mathrm{}\omega _j}{\mathrm{exp}(\mathrm{}\omega _j/kT)1}.$$
(2)
where $`k`$ is Boltzmann’s constant. The heat capacity

$$C(T)\frac{d\overline{E}}{dT}=k\underset{j=1}{\overset{N_m}{}}e^{\mathrm{}\omega _j/kT}\left[\frac{\mathrm{}\omega _j/kT}{1\mathrm{exp}(\mathrm{}\omega _j/kT)}\right]^2.$$
(3)
Ideally, we would have accurate knowledge of the $`N_m`$ frequencies $`\omega _j`$. Since the complete mode spectrum is known only for a small number of grain candidates such as C24H12 coronene (see below), we need to estimate the normal mode spectrum for candidate grains of different sizes and materials.

### 2.1 Normal Modes: Polycyclic Aromatic Hydrocarbons

A polycyclic aromatic hydrocarbon molecule containing $`N_\mathrm{C}`$ C atoms and $`N_\mathrm{H}`$ H atoms has $`N_m=3(N_\mathrm{H}+N_\mathrm{C}2)`$ distinct vibrational modes. These can be separated into $`3(N_\mathrm{C}2)`$ modes of the C-C skeleton, and $`3N_\mathrm{H}`$ modes associated with the C-H bonds at the periphery. The frequencies of these normal modes have been computed for a small number of polycyclic aromatic hydrocarbons, with some frequencies determined experimentally, but mode spectra are not yet available for most PAHs of interest. Here we estimate the mode spectrum for “generic” PAHs. We will treat 5 different types of vibration separately: out-of-plane C-C modes, in-plane C-C modes, out-of-plane C-H bending, in-plane C-H bending, and C-H stretching.

Graphite is a useful guide to the C-C modes. Krumhansl & Brooks (1953) found that the lattice vibration specific heat of graphite could be approximated by a model where the out-of-plane and in-plane C-C vibrations each had a mode spectrum given by a two-dimensional Debye model, with Debye temperatures $`\mathrm{\Theta }_{op}950\mathrm{K}`$ for the out-of-plane modes, and $`\mathrm{\Theta }_{ip}2500\mathrm{K}`$ for the in-plane modes. This suggests that the C-C vibrational modes for a PAH might be similarly approximated.

A $`n`$-dimensional Debye spectrum has modes uniformly distributed in $`E^n`$ up to a maximum energy $`k\mathrm{\Theta }`$. Consider a mode spectrum of the form

$$\mathrm{}\omega _j=k\mathrm{\Theta }_D\left[\frac{(1\beta _n)}{N_m}(j\delta _j)+\beta _n\right]^{1/n}j=1,\mathrm{},N_m.$$
(4)
For the PAH C-C modes we set the dimensionality $`n=2`$, with $`N_m=N_\mathrm{C}2`$ for the out-of-plane modes, and $`N_m=2(N_\mathrm{C}2)`$ for the in-plane modes.

If $`\delta _j`$ is independent of $`j`$, the modes are uniformly distributed in $`E^2`$. We will in fact take

$`\delta _j`$ $`=`$ $`1/2\mathrm{for}j2\mathrm{or}3`$ (5)
$`=`$ $`1\mathrm{for}j=2\mathrm{or}3`$ (6)
because this shift of modes 2 and 3 to lower frequency improves the match to the actual mode spectrum of coronene.

As discussed in Appendix A we assume PAHs to be planar for $`N_\mathrm{C}54`$, approximately spherical for $`N_\mathrm{C}>102`$, with intermediate shape for $`54<N_\mathrm{C}<102`$. Thus for $`N_\mathrm{C}54`$ we expect the lowest frequency mode $`\omega _1N_\mathrm{C}^{1/2}`$, approximately constant $`\omega _1`$ for $`54<N_\mathrm{C}<102`$, and $`\omega _1N_\mathrm{C}^{1/3}`$ for $`N_\mathrm{C}>102`$. We achieve this by taking

$`\beta _2`$ $`=`$ $`0\mathrm{for}N_\mathrm{C}54`$ (7)
$`=`$ $`{\displaystyle \frac{1}{(2N_m1)}}\left({\displaystyle \frac{N_\mathrm{C}54}{52}}\right)\mathrm{for}54<N_\mathrm{C}102`$
$`=`$ $`{\displaystyle \frac{1}{(2N_m1)}}\left[{\displaystyle \frac{(N_\mathrm{C}2)}{52}}\left({\displaystyle \frac{102}{N_\mathrm{C}}}\right)^{2/3}1\right]\mathrm{for}N_\mathrm{C}>102`$
We use equations (4-7) with $`N_m=N_\mathrm{C}2`$ and $`k\mathrm{\Theta }_{op}/hc=600\mathrm{cm}^1`$ ($`\mathrm{\Theta }_{op}=863\mathrm{K}`$) for the out-of-plane C-C modes, and $`N_m=2(N_\mathrm{C}2)`$ and $`k\mathrm{\Theta }_{ip}/hc=1740\mathrm{cm}^1`$ ($`\mathrm{\Theta }_{ip}=2504\mathrm{K}`$) for the in-plane C-C modes.

We can adequately approximate the C-H stretching and bending modes by $`N_\mathrm{H}`$ out-of-plane bending modes at $`\lambda _{\mathrm{CH},op}^1=(11.3\mathrm{\mu m})^1=886\mathrm{cm}^1`$, $`N_\mathrm{H}`$ in-plane bending modes at $`\lambda _{\mathrm{CH},ip}^1=(8.6\mathrm{\mu m})^1=1161\mathrm{cm}^1`$, and $`N_\mathrm{H}`$ C-H stretching modes at $`\lambda _{\mathrm{CH},str}^1=(3.3\mathrm{\mu m})^1=3030\mathrm{cm}^1`$.

For our model “astronomical PAH” sequence, we will assume

$`N_\mathrm{H}`$ $`=`$ $`\mathrm{int}\left(0.5N_\mathrm{C}+0.5\right)N_\mathrm{C}25`$ (8)
$`=`$ $`\mathrm{int}\left(2.5\sqrt{N_\mathrm{C}}+0.5\right)25N_\mathrm{C}100`$
$`=`$ $`\mathrm{int}\left(0.25N_\mathrm{C}+0.5\right)N_\mathrm{C}100`$
where $`\mathrm{int}(x)`$ is the integer part of $`x`$. The pericondensed species coronene C24H12, circumcoronene C54H18, and dicircumcoronene C96H24 are members of the sequence prescribed by eq. (8). Beyond $`N_\mathrm{C}=100`$ we assume H/C $`0.25`$. For a carbonaceous grain containing $`N_\mathrm{C}`$ carbon atoms, we associate a representative “radius” $`a=10\mathrm{\AA }(N_\mathrm{C}/468)^{1/3}`$, corresponding to the carbon density $`\rho =2.24\mathrm{g}\mathrm{cm}^3`$ of graphite.

We use C24H12 coronene, for which the mode spectrum is known (Cyvin 1982; Cyvin et al. 1984), to test the above model. As shown in Figure 1, our “synthetic” mode spectrum is in excellent agreement with the actual mode spectrum of coronene.

Applying our synthetic mode spectrum for PAHs in the limit $`N_\mathrm{C}\mathrm{}`$ and $`N_\mathrm{H}=0`$, the mode spectrum of the C-C modes reduces to the 2-dimensional Debye model. The heat capacity becomes

$$C_{\mathrm{graph}}=(N_\mathrm{C}2)k\left[f_2^{}(T/863\mathrm{K})+2f_2^{}(T/2504\mathrm{K})\right]$$
(9)
$$f_n(x)\frac{1}{n}_0^1\frac{y^ndy}{\mathrm{exp}(y/x)1},f_n^{}(x)\frac{d}{dx}f_n(x).$$
(10)
In Figure 2 we compare eq. (9) with experimental data for (bulk) graphite. The excellent agreement confirms that our model spectrum for the C-C vibrational modes applies from small PAH molecules up through macroscopic samples of graphite.

### 2.2 Normal Modes: Silicate Grains

In the case of silicate grains, we use experimental specific heats as a guide to the vibrational mode spectrum. The specifics heat of SiO2 glass, obsidian glass (75% SiO2 and 25% metal oxides by mass), and basalt glass (50% SiO2, 50% metal oxides by mass) reported by Léger, Jura, & Omont (1985) are plotted Figure 2. As seen in Figure 2, the experimental results for bulk silicates (basalt and obsidian) can be satisfactorily reproduced if 2/3 of the vibrational modes are distributed according to a Debye model with $`n=2`$ and Debye temperature $`\mathrm{\Theta }=500\mathrm{K}`$, and 1/3 of the modes are described by a Debye model with $`n=3`$ and $`\mathrm{\Theta }=1500\mathrm{K}`$:

$$C_{\mathrm{sil}}=(N_a2)k\left[2f_2^{}(T/500\mathrm{K})+f_3^{}(T/1500\mathrm{K})\right]$$
(11)
where $`N_a`$ is the number of atoms in the cluster, and $`f_n`$ is given by eq. (10). Eq. (11) agrees quite well with the measured specific heats for obsidian and basalt, so we use this model to estimate the mode spectrum for silicate grains.

We expect the lowest frequency mode $`\omega _1N_a^{1/3}`$. To achieve this scaling for modes distributed according to an $`n`$dimensional Debye model \[eq. (4)\], we set

$$\beta _n=\frac{N_m^{1n/3}1}{2N_m1};$$
(12)
note that for $`n=3`$ this gives $`\beta _3=0`$.

For small silicate grains, two-thirds of the vibrational modes are then approximated by eqs. (4-6) with $`n=2`$, $`\mathrm{\Theta }_D=500\mathrm{K}`$, and $`N_m=2(N_a2)`$, and the remaining third by eqs. (4-6) with $`n=3`$, $`\mathrm{\Theta }_D=1500\mathrm{K}`$, and $`N_m=N_a2`$. The resulting normal mode spectrum is shown in Figure 3.

### 2.3 Density of States

Let there be a total of $`N_m=3(N_a2)`$ distinct normal modes. The general vibrational state can then be specified by the $`N_m`$-tuple $`(v_1,v_2,\mathrm{},v_{N_m})`$ giving the vibrational quantum numbers for each of these modes. If each mode is approximated by a harmonic oscillator, then $`N(E)`$, the number of distinct vibrational states with total vibrational energy less than or equal to $`E`$, can be calculated using the Beyer-Swinehart algorithm (Beyer & Swinehart 1973; Stein & Rabinovitch 1973). In Figure 4 we show $`N(E)`$ computed for C24H12, using both the actual normal mode spectrum for coronene and our model normal mode spectrum for C24H12 (see Figure 1). The resulting $`N(E)`$ are essentially identical for $`E/hc300\mathrm{cm}^1`$.

The Beyer-Swinehart algorithm is remarkably efficient, and with IEEE standard 64-bit arithmetic can be applied to calculate $`N(E)`$ up to $`E/hc=10^5\mathrm{cm}^1`$ for clusters with up to $``$600 atoms. For a given cluster, we employ the Beyer-Swinehart algorithm for energies up to an energy $`E_t`$ where $`N(E_t)10^{300}`$, beyond which we are limited by inability to calculate floating point numbers exceeding $`2^{1024}10^{308}`$. For $`E>E_t`$ we calculate the density of states by integrating $`dS=dQ/T`$:

$$\mathrm{ln}N(E)=\mathrm{ln}N(E_t)+_{E_t}^E\frac{dE^{}}{kT(E^{})},$$
(13)
where $`T(E)`$ is the temperature at which the grain has energy $`E`$ (see §5). Figures 4 and 5 show $`N(E)`$ evaluated for C24H12, C400H100, C800H200, and C4000H1000 using our model density of states.

## 3 Absorption Cross Sections

Following previous workers (e.g., Schutte et al. 1993), we estimate the absorption cross sections for both neutral and ionized PAHs from far ultraviolet to far infrared (see Li & Draine 2001a for details) based on available experimental data (Allamandola et al. 1999; Hudgins & Allamandola 1999; and references therein) and guided by astronomical observations (e.g., Boulanger et al. 1998). In the ultraviolet and infrared the resulting cross sections are mainly a collection of Drude profiles: the $`\sigma \sigma ^{}`$ transition mode ($`\lambda ^114\mu `$m-1); the $`\pi \pi ^{}`$ transition ($`\lambda ^14.6\mu `$m-1); the C-H stretching mode ($`\lambda =3.3\mu `$m); the C-C stretching modes ($`\lambda =6.2,7.7\mu `$m); the C-H in-plane bending mode ($`\lambda =8.6\mu `$m); the C-H out-of-plane bending modes ($`\lambda =11.3,11.9,12.7\mu `$m); and a few weak features attributed to the C-C bending modes ($`\lambda =16.4,18.3,21.2,23.1\mu `$m). In addition to these features, we have included continuous absorption in the far-UV, near-UV/visible, as well as a small amount of continuous absorption in the infrared. \[For further information see Li & Draine (2001)\]. In Figure 6 we display the adopted absorption cross sections for neutral and ionized PAHs.

Following Li & Draine (2001), we make a smooth transition from PAH optical properties for $`N_\mathrm{C}6\times 10^4`$ ($`a<50\mathrm{\AA }`$) to graphite properties for $`N_\mathrm{C}5\times 10^5`$ ($`a>100\mathrm{\AA }`$), calculated using Mie theory with the 1/3-2/3 approximation (Draine & Malhotra 1993) and dielectric functions from Draine & Lee (1984).

For silicates, we use Mie theory with dielectric constants for astronomical silicate (Draine & Lee 1984).111 We use the “smoothed UV” silicate dielectric function discussed by Weingartner & Draine 2001a, with modifications to the far-infrared emissivity proposed by Li & Draine 2001a, available at http://www.astro.princeton.edu/$``$draine .

## 4 Excitation of Very Small Grains: Exact Statistical Treatment

### 4.1 Assumptions

We characterize the state of the grain by its vibrational energy $`E`$. There are too many energy levels to consider individually, so we group them into $`M+1`$ “bins” $`j=0,\mathrm{},M`$, where the $`j`$-th bin is $`[E_{j,\mathrm{min}},E_{j,\mathrm{max}})`$, with representative energy $`E_j(E_{j,\mathrm{min}}+E_{j,\mathrm{max}})/2`$, and width $`\mathrm{\Delta }E_jE_{j,\mathrm{max}}E_{j,\mathrm{min}}`$. The “degeneracy” $`g_j`$ is the number of distinct quantum states included in bin $`j`$. The procedure used for specifying the bins is described in Appendix B.

Let $`P_j`$ be the probability of finding the grain in bin $`j`$. The probability vector $`P_j`$ evolves according to

$$\frac{d}{dt}P_i=\underset{ji}{}𝐓_{ij}P_j\underset{ji}{}𝐓_{ji}P_ii=0,1,\mathrm{},M,$$
(14)
where the transition matrix element $`𝐓_{fi}`$ is the probability per unit time for a grain in bin $`i`$ to make a transition to one of the levels in bin $`f`$. All of the physics is embodied in the transition matrix $`𝐓`$.

A truly exact treatment of the grain excitation would require knowledge of all the state-to-state transition probabilities between energy levels. Given the large number of energy levels to be considered, this is utterly infeasible at this epoch. In order to estimate the required state-to-state rates, we make only the following three assumptions:

* The grain absorption cross section $`C_{\mathrm{abs}}(h\nu )`$ depends on the photon energy $`h\nu `$, but does not depend on the vibrational energy of the grain prior to the absorption.
* The vibrational modes of the grain can be approximated by harmonic oscillators.
* The energy of absorbed photons is distributed ergodically among the $`3N_a6`$ vibrational degrees of freedom before any infrared emission takes place.222 Since the grain will in general have angular momentum, and the principal values of the moment of inertia tensor will generally be nondegenerate, in principle some of the heat can appear in rotational excitation (Lazarian & Draine 1999a,b). However, since there are only three rotational degrees of freedom, this is at most a small effect, and will be neglected.

We now show that with only these three assumptions, together with knowledge of

* the spectrum of fundamental vibrational modes,
* the grain absorption section $`C_{\mathrm{abs}}(h\nu )`$, and
* the starlight energy density $`u_E`$,

we can fully determine the emission from interstellar dust grains. Note that the notion of “grain temperature” does not appear in what we will refer to as the “exact-statistical treatment”.

### 4.2 Upward Transitions

In Appendix C we show that the rate for upward transitions $`lu`$ is

$`𝐓_{ul}`$ $`=`$ $`{\displaystyle \frac{c\mathrm{\Delta }E_u}{E_uE_l}}{\displaystyle _{W_1}^{W_4}}G_{ul}(E)C_{\mathrm{abs}}(E)u_E𝑑E\mathrm{for}u<M`$ (15)
$`𝐓_{Ml}`$ $`=`$ $`{\displaystyle \frac{c}{E_ME_l}}\left[{\displaystyle _{W_1}^{W_c}}\left({\displaystyle \frac{EW_1}{W_cW_1}}\right)C_{\mathrm{abs}}(E)u_E𝑑E+{\displaystyle _{W_c}^{\mathrm{}}}C_{\mathrm{abs}}(E)u_E𝑑E\right]`$ (16)
$`G_{ul}(E)`$ $``$ $`{\displaystyle \frac{(EW_1)}{\mathrm{\Delta }E_u\mathrm{\Delta }E_l}}\mathrm{for}W_1<E<W_2`$ (17)
$``$ $`{\displaystyle \frac{\mathrm{min}[\mathrm{\Delta }E_l,\mathrm{\Delta }E_u]}{\mathrm{\Delta }E_u\mathrm{\Delta }E_l}}\mathrm{for}W_2<E<W_3`$ (18)
$``$ $`{\displaystyle \frac{(W_4E)}{\mathrm{\Delta }E_u\mathrm{\Delta }E_l}}\mathrm{for}W_3<E<W_4`$ (19)
$``$ $`0\mathrm{for}E<W_1\mathrm{or}E>W_4`$ (20)
$`W_1`$ $``$ $`E_{u,\mathrm{min}}E_{l,\mathrm{max}}`$ (21)
$`W_2`$ $``$ $`\mathrm{min}[(E_{u,\mathrm{min}}E_{l,\mathrm{min}}),(E_{u,\mathrm{max}}E_{l,\mathrm{max}})]`$ (22)
$`W_3`$ $``$ $`\mathrm{max}[(E_{u,\mathrm{min}}E_{l,\mathrm{min}}),(E_{u,\mathrm{max}}E_{l,\mathrm{max}})]`$ (23)
$`W_4`$ $``$ $`E_{u,\mathrm{max}}E_{l,\mathrm{min}}`$ (24)
$`W_c`$ $``$ $`E_{M,\mathrm{min}}E_{l,\mathrm{min}}`$ (25)
and $`u_EdE`$ is the energy density due to photons with energies in $`[E,E+dE]`$. Note that eq. (16) takes energy absorbed in transitions to levels beyond the highest bin and allocates it to the highest bin ($`M`$). In practice, we ensure that the highest bin is set high enough that the rate of such transitions is negligible, and the population $`P_M`$ is negligibly small.

The factor $`G_{ul}(E)`$ represents the correction for finite bin width; if $`\mathrm{max}(\mathrm{\Delta }E_l,\mathrm{\Delta }E_u)(E_uE_l)`$, then eq. (15) becomes

$`𝐓_{ul}`$ $``$ $`cC_{\mathrm{abs}}(E){\displaystyle \frac{u_E\mathrm{\Delta }E_u}{E_uE_l}}\mathrm{for}u<M`$ (26)
$`𝐓_{Ml}`$ $``$ $`cC_{\mathrm{abs}}(E){\displaystyle \frac{u_E\mathrm{\Delta }E_M}{E_ME_l}}+{\displaystyle \frac{c}{E_ME_l}}{\displaystyle _{E_ME_l}^{\mathrm{}}}C_{\mathrm{abs}}(E)u_E𝑑E.`$ (27)
Most previous studies (e.g., Guhathakurta & Draine 1989) have used this approximate form, but eq. (15-24) make proper allowance for finite bin width and are to be preferred. Correction for finite bin width is important when the treatment is applied to grains with radii $`a50\mathrm{\AA }`$, since for tractable numbers of bins $`M10^3`$, many bins will necessarily be broad, with $`\mathrm{\Delta }E_ukT_u`$.

For the special case of transitions $`u1u`$ we include “intrabin” absorptions:

$$𝐓_{u,u1}=\frac{c}{E_uE_{u1}}\left[_0^{W_4}\mathrm{\Delta }E_uG_{u,u1}(E)C_{\mathrm{abs}}(E)u_E𝑑E+_0^{\mathrm{\Delta }E_{u1}}\left(1\frac{E}{\mathrm{\Delta }E_{u1}}\right)C_{\mathrm{abs}}(E)u_E𝑑E\right]$$
(28)
but the “intrabin” contribution (second integral in eq. 28) is negligible in all cases of interest.

### 4.3 Downward Transitions

The radiative coupling between levels is determined by the absorption cross section $`C_{\mathrm{abs}}(h\nu )`$. The Einstein $`A`$ coefficient for radiative decay $`ul`$, averaged over the sublevels in bin $`u`$, is (see Appendix C)

$$𝐓_{lu}=\frac{8\pi }{h^3c^2}\frac{g_l}{g_u}\frac{\mathrm{\Delta }E_u}{E_uE_l}_{W_1}^{W_4}G_{ul}(E)E^3C_{\mathrm{abs}}(E)\left[1+\frac{h^3c^3}{8\pi E^3}u_E\right]𝑑E\mathrm{for}l<u1.$$
(29)
The degeneracies $`g_u`$ and $`g_l`$ are the numbers of energy states in bins $`u`$ and $`l`$, respectively (see §2.3):

$$g_jN(E_{j,\mathrm{max}})N(E_{j,\mathrm{min}})\left(\frac{dN}{dE}\right)_{E_j}\mathrm{\Delta }E_j.$$
(30)
The term proportional to $`u_E`$ in eq. (29) is the contribution of stimulated emission. For transitions $`uu1`$ to the adjacent bin we have $`W_1=0`$: even very long wavelength emission can contribute if $`C_{\mathrm{abs}}>0`$. For this case we must also include the power radiated in “intrabin” transitions within bin $`u`$:

$`𝐓_{u1,u}`$ $`=`$ $`{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{g_{u1}}{g_u}}{\displaystyle \frac{\mathrm{\Delta }E_u}{E_uE_{u1}}}{\displaystyle _0^{W_4}}G_{u,u1}(E)E^3C_{\mathrm{abs}}(E)\left[1+{\displaystyle \frac{h^3c^3}{8\pi E^3}}u_E\right]𝑑E+`$ (31)
$`{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{1}{E_uE_{u1}}}{\displaystyle _0^{\mathrm{\Delta }E_u}}\left(1{\displaystyle \frac{E}{\mathrm{\Delta }E_u}}\right)E^3C_{\mathrm{abs}}(E)\left[1+{\displaystyle \frac{h^3c^3}{8\pi E^3}}u_E\right]𝑑E`$
Thus we see that we require only $`C_{\mathrm{abs}}(E)`$, the degeneracies $`g_j`$, and the starlight spectrum $`u_E`$ to completely determine the transition matrix $`𝐓_{fi}`$. We stress that “grain temperature” does not enter the exact-statistical treatment (e.g., eq. 15,29). In §6 below we will discuss “thermal” approximations to estimate the downward transition rates $`𝐓_{lu}`$. We first discuss an estimate $`T(E)`$ for the “temperature” of a state of specified vibrational energy $`E`$.

## 5 Energy vs. Temperature

### 5.1 “Exact” Mode Spectrum

Consider a grain with vibrational energy $`E_u`$. If the grain is large, we know that the emission can be estimated from the grain “temperature”. We can use eq. (2) to define a “temperature” $`\stackrel{~}{T}(E_u)`$ such that $`\overline{E}(\stackrel{~}{T})=E_u`$, where $`\overline{E}(T)`$ is the expectation value for the vibrational energy $`E`$ when in contact with a heat bath at temperature $`T`$.

When there are only a few vibrational quanta in the grain, however, “temperature” is not a well-defined concept. As the most extreme example, consider a grain where only the first vibrational mode is excited, with exactly one quantum: $`E_u=\mathrm{}\omega _1`$. If the grain has many vibrational degrees of freedom, then the temperature $`\stackrel{~}{T}`$ may be quite low (since it is defined in terms of the expectation value of the vibrational energy summed over all the modes). However, our grain is vibrating in the fundamental vibrational mode and we want to estimate the emission from this mode. One estimate for the temperature characterizing the excitation of this mode would be the temperature for which the occupation number of this mode is equal to 1. We will see in §6 that when this “temperature” is used in the “thermal approximation”, one obtains the correct rate of emission in the $`10`$ transition. In fact, we will use this value of the temperature for all energies up to the energy of the 20th vibrational mode. For large grains and large energies, however, $`\stackrel{~}{T}`$ is appropriate. To encompass these two limits, we take the temperature $`T_u`$ of vibrational level $`u`$ to be333 For $`T=\mathrm{}\omega /k\mathrm{ln}2`$ a quantized harmonic oscillator has occupation number equal to 1.

$`T_u`$ $`=`$ $`{\displaystyle \frac{\mathrm{}\omega _1}{k\mathrm{ln}2}}\mathrm{for}E_u\mathrm{}\omega _{20}`$ (32)
$`=`$ $`\stackrel{~}{T}(E_u)\mathrm{for}E_u>\mathrm{}\omega _{20}`$
While the choice of $`\mathrm{}\omega _{20}`$ as the energy separating these two approximations is arbitrary, we remind the reader that the very notion of “temperature” is questionable for such low degrees of excitation. Our use of eq. (32) is in fact guided by the fact that the “thermal approximation” of §6, when used with this prescription for $`T_u`$, gives emission rates in fair agreement with the “exact-statistical” treatment (29). Figures 7 and 8 show our adopted $`T(E)`$ for selected PAHs and silicates.444 The discontinuity in $`T(E)`$ at $`E=\mathrm{}\omega _{20}`$ is obviously unphysical, but we note that the discontinuity is only appreciable for large grains at very low energies, and these make a negligible contribution to the overall emission spectrum.

### 5.2 Debye Mode Spectrum

When the number of modes is large, the summation in equation (2) contains many terms. In this case, for PAHs we replace the sums over the C-C modes by the (continuum) Debye model discussed above:

$$\overline{E}_{\mathrm{PAH}}(T)N_\mathrm{H}\underset{j=1}{\overset{3}{}}\frac{\mathrm{}\omega _j}{\mathrm{exp}(\mathrm{}\omega _j/kT)1}+(N_\mathrm{C}2)\left[k\mathrm{\Theta }_{op}f_2(T/\mathrm{\Theta }_{op})+2k\mathrm{\Theta }_{ip}f_2(T/\mathrm{\Theta }_{ip})\right],$$
(33)
where $`f_2`$ is defined by eq. (10), and the summation index $`j=13`$ runs over the C-H out-of-plane bending, in-plane bending, and stretching modes (see §2.1). Similarly, for silicate grains we take

$$\overline{E}_{\mathrm{sil}}(T)(N_a2)\left[2k\mathrm{\Theta }_2f_2(T/\mathrm{\Theta }_2)+k\mathrm{\Theta }_3f_3(T/\mathrm{\Theta }_3)\right]$$
(34)
where $`\mathrm{\Theta }_2=500\mathrm{K}`$ and $`\mathrm{\Theta }_3=1500\mathrm{K}`$. We use eq. (33) and (34) for $`E>N_aE_1`$.

## 6 Thermal Approximations for Spontaneous Emission

Given the above assumptions about the vibrational mode spectrum (§2), the assumption that the absorption cross section depends only on the photon energy (§3), and the assumption of rapid internal vibrational redistribution, the “statistical” approach described above (§4) is “exact”. It is often desirable, however, to use an alternative “thermal” approximation which is less computationally demanding.

The thermal approach has been frequently used (Greenberg 1976; Purcell 1976; Aannestad & Kenyon 1979; Léger & Puget 1984; Draine & Anderson 1985; Puget et al. 1985; Désert et al. 1986; Dwek 1986; Guhathakurta & Draine 1989; Aannestad 1989; Manske & Henning 1998). The only difference between the thermal approach and the statistical approach concerns the calculation of the downward transition rates $`𝐓_{lu}`$ which, in contrast to the exact-statistical treatment \[eq. (29)\], uses the notion of “grain temperature”.

We will consider two different approximations which make use of the thermal approximation, and compare results for energy distributions and infrared emission spectra.

### 6.1 “Thermal-Discrete” Approximation

In the “thermal-discrete” approximation we use the same energy bins used in the exact statistical treatment discussed in §4, and allow the same upward and downward discrete transitions. However, instead of using eq.(29) to evaluate the downward transition probabilities $`𝐓_{lu}`$, we replace

$$\frac{g_l}{g_u}\frac{\mathrm{\Delta }E_l}{\mathrm{\Delta }E_u}\frac{1}{\mathrm{exp}(h\nu /kT_u)1}$$
(35)
where $`T_u`$ is the “temperature” of the upper level. With this replacement, the downward transition rate becomes, for $`l<u1`$:555 Siebenmorgen et al. (1992) discussed both discrete heating and cooling. However, their downward transition elements were evaluated by integrating over the initial, rather than final, enthalpy interval, and hence are in error by a factor $`\mathrm{\Delta }E_u/\mathrm{\Delta }E_l`$, which is large when $`\mathrm{\Delta }E_u\mathrm{\Delta }E_l`$, as is often the case. Manske & Henning (1998) also discussed discrete cooling. In the present notation their equation (4) becomes    $`𝐓_{lu}`$ $`=`$ $`{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle _{E_uE_{l,\mathrm{max}}}^{E_uE_{l,\mathrm{min}}}}{\displaystyle \frac{E^3/\mathrm{\Delta }E_u}{\mathrm{exp}\left(E/kT_u\right)1}}C_{\mathrm{abs}}(E)𝑑E.`$   Comparison with eq. (36) reveals that Manske & Henning’s integrand is off by a factor $`(E_uE_l)/[\mathrm{\Delta }E_u\mathrm{\Delta }E_lG_{ul}(E)](E_uE_l)/\mathrm{min}(\mathrm{\Delta }E_l,\mathrm{\Delta }E_u)`$. This is a large error for the common case where the photon energy is large compared to the energy bin width, $`E\mathrm{min}(\mathrm{\Delta }E_l,\mathrm{\Delta }E_u)`$.

$$𝐓_{lu}^{\mathrm{td}}=\frac{8\pi }{h^3c^2}\frac{\mathrm{\Delta }E_l}{E_uE_l}_{W_1}^{W_4}G_{ul}(E)\frac{E^3C_{\mathrm{abs}}(E)}{\mathrm{exp}(E/kT_u)1}\left[1+\frac{h^3c^3}{8\pi E^3}u_E\right]𝑑E\mathrm{for}0<l<u1.$$
(36)
For $`l=u1`$ we include the energy-loss in “intrabin” transitions, as in eq. (31):

$`𝐓_{u1,u}^{\mathrm{td}}`$ $`=`$ $`{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{\mathrm{\Delta }E_l}{E_uE_l}}{\displaystyle _0^{W_4}}G_{u,u1}(E){\displaystyle \frac{E^3C_{\mathrm{abs}}(E)}{\mathrm{exp}(E/kT_u)1}}\left[1+{\displaystyle \frac{h^3c^3}{8\pi E^3}}u_E\right]𝑑E+`$ (37)
$`{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{1}{E_uE_{u1}}}{\displaystyle _0^{\mathrm{\Delta }E_u}}\left(1{\displaystyle \frac{E}{\mathrm{\Delta }E_u}}\right){\displaystyle \frac{E^3C_{\mathrm{abs}}(E)}{\mathrm{exp}(E/kT_u)1}}\left[1+{\displaystyle \frac{h^3c^3}{8\pi E^3}}u_E\right]𝑑E.`$
We refer to eq. (36,37) as the “thermal-discrete” approximation. The advantage of the “thermal-discrete” approximation is that one does not require explicit knowledge of the degeneracies $`g_u`$ and $`g_l`$ appearing in equation (29,31). We obtain the temperature $`T_u`$ from eq. (32).666Note that the thermal approximation (36) reduces to the statistical treatment (29) if the temperature is taken to be    $`T_u`$ $`=`$ $`{\displaystyle \frac{h\nu /k}{\mathrm{ln}\left[1+(dN/dE)|_{E_u}/(dN/dE)|_{E_uh\nu }\right]}}`$   but this definition of $`T_u`$ depends on both $`u`$ and the emitted photon energy $`h\nu `$ – only in the limit of infinite degrees of freedom does $`T_u`$ no longer depend on the value of $`h\nu `$.

As discussed in §4.1, the vibrational ground state is in a “bin” with energy $`E_0=0`$ and width $`\mathrm{\Delta }E_0=0`$, for which eq. (36) would give $`𝐓_{0u}^{\mathrm{td}}=0`$. There is really no “correct” way to approximate the energy levels as continuous when dealing with the ground state. We simply replace $`\mathrm{\Delta }E_0`$ in eq. (36) by $`\mathrm{\Delta }E_1`$.777If this seems rather arbitrary, we remind the reader that the “correct” approach to the problem is the statistical treatment given in §4. The “thermal approximation” considered here is inappropriate when discussing low degrees of excitation.  Thus we do not use eq. (36) for transitions to the ground state, and instead take888 It can be verified that if we use $`\mathrm{\Delta }E_1`$ for the width of bins 0 and 1, and $`g_0=g_1=1`$, then $`T_1=E_1/k\mathrm{ln}2`$ from eq. (32) gives the same transition rate $`𝐓_{01}`$ as the exact-statistical treatment.

$$𝐓_{0u}^{\mathrm{td}}=\frac{8\pi }{h^3c^2}\frac{\mathrm{\Delta }E_1}{\mathrm{\Delta }E_u}\frac{1}{E_u}_{E_{u,\mathrm{min}}}^{E_{u,\mathrm{max}}}\frac{E^3}{\mathrm{exp}(E/kT_u)1}\left[1+\frac{h^3c^3}{8\pi E^3}u_E\right]C_{\mathrm{abs}}(E)𝑑E.$$
(38)
### 6.2 “Thermal-Continuous” Cooling Approximation

In all of the above discussion, the grain is assumed to make discrete transitions to energy levels $`l<u`$ by emission of single photons – so-called “discrete cooling”. With discrete cooling included, the transition matrix $`𝐓_{fi}`$ is generally nonzero through the upper triangular portion $`i>f`$.

There are substantial computational advantages (see §8 below) if the cooling of the grains is approximated as continuous rather than discrete, so that the only downward transition from a level $`u`$ is to the adjacent level $`u1`$. We refer to this as the “thermal continuous” cooling approximation. We take

$`𝐓_{lu}^{\mathrm{tc}}`$ $`=`$ $`0\mathrm{if}l<u1`$ (39)
$`𝐓_{u1,u}^{\mathrm{tc}}`$ $`=`$ $`{\displaystyle \frac{1}{(E_uE_{u1})}}{\displaystyle \underset{l=0}{\overset{u1}{}}}(E_uE_l)𝐓_{lu}^{\mathrm{td}}`$ (40)
Using (36), eq.(E1), and eq. (E2), the transition rate (40) from levels $`u>2`$ can be approximated

$$𝐓_{u1,u}^{\mathrm{tc}}\frac{1}{(E_uE_{u1})}\frac{8\pi }{h^3c^2}_0^{E_u}\frac{E^3C_{\mathrm{abs}}(E)}{\mathrm{exp}(E/kT_u)1}𝑑E\mathrm{for}u>1.$$
(41)
Because of the computational advantages (see §8.2), the continuous cooling approximation has been often used in previous studies (e.g., Guhathakurta & Draine 1989). While clearly appropriate for large systems with many degrees of freedom (for which $`kTE`$ and emission of one photon does not result in a large reduction in temperature), validity of the continuous cooling approximation for grains with $`10^2`$ atoms is certainly not obvious a priori. Manske & Henning (1998) in fact argued that the continuous cooling approximation could lead to significant error. We will examine the accuracy of the continuous-cooling approximation in §10 below, where we will show that it is in fact surprisingly accurate even for grains with as few as $`30`$ atoms.

## 7 Cooling Time

The absorption cross section $`C_{\mathrm{abs}}(E)`$ determines both the rate of heating by starlight and the rate of cooling by infrared emission. It is useful to define a “radiative cooling time”

$`\tau _{\mathrm{rad}}(E_u)`$ $``$ $`{\displaystyle \frac{E_u}{_{l<u}𝐓_{lu}(E_uE_l)}}`$ (42)
$``$ $`\left[{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{1}{E_u}}{\displaystyle \frac{1}{(dN/dE)_u}}{\displaystyle _0^{E_u}}x^3C_{\mathrm{abs}}(x)\left[1+{\displaystyle \frac{h^3c^3}{8\pi x^3}}u_x\right]\left({\displaystyle \frac{dN}{dE}}\right)_{E_ux}𝑑x\right]^1`$ (43)
where eq. (E2) has been used to obtain eq. (43). In the thermal-discrete approximation, we see from eq. (41) that

$$\tau _{\mathrm{rad}}^{\mathrm{td}}(E_u)\left[\frac{1}{E_u}\frac{8\pi }{h^3c^2}_0^{E_u}\frac{E^3C_{\mathrm{abs}}(E)}{\mathrm{exp}(E/kT_u)1}𝑑E\right]^1\mathrm{for}u>1.$$
(44)
Cooling times for PAH and silicate grains are plotted in Figs. 9 and 10. The interstellar starlight radiation field (ISRF) of Mathis, Mezger & Panagia (1983, hereafter MMP) has been assumed for $`u_E`$; for this very dilute radiation field the correction for stimulated emission is entirely negligible ($`h^3c^3u_E/8\pi E^31`$ for $`E/hc10\mathrm{cm}^1`$).

We show $`\tau _{\mathrm{rad}}`$ evaluated using the “exact-statistical” method, as well as using the “thermal-discrete” approximation ($`\tau _{\mathrm{rad}}`$ for the continuous cooling approximation is, by construction, identical to that for the “thermal-discrete” approximation).

As illustrated in Figs. 9-10, $`\tau _{\mathrm{rad}}(E)`$ increases with grain size; for a given grain size, $`\tau _{\mathrm{rad}}`$ decreases rapidly with increasing vibrational energy $`E`$. The importance of single photon heating is apparent: very small grains emit most of their energy very rapidly (within seconds) following absorption of an energetic photon and spend most of their time in or close to the ground vibrational state.

The substantial difference between the results of the statistical model and the thermal models at low $`E`$ is mainly due to the treatment of transitions to the ground state, which are problematic in the thermal approximation (see §6.1). The jagged behavior of $`\tau _{\mathrm{rad}}`$ in the exact-statistical model result from the discreteness of the lowest energy states (see eq. 4) which affects the degeneracies (see eq. 29). The discontinuity in $`\tau _{\mathrm{rad}}`$ computed in the thermal approximation (Figures 9,10) is due to the discontinuity in the prescription (32) for the grain “vibrational temperature” at $`E=\mathrm{}\omega _{20}`$ (see Figs. 7,8). Note that with the present prescription for the temperature, our estimates for $`\tau _{\mathrm{rad}}`$ in the thermal approximation are in reasonable agreement with $`\tau _{\mathrm{rad}}`$ calculated in the “exact-statistical” method at the lowest energies. For C4000H1000 in Figure 9, and $`a=20\mathrm{\AA }`$ in Figure 10, the discontinuity in $`\tau _{\mathrm{rad}}`$ at $`3\times 10^4\mathrm{cm}^1`$ is due to a transition from direct computation of the density of states to use of a continuous Debye model.

The mean time $`\tau _{\mathrm{abs}}`$ between photon absorptions for a grain is given by

$$\tau _{\mathrm{abs}}^1_0^{\mathrm{}}C_{\mathrm{abs}}(h\nu )\frac{cu_\nu }{h\nu }𝑑\nu ,$$
(45)
In Figure 11 we show $`\tau _{\mathrm{abs}}`$ for grains in the MMP ISRF. The mean absorbed photon energy $`h\nu _{\mathrm{abs}}`$ is

$$h\nu _{\mathrm{abs}}\tau _{\mathrm{abs}}_0^{\mathrm{}}C_{\mathrm{abs}}(h\nu )cu_\nu 𝑑\nu .$$
(46)
Ionized PAHs with $`N_\mathrm{C}<10^4`$ in the MMP radiation field have $`h\nu _{\mathrm{abs}}5.2\mathrm{eV}`$; neutral PAHs have a somewhat higher value of $`h\nu _{\mathrm{abs}}`$ due to their different absorption cross section (see Figure 6).

Also plotted in Figure 11 is the cooling time $`\tau _{\mathrm{rad}}`$ for grains containing 5.2 eV of thermal energy. The critical size $`N_{\mathrm{sph}}`$ (and radius $`a_{\mathrm{sph}}`$) for single-photon heating is determined by the condition $`\tau _{\mathrm{rad}}(N_{\mathrm{sph}})=\tau _{\mathrm{abs}}(N_{\mathrm{sph}})`$. For carbonaceous grains, we find

$$N_{\mathrm{C},\mathrm{sph}}5\times 10^5\left(\frac{u}{u_{\mathrm{MMP}}}\right)^{0.60},a_{\mathrm{sph}}100\mathrm{\AA }\left(\frac{u}{u_{\mathrm{MMP}}}\right)^{0.20}.$$
(47)
For grains with $`N<N_{\mathrm{sph}}`$ atoms, single-photon heating effects are very important: the radiative cooling time is short compared to the interval between photon absorptions, and the grains spend most of the time between photon absorptions at energies $`Eh\nu _{\mathrm{abs}}`$. Grains with $`NN_{\mathrm{sph}}`$ heat up until their energy content is large compared to the energies of starlight photons, and individual photon absorptions produce small fractional changes in the thermal energy content. For such large grains, one may therefore estimate the average energy content $`\overline{E}`$ by requiring a steady balance between cooling and heating:

$$\frac{\overline{E}}{\tau _{\mathrm{rad}}(\overline{E})}=𝑑\nu cu_\nu C_{abs}(h\nu ).$$
(48)
## 8 Solution Method

### 8.1 Discrete Cooling

If we define the diagonal elements of $`𝐓`$ to be

$$𝐓_{ii}=\underset{ji}{}𝐓_{ji}$$
(49)
then equation (14) becomes

$$\underset{j=0}{\overset{M}{}}𝐓_{ij}P_j=0\mathrm{for}i=0,\mathrm{},M.$$
(50)
Combining this with the normalization condition $`_{j=0}^MP_j=1`$, we obtain a set of $`M`$ linear equations for the first $`M`$ elements of $`P_j`$:

$$\underset{j=0}{\overset{M1}{}}\left(𝐓_{ij}𝐓_{iM}\right)P_j=𝐓_{iM}\mathrm{for}i=0,\mathrm{},M1.$$
(51)
The remaining undetermined element $`P_M`$ is obtained by

$$P_M=\left(𝐓_{MM}\right)^1\underset{j=0}{\overset{M1}{}}𝐓_{Mj}P_j.$$
(52)
For small $`M`$ we can directly solve eq. (51) for $`P_j`$ using Gaussian elimination. When $`M10^2`$, however, direct solution, requiring $`O(M^3)`$ operations, is both slow and numerically ill-behaved. We therefore resort to iterative techniques, and have tried both the bi-conjugate gradient (BiCG) method (Press et al. 1992), and the stabilized bi-conjugate gradient (BiCGstab) method (see Sleijpen & Fokkema 1993 and references therein). We found the BiCG method to be quite efficient even without preconditioning, whereas the convergence rate of the BiCGstab method is dependent on the choice of preconditioner.

Ideally, the adopted highest energy bin $`E_M`$ for a given grain would be as large as possible, but large values of $`M`$ entail heavy memory and cpu requirements. In our calculations, we first set $`E_M=13.6\mathrm{eV}`$ and solve eq. (51) for $`P_j`$. If $`P_M<10^{14}`$ then we use this $`E_M`$ and solution $`P_j`$; otherwise we increase $`E_M`$ by a factor 1.5 and solve again, repeating until we obtain a solution with $`P_M<10^{14}`$.

Small ($`a20`$Å) grains exposed to the interstellar starlight background are able to cool completely betweeen photon absorptions, and $`E_M=13.6\mathrm{eV}`$ is sufficient. For larger grains, the interval between photon absorptions is shorter and the cooling time at given $`E`$ is longer (see Figs. 9 and 10), so that the grain generally does not cool off to the ground state between photon absorptions, and $`E_M`$ must be higher than 13.6eV. For a grain with $`a100`$Å only modest cooling takes place between photon absorptions and $`P_j`$ becomes strongly peaked around a steady state “equilibrium temperature”.

### 8.2 Continuous Cooling

For the continous cooling approximation it is straightforward to solve for the steady-state solution vector $`P_j`$. As shown by Guhathakurta & Draine (1989), if we define $`X_jP_j/P_0`$, then $`X_0=1`$ and the remaining $`X_j`$ may be obtained recursively:

$$X_j=\frac{1}{𝐓_{j1,j}^{\mathrm{tc}}}\underset{i=0}{\overset{j1}{}}B_{j,i}X_i,B_{j,i}\underset{u=j}{\overset{M}{}}𝐓_{u,i}^{\mathrm{tc}}$$
(53)
The recursion relation $`B_{j,i}=𝐓_{j,i}+B_{j+1,i}`$ can be used to generate the $`B_{j,i}`$ in $`O(M^2)`$ operations, and the $`X_j`$ can then be computed in $`O(M^2)`$ operations. Once the $`X_j`$ are determined for $`j=0,\mathrm{},M`$, one need only renormalize:

$$P_j=\frac{X_j}{_{j=0}^MX_j}.$$
(54)
Since $`M100`$ is required to have suitably narrow energy bins (we often use $`M500`$), this $`O(M)`$ procedure is much faster than the BiCG iterative method required if discrete cooling is included.

## 9 Infrared Emission Spectrum

Having obtained the steady-state energy distribution $`P_j`$, we wish to calculate the resulting infrared emission spectrum. Let $`F_\nu d\nu `$ be the time-averaged power per steradian per grain radiated into frequency interval $`[\nu ,\nu +d\nu ]`$. Then the “exact statistical” treatment gives

$$\nu F_\nu =\frac{2h\nu ^4}{c^2}C_{\mathrm{abs}}(\nu )\left[\underset{u}{}^{}P_u\underset{l=0}{\overset{u1}{}}\frac{g_l}{g_u}\mathrm{\Delta }E_uG_{ul}(h\nu )+\underset{u}{}^{\prime \prime }P_u\left(1\frac{h\nu }{\mathrm{\Delta }E_u}\right)\right]\left[1+\frac{\lambda ^3}{8\pi }u_E\right]$$
(55)
The single-primed summation – the contribution of “interbin” transitions to the emission – is limited to levels $`u`$ with $`E_{u,\mathrm{max}}>h\nu `$.999Note that only a few terms from $`_l^{}`$ have $`G_{ul}0`$ – energy levels $`klm`$ where $`E_{u,\mathrm{min}}h\nu [E_{k,\mathrm{min}},E_{k,\mathrm{max}}]`$ and $`E_{u,\mathrm{max}}h\nu [E_{m,\mathrm{min}},E_{m,\mathrm{max}}]`$. Therefore eq. (55) requires only $`O(M)`$ \[rather than $`O(M^2)`$\] operations per frequency $`\nu `$.  The double-primed summation – the contribution of “intrabin” transitions – is limited to levels with $`\mathrm{\Delta }E_u>h\nu `$.

In the thermal approximation, equation (55) becomes (see Appendix D)

$$\nu F_\nu \frac{2h\nu ^4}{c^2}C_{\mathrm{abs}}(\nu )\left[\underset{u}{}^{}\frac{P_u}{\mathrm{exp}(h\nu /kT_u)1}\right]\left[1+\frac{\lambda ^3}{8\pi }u_E\right],$$
(56)
where the single-primed summation is now limited to levels $`u`$ with $`E_u>h\nu `$.

## 10 Results

We have calculated the (vibrational) energy probability distribution $`P_j`$ for PAHs excited by the MMP ISRF using the above-described methods: the exact-statistical model (§4); the thermal-discrete model (§6.1); and the thermal-continuous model (§6.2).

In Figure 12 we present the cumulative energy probability distributions for selected PAHs obtained from the exact-statistical model, the thermal-discrete model, and the thermal-continuous model. Figure 13 shows the differential probability distributions as a function of grain energy $`E`$ for these same grain sizes.

The inset in Figure 13 shows the probability of being in the vibrational ground state as a function of grain size, calculated for the exact-statistical model and the thermal-discrete model. As already illustrated in Figure 12 and Figure 13, the probability of being in the ground state is very large for small grains: for example, for the MMP radiation field, grains with $`N4000`$ spend most of their time at $`E=0`$ \[the exact-statistical model gives $`P_0=0.975`$ for $`N=400`$, and $`P_0=0.316`$ for $`N=4000`$\]. The sharp drop at 13.6 eV ($`E/hc=1.1\times 10^5`$ cm-1) is due to the radiation field cutoff at 912 Å and to the fact that multiphoton events are rare. The jagged structure at low energies in Figure 13 is due to the quantization of the lowest energy states. These structures are less prominent in the thermal models because the statistical model involves the degeneracies of each “bin” (see eq. 30) which depend on the detailed distribution of energy states (see, e.g., Figure 4).

The energy distributions $`P(E)`$ found using the thermal-discrete model and the thermal-continuous model are both in good overall agreement with the results of the exact-statistical calculation.

It is not surprising that all three models are in good agreement at high energies since the transition rate from a very high energy state to the ground state is negligibly small, and the number of excited degrees of freedom is large enough that a thermal treatment is appropriate. The agreement becomes better with increasing grain size because the lowest vibrational states (i.e., the first few excited states) become closer and approach a “continuum”, and the gap between the ground state and the first excited state becomes smaller.

The resulting IR emission spectra are displayed in Figs. 14 and 15. As expected (from the energy distribution functions plotted in Figs. 12 and 13), the IR spectra of the thermal-discrete model are almost identical to those of the exact-statistical model, even for grains as small as $`N_\mathrm{C}=24`$. The thermal-continuous cooling model also results in spectra which are very close to those computed using the exact-statistical model. The discrepancies are mainly at long wavelengths; this can be understood from the energy probability functions and would not affect the interstellar IR emission spectrum modeling.

The statistical model spectra are generally very close to those of the thermal models (Figs. 14 and 15), except for the sawtooth features at long wavelengths. These sawtooth features are due to our treatment of transitions from the lower excited energy bins to the ground state and first few excited states. However, we stress that the overall spectra are quite similar for all three models – the differences involve only a negligible fraction of the total emission for any given grain size, and thus would be unimportant when modeling the overall interstellar IR emission spectrum.

## 11 Discussion

### 11.1 Band Ratios

From Figs. 14 and 15 it is evident that the relative strength of the different PAH emission bands (3.3, 6.2, 7.7, 8.6, 11.3$`\mathrm{\mu m}`$) is a strong function of the PAH size: small PAHs radiate strongly at 6.2 and 7.7$`\mathrm{\mu m}`$, while larger PAHs emit most of their power at increasingly long wavelengths.

The 6.2, 7.7, and 11.3$`\mathrm{\mu m}`$ features are prominent in many astronomical spectra. Let $`P(\lambda _0)`$ be the power radiated in the feature with central wavelength $`\lambda _0`$: for a Drude profile, $`P(\lambda _0)=(\pi /2)(\mathrm{\Delta }P_\lambda )_0(\mathrm{FWHM})_{\lambda _0}`$, where $`(\mathrm{\Delta }P_\lambda )_0`$ is the peak contribution of the feature to the emission, and $`(\mathrm{FWHM})_{\lambda _0}`$ is the full-width at half-maximum. In Figure 16 we show how the relative strengths of these emission bands vary depending upon the size and charge state of the PAHs, and upon the starlight intensity $`u_E`$, characterized by $`\chi `$ or $`G_0`$.101010 $`\chi `$ is the intensity at 1000$`\mathrm{\AA }`$ relative to the estimate of Habing (1968); $`G_0`$ is the ratio of the 6–13.6 eV energy density to the estimate of Habing (1968). The MMP spectrum has $`G_0=0.923\chi `$; a $`3\times 10^4\mathrm{K}`$ blackbody has $`G_0=0.602\chi `$.  Neutral PAHs have large ratios of $`P(11.3\mathrm{\mu m})/P(7.7\mathrm{\mu m})`$, while PAH ions have much smaller values. For both neutrals and ions, there is a regular progression of $`P(6.2\mathrm{\mu m})/P(7.7\mathrm{\mu m})`$ to smaller values with increasing $`N_\mathrm{C}`$. For a given PAH neutral or ion, the emission band ratios are essentially independent of $`\chi `$ for $`\chi <3\times 10^4(10^3/N_\mathrm{C})^{5/3}`$ \[see eq. (47)\]; for larger values of $`\chi `$, multiphoton heating effects begin to be evident. For $`\chi =10^6`$, multiphoton effects are evident even for $`N_\mathrm{C}=10^2`$ \[consistent with eq. (47)\].

In Figure 16 we plot the observed ratios of $`P(11.3\mathrm{\mu m})/P(7.7\mathrm{\mu m})`$ and $`P(6.2\mathrm{\mu m})/P(7.7\mathrm{\mu m})`$ for (1) the diffuse ISM (Onaka et al. 1996); (2) the average spectrum for a sample of 28 spiral galaxies (Helou et al. 2000); (3) the reflection nebula at the surface of the $`\rho `$ Oph molecular cloud (Boulanger et al. 1996); (4) the reflection nebula vdB 133, illuminated by F5Iab and B7II stars (Uchida et al. 1998); (5) the reflection nebula/photodissociation region (PDR) NGC 7023, illuminated by a B2.5V star (Cesarsky et al. 1996a); (6) the M17 PDR (Cesarsky et al. 1996b); (7) the starburst galaxies M82 and NGC 253 (Stürm et al. 2000); (8) the Seyfert 2 galaxy in Circinus (Stürm et al. 2000); and (9) the quiescent molecular cloud SMC B1#1 in the Small Magellanic Cloud (Reach et al. 2000). With the exception of SMC B1#1 (which we discuss below), the observed band ratios fall between the neutral and ionized PAH tracks in Figure 16, suggesting that the observed spectra can be reproduced by an appropriate mixture of PAH sizes, with an appropriate ionized fraction. The observed band ratios for the Milky Way ISM, $`\rho `$ Oph, 28 normal galaxies, and the reflection nebulae vdB 133 and NGC 7023 could be reproduced with a relatively low ionized fraction, but larger ionized fractions are required for the M17 PDR, M82, NGC 253, and the Circinus galaxy.

### 11.2 SMC B1#1

The quiescent molecular cloud SMC B1#1 observed by Reach et al. (2000) falls outside the region bounded by the neutral and ionized tracks in Figure 16. Reach et al. point out that SMC B1#1 has a 11.3/7.7 band ratio which differs significantly from that for Milky Way objects, including the $`\rho `$ Oph reflection nebula (which one might expect SMC B1#1 to resemble). Reach et al. suggest that the enhanced emission in the 11.3$`\mathrm{\mu m}`$ C-H bending mode might be due to increased hydrogenation of the PAHs in the SMC, where the gas phase H/C ratio is a factor of 10 higher than in the Milky Way. Given that Milky Way PAHs with $``$25 C atoms are expected to be essentially fully hydrogenated (Tielens et al. 1987; Allamandola, Tielens, & Barker 1989; Allain, Leach, & Sedlmayr 1996), it is not clear that greater hydrogenation would be expected in the SMC. As discussed in Paper II, our model calculations assume that interstellar neutral and ionized PAHs have the 11.3/7.7$`\mathrm{\mu m}`$ band strength ratio reduced by a factor of 2 relative to PAHs which have been studied in the lab – enhancement of the 7.7$`\mathrm{\mu m}`$ feature appears to be necessary to reproduce the low 11.3/7.7$`\mathrm{\mu m}`$ band ratios observed in objects such as NGC 253 or M82, and (with the possible exception of SMC B1#1) the reduced 11.3/7.7$`\mathrm{\mu m}`$ band strength ratio appears to be consistent with astronomical observations.

While SMC B1#1 appears to fall “outside” the “allowed” region in Figures 1618, we show elsewhere (Li & Draine 2001b) that the present model can reproduce the observed line ratios by adding emission from a population of larger PAHs – with $`N_\mathrm{C}3000`$ – which radiate strongly in the 11–13$`\mathrm{\mu m}`$ region, but produce little emission in the 6–8$`\mathrm{\mu m}`$ region. These larger PAHs would radiate strongly at 15–50$`\mathrm{\mu m}`$; SIRTF will be able to test this hypothesis.

### 11.3 PAH Ionization

The fraction of PAHs which are charged will depend upon the electron density $`n_e`$, temperature $`T`$, and the starlight intensity. We consider three different charging environments:

* “Cold neutral medium” (CNM; $`n_e=0.045\mathrm{cm}^3`$, $`T=100\mathrm{K}`$, and the MMP ISRF spectrum with $`\chi =1.23`$, $`G_0=1.14`$) and with both radiation field and gas density increased by $`10^2`$;
* “Warm ionized medium” (WIM; $`n_e=0.1\mathrm{cm}^3`$, $`T=8000\mathrm{K}`$, and the MMP ISRF with $`\chi =1.23`$), and with both radiation field and gas density increased by $`10^2`$;
* Photodissociation region conditions (PDR; $`n_e/n_\mathrm{H}=3\times 10^4`$, $`T=10^3\mathrm{K}`$, and a $`3\times 10^4\mathrm{K}`$ blackbody cut off at 912$`\mathrm{\AA }`$ with $`\chi /n_\mathrm{H}=0.1\mathrm{cm}^3`$) for $`\chi =10^2`$, $`10^4`$, $`10^5`$, and $`10^6`$. While PDRs come in a range of conditions, recent observations and modeling suggest that at least some PDRs include an extended region with $`T10^3\mathrm{K}`$, and with $`\chi /n_\mathrm{H}0.1\mathrm{cm}^3`$, as in a recent model for the NGC 2023 PDR (Draine & Bertoldi 2000).

For each set of conditions, we consider PAHs with $`N_\mathrm{C}16`$ C atoms, and include emission from both neutral and ionized forms, with the neutral fraction calculated assuming a balance between collisional charging by electrons and photodetachment and photoionization (Weingartner & Draine 2001b). The charging of small grains is largely a function of the quantity $`\chi \sqrt{T}/n_e`$ (Bakes & Tielens 1994; Weingartner & Draine 2001b). For the CNM, WIM, and PDR conditions considered here, $`\chi \sqrt{T}/n_e`$ = 270, 1100, and $`10^4`$ $`\mathrm{cm}^3\mathrm{K}^{1/2}`$, so we expect grains of a given size to be increasingly positively charged as we move from CNM to WIM to PDR.

The “diffuse ISM” and “28 Galaxies” points fall between the CNM and WIM tracks in Figure 16, implying that one could reproduce the overall emission from normal spiral galaxies by a weighted sum of interstellar conditions ranging from CNM to WIM, with a mixture of grain sizes such that the 6.2 – 11.3$`\mathrm{\mu m}`$ emission is dominated by grains with $`20N_\mathrm{C}500`$.

The observed flux ratios for NGC 7023 and M17 fall between the “WIM” and “PDR” tracks in Figure 17. The observed flux ratios for the starburst galaxies M82 and NGC 253, as well as the Circinus Seyfert 2 galaxy, fall close to the “PDR” tracks, consistent with the expectation that the spectra of these galaxies includes a substantial contribution from PDRs.

It is important to note that the “PDR” models in Figure 17 have charging conditions such that a substantial fraction of the small grains are still neutral: 22% of the $`N_\mathrm{C}=100`$ grains are neutral. Larger charged fractions – and smaller values of $`P(11.3\mathrm{\mu m})/P(7.7\mathrm{\mu m})`$ – can be achieved by increasing the value of $`\chi /n_\mathrm{H}`$. For example, if $`\chi /n_\mathrm{H}`$ is increased by a factor 2, only 10% of the $`N_\mathrm{C}=100`$ grains are neutral.

In Figure 18 we show the emission expected from grains of various sizes, with 43% in the CNM, 14% in the WIM, and 43% in the WNM ($`n_e=0.03\mathrm{cm}^3`$, $`T=6000\mathrm{K}`$, MMP ISRF, $`\chi \sqrt{T}/n_e=3200`$). The shaded region shows the flux ratios which could be obtained by grain mixtures containing only grains $`20<N_\mathrm{C}<1350`$ (from Figure 15 it is apparent that grains with $`N_\mathrm{C}>10^3`$ are not likely to contribute significant emission at $`\lambda 11.3\mathrm{\mu m}`$). The observed flux ratios for the diffuse ISM in the Milky Way, and for the average of 28 spiral galaxies, both fall in the shaded region, suggesting that with an appropriate size distribution, one might be able to reproduce the observed emission. The question of size distributions will be addressed in a future paper (Li & Draine 2001a), where it will be shown that the observed infrared emission can be approximately reproduced by grain size distributions similar to those invoked to explain the observed microwave emission from interstellar dust (Draine & Lazarian 1998a,b).

## 12 Summary

We have presented a method for calculating the infrared emission from dust grains, including very small PAH molecules or silicate clusters. The principal results of this paper are the following:

1. The vibrational mode spectrum and vibrational density of states of very small grains are discussed in §2, where we explicitly consider both polycyclic aromatic hydrocarbons and silicate grains.
2. We use the adopted vibrational mode spectrum to calculate the specific heat of carbonaceous and silicate grains; the results are in agreement with experimental specific heats for graphite and amorphous silicates (Figure 2).
3. Using estimates for the absorption cross sections of carbonaceous and silicate grains, we obtain (in §4) radiative transition rates for emission from vibrationally-excited states of carbonaceous and silicate grains. We refer to this as the “exact-statistical” treatment.
4. We discuss (in §6) how the “thermal approximation” can be used to estimate rates for spontaneous emission. With a judicious estimate for the “temperature”, the thermal approximation provides transition rates which are close to those obtained from the “exact-statistical” treatment. Therefore level populations estimated using the thermal approximation are close to the actual level populations given by the “exact-statistical” treatment.
5. We discuss the “continuous cooling” approximation for modeling the deexcitation of the grain. We show that level populations and emission spectra computed using this approximation are sufficiently accurate for most astrophysical applications, even for grains containing as few as $`30`$ atoms.
6. Methods for numerical solution for the vibrational energy distribution are presented in §8, and the numerical advantages of the continuous cooling approximation are stressed.
7. Relative strengths of the 6.2, 7.7, and 11.3$`\mathrm{\mu m}`$ features depend on the grain size, on the charging conditions, and on the starlight intensity. Observed 6.2/7.7 and 11.3/7.7 band ratios for the diffuse ISM in the Milky Way, for the $`\rho `$ Oph and vdB133 reflection nebulae, and for the 28 “normal” galaxy sample of Helou et al. (2000), appear consistent with grains in the diffuse ISM heated by the average starlight background. The band ratios observed for the starburst galaxies M82 and NGC 253, for the Seyfert 2 “Circinus galaxy”, for the M17 PDR, and for the NGC 7023 PDR appear to require conditions close to our “PDR” conditions, with $`\chi /n_\mathrm{H}0.050.1\mathrm{cm}^3`$.
8. The unusual band strengths observed for the quiescent molecular cloud SMC B1#1 in the SMC (Reach et al. 2000) cannot be reproduced by a mixture of small $`N_\mathrm{C}10^3`$ PAHs. Either the PAHs in this cloud differ from Milky Way PAHs in their properties, or have a different size distribution.

The present paper has concentrated on the statistical mechanics of the stochastic heating process, and on the resulting emission spectra as a function of grain size. In a separate paper (Li & Draine 2001a) we calculate the emission from grain models with realistic size distributions, and show that the observed infrared emission from interstellar dust can be reproduced by a grain model with $``$15% of the interstellar carbon in PAH particles containing less than $`10^3`$ C atoms.

We thank L.J. Allamandola, J.M. Greenberg, and J.C. Weingartner for helpful discussions. We thank D.M. Hudgins and L.J. Allamandola for providing us with their PAH database. We thank R.H. Lupton for availability of the SM plotting package. This research was supported in part by NSF grant AST-9619429, and NASA grant NAG5-7030.

## Appendix A PAH Geometry

Small PAHs are expected to be planar, with $`DN_\mathrm{C}^{1/2}`$, but large PAHs are expected to be approximately spherical, with $`DN_\mathrm{C}^{1/3}`$. The size at which the transition from planar to spherical geometry occurs is uncertain. We will assume that PAHs have a planar geometry up to $`N_\mathrm{C}=54`$ (e.g., the pericondensed species circumcoronene C54H18). We assume that beyond $`N_\mathrm{C}=54`$ PAHs consist of multiple sheets. Beyond $`N_\mathrm{C}=102`$ (e.g., C54H18 plus one coronene C24H12 on either side) we assume spherical geometry, with $`DN_\mathrm{C}^{1/3}`$.

## Appendix B Energy Bins

We specify the energy bins so that the total number $`M`$ of bins is manageable ($`300M1500`$), the number of vibrational states per bin is at least 2 for each of bins 3–10, the bin widths vary smoothly for $`j>10`$, and the last bin is at a high enough energy that its population $`P_M<10^{14}`$ (see §8.1 below).

Let $`\omega _j,j=1,\mathrm{},3(N_a2)`$ be the fundamental vibrational modes (see §2), in order of increasing frequency. Bin 0 is the ground state: $`E_{0,\mathrm{min}}=E_{0,\mathrm{max}}=E_0=0`$, to which we assign a “width” $`\mathrm{\Delta }E_0=0`$. Bins 1 and 2 contain the first two vibrational modes, and bins 3 through $`10`$ each contain two normal modes:

$`E_{1,\mathrm{min}}`$ $`=`$ $`{\displaystyle \frac{3}{2}}\mathrm{}\omega _1{\displaystyle \frac{1}{2}}\mathrm{}\omega _2`$ (B1)
$`E_{j,\mathrm{max}}=E_{j+1,\mathrm{min}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\mathrm{}\omega _j+\mathrm{}\omega _{j+1}\right)\mathrm{for}j=1,2`$ (B2)
$`E_{j,\mathrm{max}}=E_{j+1,\mathrm{min}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\mathrm{}\omega _{2j2}+\mathrm{}\omega _{2j1}\right)\mathrm{for}3j10`$ (B3)
For $`j=11,\mathrm{},K`$ we use uniform intervals $`\mathrm{\Delta }E`$, and for $`jK`$ uniform intervals in $`\mathrm{ln}E`$:

$$E_{j,\mathrm{max}}=E_{10,\mathrm{max}}+(j10)(E_{10,\mathrm{max}}E_{9,\mathrm{max}})\mathrm{for}10jK$$
(B4)
$$E_{j,\mathrm{max}}=E_{K,\mathrm{max}}\left(\frac{E_{K,\mathrm{max}}}{E_{K1,\mathrm{max}}}\right)^{jK}\mathrm{for}KjM,$$
(B5)
$`K`$ is chosen so that the total number of bins $`500`$.

## Appendix C Transition Matrix

Consider energy “bins” $`u`$ and $`l`$, with $`E_u>E_l`$. Suppose bin $`l`$ consists of many subbins $`\alpha =1,\mathrm{},N_\alpha `$, each of width $`\mathrm{\Delta }E_\alpha =\mathrm{\Delta }E_l/N_\alpha `$, and bin $`u`$ consists of many subbinss $`\beta =1,\mathrm{},N_\beta `$, each with width $`\mathrm{\Delta }E_\beta =\mathrm{\Delta }E_u/N_\beta `$. The upward Einstein $`A`$ coefficient from subbin $`\alpha `$ to subbin $`\beta `$ is

$$A_{\beta \alpha }=\frac{8\pi \nu _{\alpha \beta }^2}{hc^2}\mathrm{\Delta }E_\beta C_{\mathrm{abs}}(\nu _{\alpha \beta }),$$
(C1)
and the spontaneous decay rate from subbin $`\beta `$ to subbin $`\alpha `$ is

$$A_{\alpha \beta }=\frac{(dN/dE)_\alpha }{(dN/dE)_\beta }\frac{8\pi \nu _{\alpha \beta }^2}{hc^2}\mathrm{\Delta }E_\alpha C_{\mathrm{abs}}(\nu _{\alpha \beta }),$$
(C2)
where $`\nu _{\alpha \beta }(E_\beta E_\alpha )/h`$, and $`dN/dE`$ is the vibrational density of states. We now assume that $`(dN/dE)_\alpha (dN/dE)_l`$, and $`(dN/dE)_\beta (dN/dE)_u`$.

To obtain the effective bin-to-bin transition rate $`𝐓_{ul}`$, we equate the total rate of energy absorption to $`𝐓_{ul}(E_uE_l)`$. Thus, averaged over subbins $`\alpha `$, the transition rate from $`lu`$ is

$`𝐓_{ul}`$ $`=`$ $`{\displaystyle \frac{1}{E_uE_l}}{\displaystyle \frac{c}{N_\alpha }}{\displaystyle \underset{\alpha }{}}{\displaystyle \underset{\beta }{}}C_{\mathrm{abs}}(\nu )u_E\mathrm{\Delta }E_\beta `$ (C3)
$``$ $`{\displaystyle \frac{c}{\mathrm{\Delta }E_l}}{\displaystyle \frac{1}{E_uE_l}}{\displaystyle _{E_{l,\mathrm{min}}}^{E_{l,\mathrm{max}}}}𝑑x{\displaystyle _{E_{u,\mathrm{min}}}^{E_{u,\mathrm{max}}}}𝑑y\left[C_{\mathrm{abs}}(E)u_E\right]_{E=yx}.`$ (C4)
The integrand depends only on $`E=yx`$. One can show that

$$𝐓_{ul}=\frac{c\mathrm{\Delta }E_u}{E_uE_l}_{W_1}^{W_4}𝑑EG_{ul}(E)C_{\mathrm{abs}}(E)u_E,$$
(C5)
where $`W_1`$, $`W_4`$, and $`G_{ul}(E)`$ are defined in equations (19-24).

Similarly, averaged over subbins $`\beta `$, the energy-weighted transition rate from $`ul`$ is

$`𝐓_{lu}`$ $`=`$ $`{\displaystyle \frac{1}{N_\beta }}{\displaystyle \frac{1}{E_uE_l}}{\displaystyle \underset{\beta }{}}{\displaystyle \underset{\alpha }{}}A_{\alpha \beta }(E_\beta E_\alpha )`$ (C6)
$`=`$ $`{\displaystyle \frac{(dN/dE)_l}{(dN/dE)_u}}{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{\mathrm{\Delta }E_\beta }{E_uE_l}}{\displaystyle \underset{\beta }{}}{\displaystyle \underset{\alpha }{}}E^3C_{\mathrm{abs}}(E)`$ (C7)
$``$ $`{\displaystyle \frac{(dN/dE)_l}{(dN/dE)_u}}{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{\mathrm{\Delta }E_l}{\mathrm{\Delta }E_u}}{\displaystyle \frac{1}{E_uE_l}}{\displaystyle _{W_1}^{W_4}}𝑑EG_{ul}(E)E^3C_{\mathrm{abs}}(E)`$ (C8)
$`=`$ $`{\displaystyle \frac{g_l}{g_u}}{\displaystyle \frac{8\pi }{h^3c^2}}{\displaystyle \frac{1}{E_uE_l}}{\displaystyle _{W_1}^{W_4}}𝑑EG_{ul}(E)E^3C_{\mathrm{abs}}(E).`$ (C9)
In the case of transitions to adjacent bins, we augment $`𝐓_{u,u1}`$ by the rate of energy absorption in “intrabin” transitions divided by $`E_uE_{u1}`$, and $`𝐓_{u1,u}`$ by the power radiated in intrabin transition, divided by $`E_uE_{u1}`$ (see eq. 28, 31).

## Appendix D Emission Spectrum

$`F_\nu `$, the power radiated per grain per unit solid angle per unit frequency, can be evaluated by noting from eq. (29) that the contribution from transitions $`ul`$ is just

$$\delta F_\nu d\nu =P_u\frac{g_l}{g_u}\frac{2}{h^3c^2}G_{ul}(E)E^3C_{\mathrm{abs}}(E)\left[1+\frac{\lambda ^3}{8\pi }u_E\right]dE.$$
(D1)
If the width of bin $`u`$ exceeds $`h\nu `$, then there is an additional contribution from intrabin transitions originating in the fraction $`(1h\nu /\mathrm{\Delta }E)`$ of the subbins which are more than $`h\nu `$ above $`E_{u,\mathrm{min}}`$. Using eq. (C2) with $`(dN/dE)_\alpha (dN/dE)_\beta `$, the intrabin contribution to $`F_\nu d\nu `$ is

$`\delta F_\nu d\nu `$ $`=`$ $`P_u\left(1{\displaystyle \frac{h\nu }{\mathrm{\Delta }E_u}}\right){\displaystyle \frac{E}{4\pi }}{\displaystyle \frac{dA}{dE}}\left[1+{\displaystyle \frac{\lambda ^3}{8\pi }}u_E\right]dE`$ (D2)
$`=`$ $`P_u\left(1{\displaystyle \frac{h\nu }{\mathrm{\Delta }E_u}}\right){\displaystyle \frac{E}{4\pi }}{\displaystyle \frac{8\pi \nu ^2}{hc^2}}C_{\mathrm{abs}}(E)\left[1+{\displaystyle \frac{\lambda ^3}{8\pi }}u_E\right]dE,`$ (D3)
where $`dA/dE`$ is obtained by dividing eq. (C2) by $`\mathrm{\Delta }E_\alpha `$. With $`E=h\nu `$, it is straightforward to use eq. (D1-D3) to obtain eq. (55).

## Appendix E Lemmas

From the definition (19) of $`G_{ul}(h\nu )`$ one can show that

$$\underset{l=0}{\overset{u1}{}}\mathrm{\Delta }E_u\mathrm{\Delta }E_lG_{ul}(h\nu )=\mathrm{min}(h\nu ,\mathrm{\Delta }E_u)\mathrm{for}h\nu <E_{u,\mathrm{min}}E_{1,\mathrm{min}},$$
(E1)
$$\underset{l=0}{\overset{u1}{}}\mathrm{\Delta }E_l(E_uE_l)G_{ul}(E)E,$$
(E2)
which are used to obtain eq. (40) and (43).


format: text

text:
# Effect of Grain Geometry on Angle of Repose and Dynamics

\[

## Abstract

With simulation and experiment we examine the effect of grain anisotropy on the angle of repose in a two-dimensional geometry. We find that the structure of the granular packing can be directly controlled by changing the aspect ratio of the grains, and that the specific type of order strongly affects the angle of repose and the resulting dynamics. The angle of repose is largest for dimer grains, and is smaller for both monomer and trimer grains. The motion of a collapsing grain heap consists of large scale collective shear when monomer grains are considered, but occurs via tumbling of grains on the surface for dimers, and via disorderly bulk tumbling for trimer grains.

\]

Granular assemblies have attracted much recent attention since they are a model system in which collective nonequilibrium dynamics, compaction, pattern formation, frustration, nonlinear flows, and jamming can be readily examined . The behavior they display is relevant to a diverse array of physical systems, including powders, colloids, liquids, and vortex matter in superconductors. In addition, an understanding of granular mechanics will be valuable for increasing the efficiency of a wide range of industrial processes and applications. A fundamental question that arises is the relation between the mechanical response of the granular media and the ordering of the individual grains. The nature of the close-packed ordering in 3D of different types of granular media is itself an open issue. In this work we study a simple system of granular assemblies confined in a 2D plane and show that specific types of ordering can be controllably produced by tuning the anisotropy of the grains. The resulting structure of the granular pack has a profound effect on the angle of repose and the granular dynamics.

Most previous work on the effects of grain geometry has focused on angles of repose for grains that are roughly spherical. The effect of dispersion in the sizes of individual spherical grains has been considered , as well as the effect of varying the roughness of the grain surface  or of having nonspherical grains . Anisotropic or elongated grains have attracted interest previously, but it has generally been in the context of the flow of feed grains through a hopper . Only one systematic study has considered the effect of gradual elongation of the grains on the angles of repose. In Ref. , it was observed that as the grain is elongated continuously from a monomer to a dimer, the angle of repose increases monotonically. The grain shape was continuously varied, so it would not be possible to form an ordered lattice with the same lattice constant using any two grain geometries considered.

In this work we consider with experiments and simulations a system composed of grains that can always be packed into a triangular lattice of monomers, but are individually composed of monomers, dimers, or trimers. Some previous simulation studies have focused on fully three-dimensional systems . This introduces considerable complexity to the problem since there are many more degrees of freedom, and it is difficult to determine experimentally the behavior of the grains below the surface of the pile . For this reason, we consider a system constrained to move in two dimensions. Here the microscopic behavior of individual grains can be directly observed experimentally and compared with computer simulations.

With experiments and simulations we find that the maximum angle of stability depends on the grain aspect ratio in a non-monotonic way. Both monomers, which are fully ordered, and trimers, which are both positionally and bond orientationally disordered, are less stable than dimers, which have positional order but not bond-orientational order. We directly link these specific types of ordering to the grain dynamics, which range from orderly large scale shear flow for monomers, to tumbling surface flow for dimers, to disordered bulk flow for dimers. Our results in 2D represent a first step towards the general understanding of how the microscopic structure affects the mechanical response in 3D granular assemblies.

Experiment—We consider a system of steel balls of radius 1/16” constrained to move in 2D by two Plexiglas sheets . The steel balls are either used individually, or are welded together to form dimers or linear trimers (see inset to Fig. 1). All three shapes can tile space in triangular arrays with the same lattice constant. The grains are poured into a vertical enclosed area 9” wide by 14” high, and the container is briefly shaken vertically to prepare a more ordered initial starting condition. Monomers and dimers form nearly perfect triangular lattices under these conditions, but there is always a large amount of disorder present in the trimers . The container is slowly tilted in the plane of the Plexiglas sheets to change the angle between the surface of the granular packing and gravity . We record the angle at which the first grain motion occurs, giving us the maximum angle of stability $`\mathrm{\Theta }_s`$ for the system, which is similar to the angle of repose.

In Fig. 1 we show the appearance of a system of dimers after an avalanche has occurred. The grains remained motionless until the container had been tipped by $`60^{}`$,

and then a large avalanche occurred. A significant amount of order remains in the system even after the avalanche, as seen in Fig. 1. The maximum angles of stability $`\mathrm{\Theta }_s`$ measured experimentally in this way for grains of differing anisotropy are listed in Table I. We find the lowest angle $`\mathrm{\Theta }_s=43^{}`$ for the disordered trimer system. The well-ordered monomers have a slightly higher value of $`\mathrm{\Theta }_s=45^{}`$, which is in reasonable agreement with previous experiments on spherical or nearly spherical grains . In contrast dimers display a much larger angle of stability of $`\mathrm{\Theta }_s=62^{}`$. To better understand the reasons that the dimer shape is so stable, we turn to a simulation of the system.

Simulation—We consider a system of $`N_g=705`$ grains constrained to move in a 2D plane in a system of size $`10\times 18`$ to $`60\times 18`$. We use a granular dynamics method similar to the one employed in  to integrate the equations of motion for each grain, given by: $`m_i𝐯_i=𝐟_{\mathrm{el}}^{(i)}+𝐟_{\mathrm{diss}}^{(i)}+𝐟_{\mathrm{shear}}^{(i)}+𝐟_g+𝐟_{\mathrm{fric}}`$. Here, $`𝐟_g=0.0025`$ is the force of gravity and $`𝐟_{\mathrm{fric}}=0.3`$ is the friction between the grain and the plane in which the grains move (representing the Plexiglas plate). Two grains interact only when their relative distance is smaller than the sum of their radii, $`r_g=0.4`$. In this regime, three forces are active. The elastic restoration force between two grains is $`𝐟_{\mathrm{el}}^{(i)}=_{ij}k_gm_i(|𝐫_{ij}|2r_g)𝐫_{ij}/|𝐫_{ij}|`$, where $`k_g=20`$ is the strength of the restoring spring, $`2r_g`$ is the grain diameter, $`m_i=1`$ is the grain mass and $`𝐫_{ij}=𝐫_j𝐫_i`$, the distance between grains located at $`𝐫_i`$ and $`𝐫_j`$. The dissipation force due to the inelasticity of the collision is $`𝐟_{\mathrm{diss}}^{(i)}=_{ij}\gamma m_i(𝐯_{ij}\dot{𝐫}_{ij})𝐫_{ij}/|𝐫_{ij}|^2`$, where $`\gamma =2.4`$ is a phenomenological dissipation coefficient and $`𝐯_{ij}=𝐯_i𝐯_j`$ is the relative velocity. The shear friction force mimicking solid friction is $`𝐟_{\mathrm{shear}}=\gamma _sm_i(𝐯_{ij}\dot{𝐭}_{ij})𝐭_{ij}/|𝐫_{ij}|^2`$, where $`\gamma _s=1.2`$ is the shear friction coefficient and $`𝐭_{ij}=(r_{ij}^y,r_{ij}^x)`$ is the vector $`𝐫_{ij}`$ rotated by $`90^{}`$. The behavior we observe is not sensitive to our choice of parameters.

In this model, we neglect Coulomb friction and the rotation of individual particles. In some cases, we constrain the grains to form dimers or trimers by rigidly fixing the grains together. These composite dimers and trimers are free to rotate. The walls and floor are evenly lined with immobile grains of the same size as the simulated monomers. Such a rough floor is necessary to produce a finite angle of repose in the case of monomers; without it, the monomers continue to slump until their height is reduced to a single grain. Dimers and trimers, however, produce a finite angle of repose even on a smooth floor due to their tendency to interlock.

The system is prepared in one of two ways. In studies of $`\mathrm{\Theta }_s`$, the grains are prepared in a box of size $`10\times 18`$ such that the top surface of the grains is roughly flat, either by dropping individual grains from above, or by placing them into an ordered arrangement directly. In studies of the dynamics of collapsing grains, the grains are first introduced into a box that is 1/6 as wide as the full simulation area of $`60\times 18`$, either by dropping or by direct placement. The right wall is then instantaneously moved out to the edge of the system, and the grains collapse to the right. This method has been previously employed in Refs.  to study angles of repose.

When the grains are dropped slowly from above, we find that monomers and dimers fall into well-ordered arrangements but trimers are disordered. This is illustrated in Fig. 2, where the positions of the grains are shown in Fig. 2(a,c,e), and the Delaunay triangulations of the lattice are shown in Fig. 2(b,d,f). The monomers of Fig. 2(a,b) are perfectly ordered, with no dislocations, as indicated in Table I. For dimers as in Fig. 2(c,d) the positional order is nearly complete but there is a slight chance for a vacancy to occur. Note that although the individual elements composing the dimers are ordered, the bond angles have no overall order. For trimers, shown in Fig. 2(e,f), there is some local ordering of the bond angles, but no overall order, and there is a large amount of positional disorder in the lattice.

There is a simple geometric reason for the difference in packing of the three grain shapes. If a layer of previously dropped grains is present, and there is a single monomer-sized vacancy in the layer, then an additional monomer or dimer dropped over the vacancy will fall into the vacancy and fill it. The dimer is not long enough to span the vacancy and one side of it will slip in. In contrast, the trimer is long enough to span the vacancy and block it from being filled. Additionally, a trimer is long enough to assume a continuum of angles if dropped near a second trimer lying flat, depending on the exact relative spacing of the two trimers, whereas a dimer will assume only the angles allowed by a triangular lattice. Thus with trimers only very local ordering is possible and long-range order cannot be achieved through random dropping.

We next consider $`\mathrm{\Theta }_s`$, measured with the same

procedure followed in the experiment. We simulate tipping the box slowly by rotating the direction of gravity, and record the angle $`\mathrm{\Theta }_s`$ at which the first avalanche occurs. As listed in Table I, we find the same trend as experiment: low values of $`\mathrm{\Theta }_s=27^{}`$ for ordered monomers and $`\mathrm{\Theta }_s=13^{}`$ for disordered trimers, and a much larger angle of $`\mathrm{\Theta }_s=40^{}`$ for the dimers. It is clear that the different types of order in these three geometries lead to the different angles of stability. The disordered trimers are very unstable to rearrangements since there are many dislocations in the granular packing where motion can occur. The ordered monomers are more stable due to

the absence of dislocations, but are still unstable to shearing motion along the lattice directions. In contrast, the dimers are stabilized by the absence of dislocations like the monomers, but shearing motion is suppressed due to the randomly arranged bonds between the dimers, which lead to interlocking. The larger $`\mathrm{\Theta }_s`$ we observe for dimers compared to monomers agrees with the results of Ref. . Our study of the even more anisotropic trimers indicates, however, that the increase in $`\mathrm{\Theta }_s`$ with anisotropy found in Ref.  does not continue indefinitely but instead decreases again as the grains become longer than dimers.

To further explore how the structure of the granular packing affects its stability, we directly examine the dynamics of collapsing grains when one wall is removed, and compare the three geometries in Fig. 3. We see a striking difference in the dynamics. Collapsing monomers move in large collective motions involving much of the bead pack, which follow the lattice vectors. The motion is very orderly and occurs in the form of large pulses, during which triangular wedges of grains displace along their

lattice vectors. This wedge-like motion is similar to motion seen in experiments with flowing steel balls .

In collapsing dimers, the motion is limited to the surface of the pile due to the interlocking of the dimers which prevents shearing motion along the lattice vectors. The motion that does occur is tumbling in nature and very disorderly. A relatively small portion of the bead pack is moving at any given time.

In collapsing trimers, the large number of dislocations present in the system serve as “lubrication” to the motion, allowing easy disorderly sliding of a large portion of the bead pack. There is no significant ordering present on a scale of more than a few grains so coherent motion is impossible. Motion occurs through a combination of gliding along the dislocated area, and tumbling down the topmost portion of the pile. Although the motion is widespread, it occurs more slowly than the collective motion seen in the monomers. This is due to the fact that a considerable amount of rearrangement within the bulk must occur as the movement continues, unlike the case of monomers when no rearrangement was necessary and the movement could follow the lattice vectors.

The different dynamics can also be observed by measuring the net velocity in the x-direction, $`V_x=v_x`$, of the grains over time, as shown in Fig. 4. Monomers produce extremely large pulses as the collective motion occurs, whereas dimers and trimers produce much smaller amplitude signals with their more disperse responses.

Summary— We have examined the ordering and dynamics of monomer, dimer and trimer assemblies of grains confined in 2D and relate the specific types of ordering to the dynamics of granular collapse. Dimers, which have translational order but not bond-orientational order, show the largest angle of repose. Monomers, which are completely ordered, show a smaller angle of repose that is still higher than trimers which are completely disordered. The dynamics of the monomers during collapse occurs by large scale collective shear motions. The dimers, which do not allow for easy shear motion show a surface tumbling motion. Trimers also show a tumbling motion that can occur in the bulk as well as on the surface.

Acknowledgements: We thank H. Jaeger and N. Grønbech-Jensen for helpful discussions. This work was supported in part by CLC and CULAR (LANL/UC) and by NSF, DMR-9733898.