# Contribution of polycyclic aromatic hydrocarbon ionization to neutral gas heating in galaxies: model versus observations

O. Berné, S. Foschino, F. Jalabert, C. Joblin

► To cite this version:

O. Berné, S. Foschino, F. Jalabert, C. Joblin. Contribution of polycyclic aromatic hydrocarbon ionization to neutral gas heating in galaxies: model versus observations. Astronomy and Astrophysics - A&A, 2022, 667, pp.A159. 10.1051/0004-6361/202243171. hal-03871916

HAL Id: hal-03871916

https://hal.science/hal-03871916

Submitted on 25 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Contribution of polycyclic aromatic hydrocarbon ionization to neutral gas heating in galaxies: model versus observations*

O. Berné,¹ S. Foschino¹, F. Jalabert¹, and C. Joblin¹

Institut de Recherche en Astrophysique et Planetologie, Université Toulouse III - Paul Sabatier, CNRS, CNES,
9 Av. du colonel Roche, 31028 Toulouse Cedex 04, France
e-mail: olivier.berne@irap.omp.eu

Received ?? ; accepted ??

## ABSTRACT

The ionization of polycyclic aromatic hydocarbons (PAHs), by ultraviolet (UV) photons from massive stars is expected to account for a large fraction of the heating of neutral gas in galaxies. Evaluation of this proposal, however, has been limited by our ability to directly compare observational diagnostics to the results of a molecular model describing PAH ionization. The objective of this article is to take advantage of the most recent values of molecular parameters derived from laboratory experiments and quantum chemical calculations on PAHs and provide a detailed comparison between modeled values and observational diagnostics for the PAH charge state and the heating efficiency for PAHs. Despite the use of a simple analytical model, we obtain a good agreement between model results and observational diagnostics over a wide range of radiation fields and physical conditions, in environments such as star-forming regions, galaxies, and protoplanetary disks. In addition, we found that the modeled photoelectric heating rates by PAHs are close to the observed cooling rates given by the gas emission. These results show that PAH ionization is the main source of neutral gas heating in these environments. The results of our photoelectric heating model by PAHs can thus be used to assess the contribution of UV radiative heating in galaxies (vs e.g. shock). We conclude on the importance of implementing the physics of PAH ionization in astrophysical codes, which are developed, for example, for the evaporating surfaces of protoplanetary disks, the diffuse interstellar medium, and the photodissociation regions associated with star-forming regions in the local and distant universe. We provide the empirical formulas to calculate the heating rates and heating efficiencies for PAHs in this article.

## 1. Introduction

In a seminal paper on the temperature of interstellar gas, it was proposed by Spitzer Jr (1948) that the photoelectric (PE) effect on dust grains provides an important source of heating of the neutral gas in galaxies. The central idea behind this PE heating mechanism is that the interaction with dust of ultraviolet (UV) photons from young stars, with energies up to the Lyman limit (E = 13.6 eV), releases electrons that can carry a couple eV of kinetic energy and thermalize with the gas, resulting in gas heating. Following the initial proposal of Spitzer, other authors proposed more detailed models of this mechanism (Jura 1976; de Jong 1977; Draine 1978). After the discovery that large molecules of the family of polycyclic aromatic hydrocarbons (PAHs) are ubiquitous in space (Leger & Puget 1984; Allmandola et al. 1985), it was proposed that these large molecules could contribute to the heating of the HI gas (d'Hendecourt & Léger 1987; Lepp & Dalgarno 1988). Verstraete et al. (1990) measured the photoionization yield of the PAH coronene, and estimated the contribution of a PAH population with a typical size of 80 carbons to the PE heating. They concluded that PE effect on PAHs represents a major contribution to the gas heating in cold diffuse clouds and a significant contribution in the warm phase. Bakes & Tielens (1994) and Weingartner & Draine (2001) (respectively BT94 and WD01 hereafter) developed models that take into account grain size distributions from classical dust sizes (up to typically ~0.1 µm) down to the molecular (PAH) domain. BT94 concluded that about half of the gas heating is due to grains

containing less than 1500 C atoms. While a proper description of PE heating in dust models is central to a number of astrophysical models and analysis, its implementation has only been benchmarked against observations in a limited number of studies. Deriving values of the PE heating efficiency for dust, ε, from observations requires quantifying dust emission but also the total neutral gas cooling power from all relevant infrared (IR) cooling lines in particular from atomic oxygen and ionized carbon. Several studies have used observations from Herschel, SOFIA, and Spitzer to derive the cooling budget of the neutral gas (Bernard-Salas et al. 2015; Pabst et al. 2021; Joblin et al. 2018). Recent studies of Galactic photodissociation regions (PDRs) associated with star-forming regions have compared observationally derived ε values with those calculated with the model of BT94 (Salgado et al. 2016; Salas et al. 2019; Pabst et al. 2021). There is a good qualitative agreement in terms of the shape of the curves but the calculated values are found to be larger than the observed ones by at least a factor of three. Since the earlier proposal that PAHs may provide a major contribution to the PE heating (d'Hendecourt & Léger 1987; Lepp & Dalgarno 1988; Verstraete et al. 1990), a number of observations have revealed a tight correlation between the [CII] flux and the PAH emission flux, which is in agreement with this proposal (Helou et al. 2001; Joblin et al. 2010; Lebouteiller et al. 2012). Joblin et al. (2010) mentioned that the PAH charge varies over the NGC 7023 nebula (Rapacioli et al. 2005; Berné et al. 2007) and this is expected to affect the contribution of PAHs to PE heating. Following this idea, Okada et al. (2013) provided the first study attempting to assess the connection between PAH charge and PE heating combining Spitzer and Herschel observations. The authors derived values of the ratio of gas emission to the sum of gas and PAH

<footnote>* The python code for the full model is available at https://gitlab.irap.omp.eu/COSMIC-PAH/photoelectric-heating/tree/master</footnote>
emission in several PDRs, and showed that its variation with the ionization parameter can be rationalized with a simple prescription based on previous modeling works (Tielens 2005; Weingartner & Draine 2001). This study concluded that PAHs provide a dominant contribution to the PE heating. However, it was conducted on unresolved sources which provided limited accuracy on the derived values and did not include any modeling of PAH ionization to support the conclusions. Modeling the charge state of PAHs and their contribution to the PE heating of the gas relies on the knowledge of the molecular properties of these species. Part of the relevant data have been gathered into models that describe the variation of the charge state of PAHs (e.g. Montillaud et al. 2013; Andrews et al. 2016). In addition, key molecular parameters that are needed to model the PE heating rate by PAHs of relevant sizes (i.e. > 30 carbon atoms), both in neutral and cationic forms, are the photo-absorption cross sections $\sigma_{\text{abs}}$, the ionization yields $Y$, and the ionization potentials (IPs). The earlier values of $Y$ obtained for neutral PAHs (Verstraete et al. 1990; Jochims et al. 1996) are now completed with values for cations (Wenzel et al. 2020). An additional fundamental parameter is the partition coefficient (d'Hendecourt & Léger 1987; Verstraete et al. 1990), here written $\langle\gamma_e(E)\rangle$, which determines, upon ionization by a photon of energy $E$, the fraction of the energy that is carried away in kinetic energy by the photoelectron relative to the available energy $E-IP$. In this regard, photoelectron spectral images upon UV ionization of coronene have been reported by Bréchignac et al. (2014).

In this article we present a study in which we evaluate in detail the importance of the ionization of PAHs as a source of neutral gas heating in galaxies. We first derive in Sect. 2 the observational diagnostics on the PAH charge state and the gas heating efficiency for PAHs. In Sect. 3, we present a simple analytical model of PAH ionization and gas heating by the PE effect on PAHs, which includes state-of-the-art data for molecular parameters as provided by laboratory measurements and quantum chemical calculations. In Sect. 4 we present the model results and comparisons to the models of WD01 and BT94. In Sect. 5, we compare the observational diagnostics with results obtained with our PAH model. In Sect. 6, we discuss the dominant role of PAH ionization as a source of gas heating, before concluding in Sect. 7.

## 2. Observational diagnostics of PAH charge state and efficiency of gas heating

In order to characterize the contribution of PE effect on PAHs to the gas heating, we derive two diagnostics from observations. The first diagnostic we are interested in is the ionization fraction of PAHs as derived from observations, $R_i$, which is determined by:

$$R_i = \frac{I_{\text{PAH}^+}}{I_{\text{PAH}^+} + I_{\text{PAH}^0}}, \qquad (1)$$

where $I_{\text{PAH}^+}$ and $I_{\text{PAH}^0}$ are the IR emission intensities attributed to cationic and neutral PAHs, respectively. In this study, we use the generic tool called PAHTAT (Pilleri et al. 2012) that allows one to derive $I_{\text{PAH}^+}$ and $I_{\text{PAH}^0}$ from spectroscopic IR observations<fnref n="1" />.

The second diagnostic we derive, following Okada et al. (2013), is the ratio of gas emission to the sum of gas and PAH

emission:

$$R_e = \frac{I_{\text{gas}}}{I_{\text{PAH}} + I_{\text{gas}}}, \qquad (2)$$

with $I_{\text{PAH}}$ the intensity of PAH emission and $I_{\text{gas}}$ the intensity of neutral gas cooling lines (dominated by the far-infrared [CII] and [OI] lines, see Sect. A.3). When $I_{\text{PAH}} \gg I_{\text{gas}}$, $R_e \approx \frac{I_{\text{gas}}}{I_{\text{PAH}}}$.

In addition to the two observational diagnostics described above, we derive the value of the ionization parameter:

$$\gamma = G_0 \times \sqrt{T}/n_e, \qquad (3)$$

where $T$ is the gas temperature in K, and $n_e$ is the electron density in cm⁻³. $G_0$ quantifies the intensity of the UV radiation field:

$$G_0 = \frac{1}{Z} \int_{0}^{\Omega} \int_{91.2\text{nm}}^{240\text{nm}} I_{\lambda}(\lambda) d\lambda \, d\Omega, \qquad (4)$$

where $I_{\lambda}(\lambda)$ is the intensity of the radiation field in the Solar neighborhood as given by Habing (1968) (here in W m⁻<fnref n="1" /> sr⁻<fnref n="1" /> nm⁻<fnref n="1" />), $\Omega$ is the solid angle of this radiation field and $Z = 1.68 \times 10^{-6}$ W m⁻² is the normalization factor. This definition of $G_0$ is equivalent to that of Le Petit et al. (2006). The interval of integration over wavelengths of Eq. 4 corresponds to an energy range between 13.6 and 5.17 eV. The two observational diagnostics ($R_i$ and $R_e$) and $\gamma$ are derived for the NGC 7023 nebula in spatially resolved observations, and for a sample of Galactic and extragalactic sources described in the following sections.

### 2.1. The spatially resolved NGC 7023 NW PDR

NGC 7023 is a reflection nebula illuminated by the HD 200775 double system of spectral types B3V and B5 (Finkenzeller 1985; Racine 1968) situated at 320 pc from the Sun (Benisty et al. 2013). We are interested in the PDR which is located at 42" north west (NW) of the star. NGC 7023 NW PDR is the most studied PDR and is seen almost edge-on. The dissociation front (i.e. the interface where H₂ is dissociated to 2 H atoms) is well seen (Fig. 1) through bright emission in the H2 (1-0) S(1) line at 2.12 µm, showing dense filamentary structures ($n_H \sim 10^{5-6}$cm⁻³) embedded in a more diffuse gas ($n_H \sim 10^4$cm⁻³) (Joblin et al. 2018; Lemaire et al. 1996; Fuente et al. 1996; Chokshi et al. 1988). Several estimates point to a UV radiation field intensity at the dissociation front of $G_0=2600$ (Pilleri et al. 2012; Joblin et al. 2018; Chokshi et al. 1988). We derive the values of $\gamma$ in NGC 7023 as a function of the distance to the star (see Appendix A.1) and summarize these values in Table 1. We extract the values of $R_i$ and $R_e$ (see Appendix A.2 and Appendix A.3, respectively), over the area delimited in Fig. 1.

The resulting maps are presented in Fig. 2. The derived ionization fraction for PAHs is found to decrease from regions near HD 200775 to regions beyond the dissociation front as the radiation field intensity decreases, as expected. This result is consistent with earlier studies (e.g., Pilleri et al. 2012; Boersma et al. 2016). $R_e$ follows the opposite trend, that is it increases with distance to the star. This is in line with the trend observed by Okada et al. (2013) for three position in NGC 7023. Since PAH emission is directly related to the intensity of the radiation field, this variation of $R_e$ is indicative of a better coupling between the radiation field and the gas heating in regions far from the star, where $\gamma$ is low and PAHs are mainly neutral.

<footnote n="1">¹ PAHTAT is available on the Cosmic PAH portal at https:// cosmic-pah.irap.omp.eu/</footnote>
<table><caption>Table 1: Adopted physical conditions, and derived values of the ionization parameter γ for eight angular distances from the star HD 200775 in NGC 7023 NW.</caption><thead><tr><th>d (")</th><th>20</th><th>25</th><th>30</th><th>35</th><th>40</th><th>45</th><th>50</th><th>55</th></tr></thead><tbody><tr><th>G<sub>0</sub></th><td>11004</td><td>7042</td><td>4890</td><td>3593</td><td>2751</td><td>2173</td><td>1760</td><td>733</td></tr><tr><th>T (K)</th><td>750</td><td>750</td><td>750</td><td>750</td><td>750</td><td>750</td><td>350</td><td>150</td></tr><tr><th>n<sub>H</sub> (10<sup>4</sup>cm<sup>-3</sup>)</th><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><th>γ (K<sup>1/2</sup>cm<sup>3</sup>)</th><td>125567</td><td>80363</td><td>55807</td><td>41001</td><td>23543</td><td>18602</td><td>10293</td><td>2807</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Fig. 1: General view of the NGC 7023 NW PDR. Grey image is Spitzer-IRAC data at 8 µm. Green contours show the emission in the H₂ 1-0 S(1) line at 0.5, 1 and 2×10⁻⁷ W m⁻<fnref n="2" /> sr⁻¹, from observations at CFHT presented in Lemaire et al. (1996). The orange contour defines the region considered here for the PE heating study.</figcaption></figure>

## 2.2. Sample of Galactic and extragalactic sources

Values for $R_e$ are derived from observations for a number of astrophysical environments for which we expect the PE effect on PAHs to be a major contributor to gas heating. These include star-forming regions, protoplanetary disks, star-forming galaxies, ultra luminous infrared galaxies, and a high-redshift galaxy. The methodology to derive these parameters is detailed in Appendix B. All derived values are summarized in Table 2 together with associated $\gamma$ values, and are shown in Fig. 3. Similar to the trend observed in NGC 7023 NW, $R_e$ values are found to be lower in sources with high $\gamma$ values, such as protoplanetary disks or galaxies that actively form stars. However, the trend is also highly variable. For example, the Orion Bar and Ced 201 have similar $\gamma$ values but their values of $R_e$ differ by an order of magnitude.

## 3. Molecular model of the charge state and photoelectric emission from PAHs

Here we present a simple molecular model, whose objective is to compute the PAH ionization fraction, the heating efficiency for PAHs, and the heating rates by PAHs. This model is based on the photophysics of isolated interstellar PAHs and includes the latest results from laboratory experiments for PAH molecular parameters. Upon absorption of a UV photon, three main processes are expected to be in competition to relax the energy: ionization, dissociation, and radiative cooling (see, e.g., Joblin et al. 2020).

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2: Maps of the ionization fraction of PAHs in NGC 7023 (upper panel), and ratio of gas cooling to PAH emission (lower panel).</figcaption></figure>

Photoelectric heating is related to ionization, and therefore we focus on this process in the following.

### 3.1. Selection of species and molecular parameters

A first step consists in selecting species for which molecular data is available. Here we rely on information given by the *Theoretical spectral database of PAHs*<fnref n="2" /> (Malloci et al. 2007). The chosen molecules are listed in Table 3, having compact shapes and different characteristic sizes relevant for interstellar PAHs (Leger & Puget 1984; Leger et al. 1989; Allamandola et al. 1985; Montillaud et al. 2013; Andrews et al. 2016).

In the present model, we consider four charge states that is $Z \in \{-1, 0, 1, 2\}$). Although there is no clear spectroscopic evidence for the presence of anions ($Z \le -1$) in the ISM, we include

<footnote n="2">² http://astrochemistry.oa-cagliari.inaf.it/database/pahs.html</footnote>
<table><caption>Table 2: Summary of values of $R_e$ and ionization parameter $\gamma$ derived from observations.</caption><thead><tr><th>Object</th><th><i>R</i><sub>e</sub></th><th><i>γ</i></th></tr></thead><tbody><tr><td colspan="3">Galactic photodissociation regions</td></tr><tr><td>Orion Veil</td><td>0.054 (0.024,0.091)</td><td>997 (595,1850)</td></tr><tr><td>Orion Bar</td><td>[0.07-0.11]</td><td>[0.7 × 10<sup>4</sup> - 2 × 10<sup>4</sup>]</td></tr><tr><td>IC 63</td><td>0.14 (0.11,0.17)</td><td>[500-1500]</td></tr><tr><td>Ced 201</td><td>0.018 (0.012,0.024)</td><td>[6000 - 10 000]</td></tr><tr><td>Horsehead</td><td>0.11 (0.08,0.14)</td><td>[30-1400]</td></tr><tr><td colspan="3">Other environments</td></tr><tr><td>Disks around HAeBe stars</td><td>0.0037 (9.0 × 10<sup>-4</sup>, 0.0053)</td><td>[3 × 10<sup>4</sup> - 3 × 10<sup>6</sup>]</td></tr><tr><td>Diffuse ISM</td><td>0.021 ± 0.002</td><td>[1580 - 1900]</td></tr><tr><td>LMC</td><td>[0.05-0.1]</td><td>[5460 - 10 900]</td></tr><tr><td>Star-forming Galaxies</td><td>[0.0048 - 0.0073]</td><td>[1.4 × 10<sup>3</sup> - 4 × 10<sup>3</sup>]</td></tr><tr><td>ULIRGs</td><td>0.034 (0.011, 0.12)</td><td>6.3 × 10<sup>4</sup>(3.4 × 10<sup>4</sup>, 1.9 × 10<sup>5</sup>)</td></tr><tr><td>High z Galaxy</td><td>[0.008-0.022]</td><td>[5 × 10<sup>4</sup> - 2 × 10<sup>5</sup>]</td></tr></tbody></table>

**Notes.** The listed values are either median with confidence intervals (see Appendix B details on the boundaries of these intervals), or provide the range of considered values in brackets.

<figure><img src="image_5.png" /><figcaption>Fig. 3: Measured emission ratio $R_e$ (Eq. 2), as a function of physical conditions traced by the ionization parameter $\gamma$, in a sample of astrophysical objects.</figcaption></figure>

<table><caption>Table 3: PAH species considered in the present study.</caption><thead><tr><th>Name</th><th>Chemical formula</th></tr></thead><tbody><tr><td>Ovalene</td><td>C<sub>32</sub>H<sub>14</sub></td></tr><tr><td>Tetrabenzocoronene</td><td>C<sub>36</sub>H<sub>16</sub></td></tr><tr><td>Circumbiphenyl</td><td>C<sub>38</sub>H<sub>16</sub></td></tr><tr><td>Circumanthracene</td><td>C<sub>40</sub>H<sub>16</sub></td></tr><tr><td>Circumpyrene</td><td>C<sub>42</sub>H<sub>16</sub></td></tr><tr><td>Hexabenzocoronene</td><td>C<sub>42</sub>H<sub>18</sub></td></tr><tr><td>Dicoronylene</td><td>C<sub>48</sub>H<sub>20</sub></td></tr><tr><td>Circumcoronene</td><td>C<sub>54</sub>H<sub>18</sub></td></tr><tr><td>Circumovalene</td><td>C<sub>66</sub>H<sub>20</sub></td></tr></tbody></table>

them because they could play a role in gas heating in environments with low $\gamma$ values. Cations ($Z = 1$) and neutrals ($Z = 0$) have long been proposed to exist in PDRs (see review in Tielens 2008). Theoretical (Malloci et al. 2007) and experimental (Wenzel et al. 2020) data have shown that the dication stage ($Z = 2$) is accessible in PDRs for the considered species (Table 3), because their ionization potential IP² is lower than the Lyman limit of

13.6 eV. Data for IP<fnref n="3" /> values for PAHs are scarce, but Zhen et al. (2016) have shown that they are greater than 13.6 eV for small PAHs, and Wenzel et al. (2020) demonstrate that this is probably the case up to a typical size of $N_C \sim 50$. Therefore, we do not consider trications in our model.

Values of the ionization potentials of a PAH with a given size $N_C$, are obtained using the empirical formalism of WD01, with updated parameters from Wenzel et al. (2020) as given in Appendix C.1. Several studies provide measurements of the ionization yields of neutral (Jochims et al. 1996; Verstraete et al. 1990) and cationic (Zhen et al. 2016; Wenzel et al. 2020) PAHs. From these measurements, Jochims et al. (1996) and Wenzel et al. (2020) derived empirical laws for the ionization yields of neutral and cationic PAHs, respectively, which are used in this work (see Appendix C.2). The ionization yield curves are shown as a function of incident photon energy in Fig. 4 for three PAH sizes. Similarly to Visser et al. (2007), we use $Y^-(E) = 1$. The threshold in energy for the detachment of the electron on anions is given by the electron affinity (e.g., Tschurl & Boesl 2006). The latter is very variable with the molecule but remains below 2 eV (see *Theoretical spectral database of PAHs*; Malloci et al. 2005; Wahab et al. 2022). We therefore assumed, as an approximation, that $Y^-(E) = 1$ at all energies. We present another case in Appendix C.3.

Finally, we use the photoabsorption cross sections given by the *Theoretical spectral database of PAHs*<fnref n="3" /> (see also Malloci et al. 2004). First, we retrieve the photoabsorption cross section for each molecule in Table 3 and for each value of $Z \in \{-1, 0, 1, 2\}$ and normalize it to the number of carbon atoms $N_C$. Second, for each value of $Z \in \{-1, 0, 1, 2\}$ we compute the averaged normalized photoabsorption cross section. This yields an average photoabsorption cross section per carbon atom $\sigma(E, Z)$, from which the ionization cross section per carbon atom can be derived using:

$$ \sigma_{\text{ion}}(E, Z) = Y^Z(E) \times \sigma(E, Z) \ (\text{cm}^2 \text{ C}^{-1}). \qquad (5) $$

Figure 5 represents the $\sigma(E, Z)$ and $\sigma_{\text{ion}}(E, Z)$ curves as used in the model.

<footnote n="3">³ https://www.dsf.unica.it/~gmalloci/pahs/</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4: Ionization yields $Y^Z(E)$ for charge states $Z = 0$ (solid lines) and $Z = 1$ (dashed lines), computed from Eqs. C.3 and C.4-C.5 respectively, and for the three different PAH sizes. The yield for $Z = -1$, $Y^-(E) = 1$ is not shown, for clarity of the figure.</figcaption></figure>

## 3.2. PAH charge states and ionization fraction

At equilibrium the fractions $f^Z$ of each PAH population in a Z charge state is given by:

$$ \left\{ \begin{array}{ll} \text{charge level } Z = -1: & f^{-}k_{\text{det}} = f^{0}k_{\text{att}}n_{e} \\ \text{charge level } Z = 0: & f^{0}k_{\text{pe}}^{0} + f^{0}k_{\text{att}}n_{e} = f^{+}k_{\text{rec}}^{+}n_{e} + f^{-}k_{\text{det}} \\ \text{charge level } Z = 1: & f^{+}k_{\text{rec}}^{+}n_{e} + f^{+}k_{\text{pe}}^{+} = f^{2}k_{\text{rec}}^{2}n_{e} + f^{0}k_{\text{pe}}^{+} \\ \text{charge level } Z = 2: & f^{2}k_{\text{rec}}^{2}n_{e} = f^{+}k_{\text{pe}}^{+} \end{array} \right. \qquad (6) $$

where $k_{\text{pe}}^Z$ and $k_{\text{rec}}^Z$ are respectively the ionization and recombination rates of PAHs in the $Z \in \{0, 1, 2\}$ charge state. We note that when $Z$ is an exponent in the notation, it is written $\{-, 0, +, 2+\}$ rather than $\{-1, 0, 1, 2\}$, to avoid confusion. We also introduce the attachment rate of an electron to a neutral PAH, $k_{\text{att}} \equiv k_{\text{rec}}^0$, and the detachment rate of the electron from an anion, $k_{\text{det}} \equiv k_{\text{pe}}^-$.

$$ \sum_{Z_{\text{min}}}^{Z_{\text{max}}} f^Z = 1. \qquad (7) $$

For $Z \in \{-1, 0, 1, 2\}$, Eq. 6 and Eq. 7 yield

$$ f^{-} = \left( 1 + \frac{k_{\text{det}}}{k_{\text{att}}n_{e}} + \frac{k_{\text{det}}k_{\text{pe}}^{0}}{k_{\text{att}}k_{\text{rec}}^{+}n_{e}^{2}} + \frac{k_{\text{det}}k_{\text{pe}}^{0}k_{\text{pe}}^{+}}{k_{\text{att}}k_{\text{rec}}^{+}k_{\text{rec}}^{2}n_{e}^{3}} \right)^{-1}, \qquad (8) $$

$$ f^{0} = \left( 1 + \frac{k_{\text{att}}n_{e}}{k_{\text{det}}} + \frac{k_{\text{pe}}^{0}}{k_{\text{rec}}^{+}n_{e}} + \frac{k_{\text{pe}}^{0}k_{\text{pe}}^{+}}{k_{\text{rec}}^{+}k_{\text{rec}}^{2}n_{e}^{2}} \right)^{-1}, \qquad (9) $$

$$ f^{+} = \left( 1 + \frac{k_{\text{rec}}^{+}n_{e}}{k_{\text{pe}}^{0}} + \frac{k_{\text{pe}}^{+}}{k_{\text{rec}}^{2} + \frac{k_{\text{att}}k_{\text{rec}}^{+}n_{e}^{2}}{k_{\text{det}}k_{\text{pe}}^{0}}} \right)^{-1}, \qquad (10) $$

$$ f^{2+} = \left( 1 + \frac{k_{\text{rec}}^{2+}n_{e}}{k_{\text{pe}}^{+}} + \frac{k_{\text{rec}}^{+}k_{\text{rec}}^{2+}n_{e}^{2}}{k_{\text{pe}}^{0}k_{\text{pe}}^{+}} + \frac{k_{\text{att}}k_{\text{rec}}^{+}k_{\text{rec}}^{2+}n_{e}^{3}}{k_{\text{det}}k_{\text{pe}}^{+}k_{\text{pe}}^{+}} \right)^{-1}. \qquad (11) $$

<figure><img src="image_7.png" /><figcaption>Fig. 5: Average photoabsorption ($\sigma$, solid line) and photoionization ($\sigma_{\text{ion}}$, dashed line) cross sections per C atom of PAHs adopted in the 4 energy level model.</figcaption></figure>

The ionization rates are given by:

$$ k_{\text{pe}}^{Z} = \int_{0}^{\Omega} \int_{IP(Z+1)}^{13.6} \frac{\sigma_{\text{ion}}(E, Z)I_{\text{E}}(E)}{E} dE \, d\Omega \quad (\text{s}^{-1}), \qquad (12) $$

where $I_E(E)$, is the local intensity of the UV radiation field in W m$^2$ eV$^{-1}$ sr$^{-1}$.

The values of $\sigma_{\text{ion}}(E, Z)$ have been defined above (Eq. 5, in particular). On the other hand, the recombination rate is calculated using Spitzer's formalism (Spitzer 2004), adapted for PAH cations by Verstraete et al. (1990):

$$ k_{\text{rec}} = 1.28 \times 10^{-10} N_{\text{C}} \sqrt{T}(1 + \Phi) \, (\text{cm}^3 \, \text{s}^{-1}), \qquad (13) $$

where T is the gas temperature, and

$$ \Phi = \frac{eU}{k_{\text{B}}T} = \frac{1.85 \times 10^5}{T \sqrt{N_{\text{C}}} \sqrt{T}} \quad (\text{dimensionless}). \qquad (14) $$
where e is the electron charge, U is the mean electrostatic potential evaluated at the radius of the PAH considered, defined by $U = \frac{1}{4\pi\epsilon_0} \frac{e}{0.9 \times 10^{-10} \sqrt{N_C}}$ (J C⁻¹), k_B the Boltzmann constant and T the gas temperature.

Eq. 13 can be extended to include all charge states Z > 0 in the recombination rate:

$$k_{\text{rec}}^{Z} = 1.28 \times 10^{-10} N_{\text{C}} \sqrt{T} (1 + \Phi \times (1 + Z)) \ (\text{cm}^{3} \text{ s}^{-1}). \quad (15)$$

For k\_att, we use the work of Carelli et al. (2013) who provide an empirical fit to their quantum chemistry calculations that is

$$k_{\text{att}} = a \left( \frac{T}{300\text{K}} \right)^b \exp\left(-\frac{c}{T}\right) \ (\text{cm}^3 \text{ s}^{-1}). \quad (16)$$

For the parameters a, b, c we use the values determined by Carelli et al. (2013) for coronene, C₂₄H₁₂ which is the largest molecule in this study, that is a = 2.74 × 10⁻⁹ cm³s⁻¹, b = 0.11, and c = -1.11. We note however that for all the molecules studied by Carelli et al. (2013), the attachment rate is of the order of ~10⁻⁹ cm³s⁻¹ with only a small dependence on temperature. Hence this choice is unlikely to affect the model results. Finally, we write the ensemble of charge fractions:

$$\mathbf{F} = \{f^{++}, f^{+}, f^{0}, f^{-}\}, \quad (17)$$

and the modeled PAH positive ionization fraction is defined as:

$$f_i = f^+ + f^{2+}. \quad (18)$$

### 3.3. Photoelectric heating efficiency for PAHs

The power injected into the gas via the photoelectrons, P_e, is given by:

$$P_e = f^- \times P_e^- + f^0 \times P_e^0 + f^+ \times P_e^+ \ (\text{W}), \quad (19)$$

where f⁻, f⁰ and f⁺ are the fractions of anionic, neutral and cationic PAHs, derived from Eqs 8, 9 and 10, respectively. P_e⁻, P_e⁰ and P_e⁺ are the powers injected in the gas by photoelectrons ejected from anionic, neutral and cationic PAHs, respectively, which are defined by:

$$P_e^Z = \int_0^{\Omega} \int_{IP^{Z+1}}^{13.6} \gamma_e(E) (E - IP^{Z+1}) \times \sigma_{\text{ion}}(E, Z) \frac{I_E(E)}{E} dE \ d\Omega \ (\text{W}), \quad (20)$$

where γ_e(E) is the fraction of the energy that, after ionization, goes into kinetic energy of the photoelectron.

Energy conservation upon ionization requires:

$$E - IP = E^*(E) + E_K(E) + E_K^M(E), \quad (21)$$

where E* is the internal energy of the PAH and E_K^M is its kinetic energy, which, due to the much higher mass of the PAH compared to the electron and due to conservation of momentum, is negligible compared to the kinetic energy of the electron, E_K. Next, following d'Hendecourt & Léger (1987), we define the partition coefficient 〈γ_e(E)〉 as the energy averaged γ_e(E):

$$\langle \gamma_e(E) \rangle = \left\langle \frac{E_K(E)}{E - IP} \right\rangle. \quad (22)$$

We estimate the value of 〈γ_e(E)〉 from photoelectron spectroscopy measurements performed using VUV synchrotron radiation and for the coronene molecule (Bréchignac et al. 2014).

The derived value of 〈γ_e(E)〉 = 0.46 ± 0.06 (68% confidence interval) is found to be in agreement with the value of 〈γ_e(E)〉 = 0.5, which was used in a number of articles (e.g., Verstraete et al. 1990; Tielens 2005; Bakes & Tielens 1994) and derived from older measurements on benzene by Terenin & Vilessov (1964) (see Appendix C.4). For anionic PAHs, we assume that all energy available is taken by the electron. For IP⁰ = 0, this means that all the energy of the photon is transformed into kinetic energy of the electron.

On the other hand, the power absorbed, mostly (but not only) in the UV, by each PAH charge state, P_Rad, is given by:

$$P_{\text{Rad}}^{Z} = \int_{0}^{\Omega} \int_{0}^{13.6 \text{ eV}} \sigma(E, Z) I_{E}(E) \, dE \, d\Omega \ (\text{W}). \quad (23)$$

The total radiative power absorbed by PAHs, P_Rad, is then defined by:

$$P_{\text{Rad}} = f^{-} \times P_{\text{Rad}}^{-} + f^{0} \times P_{\text{Rad}}^{0} + f^{+} \times P_{\text{Rad}}^{+} + f^{2+} \times P_{\text{Rad}}^{2+} \ (\text{W}). \quad (24)$$

We define the heating efficiency for PAHs, ε\_PAH, as the ratio between the power injected in the gas via photoelectrons from PAHs to the total power of the radiation absorbed by PAHs, that is

$$\epsilon_{\text{PAH}} = \frac{P_e}{P_{\text{Rad}}}. \quad (25)$$

We note that the derived values for the radiative power and heating efficiency depend on N_C, through absorption and ionization cross sections, and through recombination rates with electrons. For the sake of simplicity this variable has, however, not been included in the notation. Also, our definition of the PE heating efficiency for PAHs is similar to that used by WD01 to compute their heating efficiency for grains, ε_T, but differs from the definition of ε by BT94, who considered only absorbed photons in the 6-13.6 eV range. Finally, the PE heating rate by PAHs can be derived using the following relation:

$$\Gamma_{\text{PAH}} = P_e \times \frac{f_C[C]}{N_C} = \epsilon_{\text{PAH}} P_{\text{Rad}} \frac{f_C[C]}{N_C} \ (\text{W H}^{-1}), \quad (26)$$

where f_C is the fraction of cosmic carbon locked in PAHs, and [C] is the cosmic abundance of carbon relative to hydrogen. The python code to compute Γ\_PAH and ε\_PAH with the model described in this section is provided through the Cosmic PAH portal at https://cosmic-pah.irap.omp.eu/ or directly at https://gitlab.irap.omp.eu/COSMIC-PAH/ photoelectric-heating/tree/master.

## 4. Model results and comparison with previous models

### 4.1. Charge balance

Figure 6 shows the variation of the population fractions of PAHs in charge state Z, as a function of γ, derived from the model, for T = 1000 K, N_C = 54, and stellar effective temperature of T\_eff = 3 × 10⁴ K using a Kurucz (1993) star model. For γ ≤ 10³ (e.g. diffuse ISM or dark clouds), the population is dominated by neutral PAHs, with an anionic contribution up to 20%. A comparable fraction of neutrals and anions was also derived at low γ values by both BT94 and WD01 for 5 Å grains (comparable to PAHs of size N_C = 54). Neutral PAHs dominate up to γ ~ 10⁴, followed by cations with a maximum at about γ ~ 3 × 10⁴. Finally, dications are the dominant species for γ ≥ 10⁵. For
<figure><img src="image_8.png" /><figcaption>Fig. 6: Fractions in the different PAH charge state Z = -1, 0, 1, 2 as a function of the ionization parameter γ, which were derived using our PAH model and $T_{eff} = 3 \times 10^4$ K.</figcaption></figure>

$\gamma = 10^5$ and the same stellar effective temperature, WD01 find values for the charge fractions (Eq. 17) of $\mathbf{F} = \{0.7, 0.2, 0.1, 0.0\}$. We infer comparable fractions but with less dications and more cations, $\mathbf{F} = \{0.5, 0.4, 0.1, 0.0\}$ (Fig. 6). The difference results from the higher ratio of ionization rate to recombination rate in WD01 compared to our model for high $\gamma$ values. On the contrary, with the same stellar effective temperature and the same gas temperature, BT94 find at $\gamma = 4 \times 10^4$ $\mathbf{F} = \{0, 0.4, 0.6, 0.0\}$, and $\mathbf{F} = \{0.2, 0.5, 0.3, 0.0\}$, for planar and spherical particles, respectively, resulting in lower ionized fractions compared to our model with $\mathbf{F} = \{0.4, 0.5, 0.1, 0\}$ (Fig. 6).

## 4.2. Heating efficiency

<figure><img src="image_9.png" /><figcaption>Fig. 7: Photoelectric efficiency $\epsilon_{PAH}$ computed for $N_C = 54$, $T_{eff} = 30\,000$ K, and $T=100$ K and $1000$ K with the PAH heating model. For comparison, we give the photoelectric efficiency $\epsilon_T$ computed by WD01 for grains with a size $a = 4.9\AA$ equivalent to $N_C = 54$ (Eq. C.2), and the heating efficiency for PAHs given in Tielens (2021).</figcaption></figure>

In Fig. 7, we present the results of the computation of $\epsilon_{PAH}$ with the PE heating model by PAHs for $T_{eff} = 30\,000$ K, $N_C = 54$, and two values of the gas temperature, $T = 100$ K and $T = 1000$ K. We compare $\epsilon_{PAH}$ with the efficiencies $\epsilon_T$ derived

for the same physical conditions by WD01, for grains of size $a = 4.9\AA$, corresponding to $N_C = 54$ (Eq. C.2). For $\gamma = 100$ we find a good agreement between the two models, for both gas temperatures. For higher values of $\gamma$, $\epsilon_{PAH}$ is systematically higher than $\epsilon_T$. The difference increases with increasing $\gamma$, up to about an order of magnitude for values of $\gamma \gtrsim 10^5$. This can be rationalized by a slower decrease of the ionization to recombination rate ratio with $\gamma$ in our model as compared to WD01 (see previous section). Tielens (2021) provides the following formula for the PAH heating efficiency:

$$ \epsilon_{PAH} = \frac{0.06}{1 + 7 \times 10^{-5} \gamma}, \qquad (27) $$

which we also display in Fig. 7. This curve is also below the heating efficiency we infer for PAHs. At $\gamma \lesssim 10^3$ this is due to the absence of photo-detachment in Eq. 27. For larger values of $\gamma \gtrsim 10^3$, the difference stems, notably, from the use of a photoelectron yield $Y(E) = 0.3$ at 10 eV by Tielens (2021), while we use a value closer to $Y(E) = 0.5$ (see Fig. 4). Since BT94 do not provide the photoelectric heating efficiency curves for different sizes, we cannot directly compare our results to their model.

## 4.3. Heating rates

<figure><img src="image_10.png" /><figcaption>Fig. 8: Photoelectric heating rate of the gas. Dash-dotted line: average model from Weingartner & Draine (2001) (see text for details). Dashed line: model from Bakes & Tielens (1994). Continuous lines: PAH model (this work) for two PAH abundances, that is $f_C = 0.05, 0.1$.</figcaption></figure>

The PAH heating rate we derived as a function of $\gamma$ is shown in Fig. 8 for $T = 100$ K, $N_C = 54$, $T_{eff} = 3 \times 10^4$ K, and two PAH abundances. We compare these values to the total heating rate of WD01 averaged over several dust size distributions (Fig. 15 of WD01) and that of BT94, obtained for the same physical conditions. Overall, the curves agree well. This implies that, with a purely molecular model, PAHs can produce a heating rate comparable to that of earlier models including a full size distribution of dust grains. We can also compare the results of our model with those of Verstraete et al. (1990) for the diffuse ISM. The physical conditions are $G_0 = 1.7$ with the Draine (1978) field, a gas temperature of $T = 80$ K, and an electron density $n_e = 3 \times 10^{-2}$ cm$^{-3}$ yielding $\gamma \sim 500$ K$^{1/2}$ cm$^{-3}$. With a value of $f_C = 0.1$, as in Verstraete et al. (1990), our model yields $\Gamma_{PAH} = 2.7 \times 10^{-33}$ W H$^{-1}$
(Tab. 4). The latter authors find a value that is larger by a factor ~3, that is $Γ_{PAH} = 8 × 10^{-33}$ W H$^{-1}$. This is probably due to a combination of the following factors: 1) the higher ionization yield used by these authors, for example $Y \sim 0.7$ at 10 eV compared to $Y \sim 0.5$ in our model (Fig. 4), 2) their larger value for the partition coefficient $\langle γ_e(E) \rangle = 0.5$ compared to 0.46 in our model, 3) the fact that we include the ionization of cations in our model, which affects the charge balance and reduces the heating rate, because the cations do not heat the gas.

<table><caption>Table 4: Gas heating efficiency ($ε_{PAH}$) and heating rate ($Γ_{PAH}$) for the diffuse ISM, computed with the PAH model.</caption><thead><tr><th>ε<sub>PAH</sub></th><th>Γ<sub>PAH</sub> (10<sup>-33</sup> W H<sup>-1</sup>)</th><th>Radiation field</th></tr><tr><td colspan="3">This work</td></tr></thead><tbody><tr><td>0.037</td><td>2.7</td><td>Draine (1978)</td></tr><tr><td>0.059</td><td>1.6</td><td>Habing (1968)</td></tr><tr><td>0.039</td><td>2.0</td><td>Mathis et al. (1983)</td></tr></tbody></table>

**Notes.** Computed with the PAH model, for $T = 80$ K, $f_c = 0.1$ and $n_e = 3 × 10^{-2}$, using three different model spectra of UV interstellar radiation field.

### 4.4. Empirical formulas for the heating efficiency and rate

We computed $ε_{PAH}$ using the molecular model presented in Sect. 3, for a fixed size $N_C = 54$ (corresponding to circum-coronene, see Table 3), and a gas temperature $T = 500$ K, using radiation fields with stellar spectral types at various effective temperatures ($T_{eff}$ in the $10^4-4×10^4$ K range). For $T_{eff} ≥ 3×10^4$, we used stellar spectra from the Pollux database 4 (Palacios et al. 2010) computed with the CMFGEN code (Hillier & Miller 1998). The lower temperature stellar spectra are from Kurucz (1993) database. The resulting PE efficiency curves for PAHs as a function of the ionization parameter $γ = G_0 \sqrt{T}/n_e$ are shown in Fig. 9. The values of $ε_{PAH}$ can be fitted using the following empirical formula:

$$ε_{PAH} \approx \frac{a}{1+b \times γ} + \frac{c}{1+d \times γ}. \qquad (28)$$

This is similar to the expression, proposed by Tielens (2005) (see Eq. 27) but includes an additional term to take into account photo-detachment at low $γ$ values. We fit the analytical expression of Eq. 28 to the model derived values $ε_{PAH}$ (Fig. 9) for all radiation fields, and provide the obtained values of the parameters in Table 5. Using Eq. 24 we computed $P_{Rad}$ for $N_C = 54$ for $T = 500$K and $T_{eff}$ between $10^4$ and $4 × 10^4$K, as a function of $G_0$. The result of this calculation, combined to Eq. 26 using $N_C = 54$, and $[C] = 2.7 × 10^{-4}$ per H atom (Tielens 2021) yields an average heating rate for PAHs:

$$Γ_{PAH} = 3.0 \pm 0.6 \times 10^{-31} ε_{PAH} G_0 f_C \ (W H^{-1}), \qquad (29)$$

where the uncertainty on the numerical coefficient is the standard deviation among all numerical values obtained with the model.

## 5. Charge state and photo-electric heating of PAHs: Model vs observations

In this section we discuss the comparison between the observational results presented in Sect. 2 and the results of the molecular model described in Sect. 3. More specifically, for the case

<table><caption>Table 5: Parameters of the fit with analytical expression in Eq. 28 to the model derived curves of the photoelectric heating efficiency for PAHs, $ε_{PAH}$.</caption><thead><tr><th>T<sub>eff</sub>(/10<sup>4</sup>)</th><th>S. Type</th><th>a</th><th>b</th><th>c</th><th>d(/10<sup>-5</sup>)</th></tr></thead><tbody><tr><td>4.00</td><td>O8</td><td>0.949</td><td>0.641</td><td>0.130</td><td>7.45</td></tr><tr><td>3.00</td><td>O9.5</td><td>0.947</td><td>0.646</td><td>0.119</td><td>6.77</td></tr><tr><td>2.50</td><td>B0.5</td><td>0.948</td><td>0.645</td><td>0.108</td><td>6.20</td></tr><tr><td>1.50</td><td>B5</td><td>0.947</td><td>0.660</td><td>0.089</td><td>5.20</td></tr><tr><td>2.00</td><td>B2.5</td><td>0.947</td><td>0.732</td><td>0.0516</td><td>3.54</td></tr><tr><td>1.50</td><td>B5</td><td>0.968</td><td>0.880</td><td>0.0267</td><td>2.60</td></tr><tr><td>1.25</td><td>B9</td><td>0.985</td><td>1.08</td><td>0.015</td><td>2.18</td></tr><tr><td>1.00</td><td>A0</td><td>0.992</td><td>1.35</td><td>0.00889</td><td>1.91</td></tr></tbody></table>

**Notes.** Fit results are provided for radiation fields of varying effective temperatures $T_{eff}$ and corresponding spectral types. All fits are computed on 1001 model points and yield a determination coefficient $R^2 \approx 0.99$.

<figure><img src="image_11.png" /><figcaption>Fig. 9: Observed emission ratio $R_e$ (Table 2) and modeled photoelectric heating efficiency for PAHs, $ε_{PAH}$, as a function of the ionization parameter $γ$, for various stellar effective temperatures. The upper curve corresponds to $T_{eff} = 4 × 10^4$K, and the lowest curve to $T_{eff} = 10^4$K. All models use $N_C = 54$ and a gas temperature of $T = 500$ K.</figcaption></figure>

of NGC 7023, we compare $R_i$, the ionization fraction of PAHs derived from the analysis of the mid-infrared spectra of the PDR (see Appendix A.2) to $f_i$, the modeled PAH ionization fraction. We also compare the emission ratio, $R_e$, to computed values of $ε_{PAH}$. If PAHs are the main source of gas heating, then the later two values and their variations should be in close agreement. We then extend this comparison to our sample of galactic and extra-galactic sources. In addition, for NGC 7023, the Orion Bar, and the diffuse ISM, we compare the total gas cooling $Λ_{gas}$ to the modeled heating rate by PAHs $Γ_{PAH}$.

### 5.1. Comparison of model results with observations in NGC 7023

We first computed the model for three PAH sizes that is $N_C = 96, 54, 32$, using the physical parameters in Tab. 1. and the stellar spectrum of HD 200775 (see Appendix A.1). The gas temper-

<footnote n="4">4 http://npollux.lupm.univ-montp2.fr</footnote>
ature in NGC 7023 NW varies between 750 and 150 K. We there-
fore adopt a fixed value of $T = 500$ K, since this parameter has
only a small effect on the model results (Fig. 7). The comparison
of model results with the observational diagnostics is provided
for the eight positions listed in Table 1. In Fig. 10 we compare
the ionization fraction obtained with the model, $f_i$, to those de-
rived from observations in NGC 7023, $R_i$, as a function of $\gamma$. The
obtained good agreement suggests that the model is able to quan-
tify the evolution of the charge state of PAHs in NGC 7023. The
PAH size plays a small role in the values of $f_i$, but better results
are obtained for PAHs with $N_C \sim 54$. In Fig. 11, we compare $R_e$
as measured in NGC 7023 with values of the heating efficiency,
$\epsilon_{PAH}$, derived from the model for three PAH sizes. The agree-
ment in terms of curve shape and absolute values is good, with
all observed values intersecting the model, except for the closest
positions to the star. There are a number of parameters that could
explain this discrepancy including less precise values for the
physical conditions in this region, or a possible additional source
of heating in this environment, such as shocks due to the expan-
sion of a small HII region. This HII region could also contribute
to some fraction of the [CII] emission that we include in the total
cooling. Overall, however, $R_e$ compares well to $\epsilon_{PAH}$. Finally, we
can also compare the total photoelectric heating rate by PAHs
$\Gamma_{PAH}$ to the total gas cooling rate $\Lambda_{gas} = 4\pi \times I_{gas}/N_H$, with $N_H$
the total hydrogen column density. Joblin et al. (2010) provide
values for $N_H$ in the NGC 7023 PDR, which are $N_H = 10.5, 8.1$
and $3.9 \times 10^{21} \text{ cm}^{-2}$, at 50, 48 and 47" from the star, respectively.
We thus adopt $N_H = 7.5 \pm 3.3 \times 10^{21} \text{ cm}^{-2}$. At distances 45 - 50"
from the star, we derive $I_{gas} = 3.8 \pm 1.6 \times 10^{-6} \text{ W m}^{-2} \text{ sr}^{-1}$. This
yields a cooling rate $\Lambda_{gas} = 6.4 \pm 4.6 \times 10^{-31} \text{ W H}^{-1}$ for this range
of distances. The heating rates by PAHs derived from the model
for the physical conditions at 45 and 50" (Table 1) and a value
$f_C = 0.05$ are, for both positions, $\Gamma_{PAH} = 6.5 \times 10^{-31} \text{ W H}^{-1}$, in
agreement with the cooling rates derived from the observations.

<figure><img src="image_12.png" /></figure>

Fig. 10: Variation with $\gamma$ of the PAH ionization fraction $R_i$ in NGC 7023 (diamonds with errobars). Values derived from the map in Fig. 2 for the eight angular distances between 20 and 55" from HD 200775 given in Table 1. The values of $f_i$ from the PAH model are shown with lines for a gas temperature $T = 500$K, and three PAH sizes.

<figure><img src="image_13.png" /></figure>

Fig. 11: Variation with $\gamma$ of the emission ratio $R_e$ in NGC 7023 derived from the map in Fig. 2 for the eight angular distances between 20 and 55" from HD 200775 given in Table 1. The values of $\epsilon_{PAH}$ from the PAH model are shown with lines for a gas temperature $T = 500$K, and three PAH sizes.

## 5.2. Comparison of model results with observations for the sample of Galactic and extragalactic regions

### 5.2.1. Photoelectric heating efficiency

Figure 9 compares the variation with $\gamma$ of the observational diagnostic $R_e$ to the modeled values of $\epsilon_{PAH}$ for the set of objects listed in Table 2. Here we adopt $N_C = 54$ in the model, and a gas temperature $T = 500$ K, but we note that the choice of values for these parameters has a limited effect on the model results for $\gamma \gtrsim 10^3$ (see Figs. 10, 11, and 7), applicable to most sources. In contrast, since the effective temperature of the radiation field has strong effects on the value of $\epsilon_{PAH}$ (see Spaans et al. 1994), we compute the models with appropriate effective temperatures. Figure 9 shows that, in general, there is a good agreement between $R_e$ and $\epsilon_{PAH}$, in terms of curve shape but also of absolute values. The upper part of Fig. 9 presents the comparison between $R_e$ and $\epsilon_{PAH}$ for interstellar regions irradiated by the most massive stars with $20000 \le T_{\text{eff}} \le 40000$ K. In the Orion Veil, $R_e$ is somewhat lower than the predicted $\epsilon_{PAH}$, for a star with $T_{\text{eff}} = 40\,000$ K (which is appropriate for $\Theta^1$ Orionis C). This lower value can be explained by the fact that the radiation field impinging on the Veil is likely softer than that of $\Theta^1$ Orionis C, because of the UV absorption by dust and PAHs situated in between the Veil and the massive star. The Orion Bar which is much closer to $\Theta^1$ Orionis C has a value of $R_e$ which is in good agreement with that of $\epsilon_{PAH}$ for $T_{\text{eff}} = 40\,000$ K. The other nebulae illuminated by O or early B stars in the sample, that is the Horsehead and IC 63 also have values of $R_e$ in good agreement with modeled values of $\epsilon_{PAH}$ for $T_{\text{eff}} = 30\,000 - 40\,000$ K. For the LMC, $R_e$ is in good agreement with the $\epsilon_{PAH}$ curve for the most massive stars ($T_{\text{eff}} > 3 \times 10^4$K). This is compatible with the fact that the values we use to derive $R_e$ from Rubin et al. (2009) and Lebouteiller et al. (2012) concern star-forming regions. $R_e$ in star-forming galaxies falls within the range of the models with $T_{\text{eff}} \sim 25\,000$ K, implying an average radiation field which is softer than that of O stars and dominated by B-type stars. $R_e$ for ULIRGs is also in good agreement with the model. ULIRGs have a lower heating efficiency as compared to star-forming galaxies of the local universe, which is well explained
in the model by the higher value of the ionization parameter $\gamma$ in these galaxies due to intense massive star formation, increasing the average integrated intensity of the radiation field $G_0$ (Díaz-Santos et al. 2017; McKinney et al. 2021). The $R_e$ value of the high redshift galaxy GS IRS20 is also in good agreement with the model. Figure 9 also presents the comparison between $R_e$ and $\epsilon_{\text{PAH}}$ for interstellar regions irradiated by intermediate-mass stars with $10000 \le T_{\text{eff}} \le 15000$ K. This includes NGC 7023 which was presented in detail in Sect. 5.1. The case of Ced 201 is particularly interesting, because it has a low PAH heating efficiency, which is well explained by the fact that this nebula is illuminated by a B 9.5 star, whose cooler spectrum is less efficient to ionize PAHs (see Sect. C.4). The value of $R_e$ for Ced 201 indeed falls close to the model value of $\epsilon_{\text{PAH}}$ for a $T_{\text{eff}} = 10000$ K, while a B 9.5 spectral type has $T_{\text{eff}} \sim 11000$ K. The lowest efficiencies are observed in protoplanetary disks (Fig. 9). In this case too there is a good agreement between $R_e$ and the modeled value of $\epsilon_{\text{PAH}}$, which are low due to both the high $\gamma$ values and cool spectral types (A and B) for these sources. Table 4 gives the values derived for $\epsilon_{\text{PAH}}$ for diffuse ISM conditions, that is $T = 80$ K and $n_e = 3 \times 10^{-2}$ cm$^{-3}$ for three classical radiation fields. Observations yield $R_e = 2.1\%$ (Table 2 and Appendix B), close but about a factor of ~2-3 lower than the model values, that is $\epsilon_{\text{PAH}} = 3.7, 5.9$ and $3.9\%$ for the interstellar radiation fields of Draine (1978), Mathis et al. (1983) and Habing (1968), respectively.

### 5.2.2. Heating rates

We compute the heating rates in sources for which there is a good estimate of the column density of the cooling gas, that is the Orion Bar (Salgado et al. 2016) and a diffuse ISM line of sight (Boulanger et al. 1996). For the Orion bar, with $I_{\text{gas}} = 6.6 \times 10^{-5}$ W m$^{-2}$ (Appendix B) and $N_{\text{H}} = 7 \times 10^{22}$ cm$^{-2}$ (Salgado et al. 2016), $\Lambda_{\text{gas}} = 1.2 \times 10^{-30}$ W H$^{-1}$. This is close to the modeled gas heating rate by PAHs, that is $\Gamma_{\text{PAH}} = 4.2 \pm 2.0 \times 10^{-30}$ W H$^{-1}$, considering $G_0 = 2 \times 10^4$ (Joblin et al. 2018) and $f_c = 0.01$ (Castellanos et al. 2014) and the range of $\gamma$ values in Table 2. For the diffuse ISM, Boulanger et al. (1998) find a cooling rate of $\Lambda_{\text{gas}} = 3 \times 10^{-33}$ W H$^{-1}$ from FIRAS observations, in good agreement with the model values of $\Gamma_{\text{PAH}}$ presented in Table 4.

# 6. Discussion

## 6.1. The role of PAHs in neutral gas heating

The comparison presented in Sect. 5 shows that there is a good agreement between the values of $\epsilon_{\text{PAH}}$ derived from the PAH heating model and the values of the observational diagnostics $R_e$, for a wide range of physical conditions and environments. In addition, there is an excellent agreement in terms of the spatial variations of the charge state of PAHs in NGC 7023. The agreement between the modeled heating rates by PAHs and the gas cooling rates derived from the observations is also good for NGC 7023, the Orion Bar, and the diffuse ISM. Overall, this suggests that photoelectric heating by PAHs alone can explain the neutral gas heating in a number of environments. While this has been suggested by several studies (see Introduction), a direct assessment comparing the results from a molecular model including state-of-the-art molecular parameters for PAHs to observational diagnostics, had not been conducted before. BT94 pointed out the importance of PAHs, however 50% of the gas heating in their model is due to grains larger than $N_C = 1500$, while the present PAH model contains only PAHs of $N_C = 54$ C

atoms. This stresses the central role that PAHs play in the physics of the gas heating of galaxies.

## 6.2. Heating by dust vs heating by PAHs

Pabst et al. (2021) and earlier works have compared the theoretical photoelectric heating efficiency of dust from BT94, $\epsilon$, to $R_{\text{IR}} = \frac{I_{\text{gas}}}{I_{\text{IR}}}$, where $I_{\text{IR}}$ is the total power radiated by dust in the far infrared. They found that $\epsilon > R_{\text{IR}}$ by a factor ~3 or more. This discrepancy is likely due to the fact that $R_{\text{IR}}$ cannot easily be compared to $\epsilon$, for two reasons. The first one is that, in the definition of BT94, $\epsilon$ is computed by integrating the power absorbed (emitted) by dust in the 6-13.6 eV range, while some fraction of the observed radiated power $I_{\text{IR}}$ is in reality due to absorption of UV photons with energies below 6 eV, and therefore $I_{\text{IR}}$ overestimates the denominator value. Still, it appears unlikely that this effect only can explain this factor of 3 difference, because most of the energy absorbed by dust does lie above 6 eV. A most likely reason is that $I_{\text{IR}}$ includes an important fraction of emission from dust grains which do not contribute at all to photoelectric heating.

## 6.3. PAH heating vs X-rays, cosmic rays, shocks, and turbulence in neutral gas.

Since the gas heating rate by PAHs $\Gamma_{\text{PAH}}$ depends on their abundance (Eq. 29), in regions where they are depleted, such as dense molecular clouds where PAHs are in condensed form (Rapacioli et al. 2006), other mechanisms are likely to take over, such as heating by cosmic rays (e.g. Padovani et al. 2009). In regions where the UV radiation field is weak compared to other sources of radiation such as X-rays (e.g. in the inner regions of protoplanetary disks near T-Tauri stars, or near active galactic nuclei), heating by PAHs may also become negligible. In PDRs, heating by PAHs is likely to dominate, but heating by H$_2$ can also become important in dense regions near the dissociation front (Bron 2014). In some regions, shocks may become an important source of neutral gas heating. The competition between heating by PAHs and other sources of gas heating can be tested using $R_e$: when this ratio is larger than the $\epsilon_{\text{PAH}}$ values predicted by the model, this is indicative that some other processes are at play. For instance, shock heating of neutral gas at the scale of galaxies results in extreme values of $R_e$ (e.g. >0.3, Appleton et al. 2013, Table 1), well above the maximum efficiencies derived with our model ($\epsilon_{\text{PAH}}^{\text{Max}} \sim 0.11$, see Table 5). This illustrates how the PAH model efficiencies can be used as a diagnostic of UV radiative feedback vs other sources of heating in galaxies.

# 7. Conclusion

In this article, we re-evaluate the contribution of PAHs to the photoelectric heating of neutral gas in astrophysical environments, using a simple analytical model in which state-of-the-art molecular parameters from laboratory measurements and quantum chemical calculations are included. We find that, for standard abundances of PAHs, the heating rates resulting from photoelectric heating by PAHs alone is comparable to those produced using a full size distribution of dust grains in the classical models of Bakes & Tielens (1994) and Weingartner & Draine (2001). The modeled and observed charge states of PAHs and their variation with physical conditions in NGC 7023 are in excellent agreement. The values of the photoelectric heating efficiency derived for PAHs from the model for a wide range of
Galactic and extragalactic regions are in good agreement with the values of the observational diagnostics $R_e$, that is the ratio of gas emission to the sum of gas + PAH emission. The values of the heating rate by PAHs derived from the model are also in good agreement with the observed gas cooling rates for the NGC 7023 nebula, the Orion Bar, and the diffuse ISM. Overall, this allows us to conclude that PAHs can explain most of the heating of the neutral gas in a variety of astrophysical environments, ranging from protoplanetary disks around young intermediate-mass stars to starburst galaxies in the early Universe. Our study highlights the importance to implement a robust description of the PE heating by PAHs in astrophysical codes computing the thermal balance of neutral gas, such as PDR models. The formalism presented here is included in a simple analytical model but it relies on our knowledge of the molecular properties of PAHs. This knowledge is crucial and laboratory studies (experiments and quantum chemical calculations must be guided by astronomical observations. In particular, the coming spectroscopic data obtained from the MIRI and NIRspec spectrometers on board the James Webb Space Telescope will provide unprecedented details on the PAH populations (Berné, Habart, Peeters, & the PDRs4All Team 2022).

*Acknowledgements.* OB acknowledges C. Pabst and J. R. Goicoechea for providing the data on the PE heating efficiency in Orion, J. Mc Kinney for providing the data regarding the PE heating efficiency in galaxies, F. Boulanger for fruitful discussions on the heating and cooling of the diffuse ISM, E. Josselin and A. Palacios for their help on extracting the stellar spectra used in this paper, and I. Schroetter for support on the writing of the code. The authors wish to thank L. Verstraete and V. Guillet for their comments on an early version of this work. We are grateful to the anonymous referee for his detailed report which greatly improved the manuscript. Finally, the research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) ERC-2013-SyG, Grant agreement N°610256 NANOCOSMOS. It has also been supported by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES, and through an APR grant provided by CNES.

## References

Acke, B. & van den Ancker, M. E. 2004, A&A, 426, 151

Alecian, E., Catala, C., Wade, G., et al. 2008, Monthly Notices of the Royal Astronomical Society, 385, 391

Allamandola, L. J., Tielens, A. G. G. M., & Barker, J. R. 1985, ApJ, 290, L25

Andrews, H., Candian, A., & Tielens, A. G. G. M. 2016, A&A, 595, A23

Andrews, H., Peeters, E., Tielens, A., & Okada, Y. 2018, A&A, 619, A170

Appleton, P. N., Guillard, P., Boulanger, F., et al. 2013, ApJ, 777, 66

Bakes, E. L. O. & Tielens, A. G. G. M. 1994, ApJ, 427, 822

Bauschlicher, C. W., Ricca, A., Boersma, C., & Allamandola, L. J. 2018, ApJS, 234, 32

Benisty, M., Perraut, K., Mourard, D., et al. 2013, A&A, 555, A113

Bernard-Salas, J., Habart, E., Arab, H., et al. 2012, A&A, 538, A37

Bernard-Salas, J., Habart, E., Köhler, M., et al. 2015, A&A, 574, A97

Berné, O., Habart, É., Peeters, E., & the PDRs4All Team. 2022. Publications of the Astronomical Society of the Pacific, 134, 054301

Berné, O., Montillaud, J., & Joblin, C. 2015, A&A, 577, A133

Berné, O., Mulas, G., & Joblin, C. 2013, A&A, 550, L4

Berné, O., Joblin, C., Deville, Y., et al. 2007, A&A, 469, 575

Boersma, C., Bregman, J., & Allamandola, L. J. 2016, ApJ, 832, 51

Boersma, C., Bregman, J. D., & Allamandola, L. J. 2013, ApJ, 769, 117

Boulais, A. 2013, Méthodes de séparation aveugle de source et application à l'imagerie hyperspectrale en astrophysique

Boulanger, F. 1999, in Solid Interstellar Matter: The ISO Revolution, 20

Boulanger, F., Abergel, A., Bernard, J.-P., et al. 1996, A&A, 312, 256

Boulanger, F., Boissel, P., Cesarsky, D., & Ryter, C. 1998, A&A, 339, 194

Bron, E. 2014, PhD thesis, Université Paris Diderot

Bréchignac, P., Garcia, G. A., Falvo, C., et al. 2014, The Journal of Chemical Physics, 141, 164325

Carelli, F., Grassi, T., & Gianturco, F. 2013, Astronomy & Astrophysics, 549, A103

Castellanos, P., Berné, O., Sheffer, Y., Wolfire, M. G., & Tielens, A. G. 2014, ApJ, 794, 83

Chokshi, A., Tielens, A. G. G. M., Werner, M. W., & Castelaz, M. W. 1988, ApJ, 334, 803

Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103

de Jong, T. 1977, A&A, 55, 137

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32

Draine, B. T. 1978, ApJS, 36, 595

Draine, B. T. 1978, ApJS, 36, 595

d'Hendecourt, L. & Léger, A. 1987, Astron. Astrophys, 180, L9

Finkenzeller, U. 1985, A&A, 151, 340

Fuente, A., Martin-Pintado, J., Neri, R., Rogers, C., & Moriarty-Schieven, G. 1996, A&A, 310, 286

Habart, E., Abergel, A., Walmsley, C., Teyssier, D., & Pety, J. 2005, A&A, 437, 177

Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421

Helou, G., Malhotra, S., Hollenbach, D. J., Dale, D. A., & Contursi, A. 2001, ApJLetters, 548, L73

Hillier, D. J. & Miller, D. L. 1998, ApJ, 496, 407

Iida, S., Kuma, S., Kuriyama, M., et al. 2021, Phys. Rev. A, 104, 043114

Joblin, C., Bron, E., Pinto, C., et al. 2018, A&A, 615, A12

Joblin, C., Pilleri, P., Montillaud, J., et al. 2010, A&A, 521, L25

Joblin, C., Wenzel, G., Castillo, S. R., et al. 2020, J. Phys. Conf. Ser., 1412, 062002

Jochims, H. W., Baumgaertel, H., & Leach, S. 1996, A&A, 314, 1003

Jura, M. 1976, ApJ, 204, 12

Kurucz, R. L. 1993, VizieR Online Data Catalog, VI/39

Le Petit, F., Nehme, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506

Lebouteiller, V., Cormier, D., Madden, S. C., et al. 2012, A&A, 548, A91

Leger, A., D'Hendecourt, L., & Defourneau, D. 1989, A&A, 216, 148

Leger, A. & Puget, J. L. 1984, A&A, 500, 279

Lemaire, J. L., Field, D., & Maillard, J. P. 1996, A

Lepp, S. & Dalgarno, A. 1988, ApJ, 335, 769

Malloci, G., Joblin, C., & Mulas, G. 2007, Chemical Physics, 332, 353

Malloci, G., Mulas, G., Cappellini, G., Fiorentini, V., & Porceddu, I. 2005, A&A, 432, 585

Malloci, G., Mulas, G., & Joblin, C. 2004, A&A, 426, 105

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 500, 259

McKinney, J., Armus, L., Diaz-Santos, T., et al. 2021, arXiv preprint arXiv:2101.01182

McKinney, J., Pope, A., Armus, L., et al. 2020, ApJ, 892, 119

Meeus, G., Montesinos, B., Mendigutía, I., et al. 2012, A&A, 544, A78

Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2015, 151, 5

Montillaud, J., Joblin, C., & Toublanc, D. 2013, A&A, 552, A15

Mulas, G., Malloci, G., Joblin, C., & Toublanc, D. 2006, A&A, 456, 161

Okada, Y., Pilleri, P., Berné, O., et al. 2013, A&A, 553, A2

Pabst, C., Goicoechea, J., Hacar, A., et al. 2021, arXiv preprint arXiv:2111.12363

Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619

Palacios, A., Gebran, M., Josselin, E., et al. 2010, A&A, 516, A13

Pilleri, P., Montillaud, J., Berné, O., & Joblin, C. 2012, A&A, 542, A69

Racine, R. 1968, AJ, 73, 233

Rapacioli, M., Calvo, F., Joblin, C., et al. 2006, A&A, 460, 519

Rapacioli, M., Joblin, C., & Boissel, P. 2005, A&A, 429, 193

Ricca, A., Bauschlicher, C. W., Boersma, C., Tielens, A. G. G. M., & Allamandola, L. J. 2012, The Astrophysical Journal, 754, 75

Rubin, D., Hony, S., Madden, S., et al. 2009, A&A, 494, 647

Salas, P., Oonk, J., Emig, K., et al. 2019, A&A, 626, A70

Salgado, F., Berné, O., Adams, J. D., et al. 2016, ApJ, 830, 118

Smith, J. D. T., Armus, L., Dale, D. A., et al. 2007a, PASP, 119, 1133

Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007b, ApJ, 656, 770

Spaans, M., Tielens, A. G. G. M., van Dishoeck, E. F., & Bakes, E. L. O. 1994, ApJ, 437, 270

Spitzer, L. 2004, Physical processes in the interstellar medium, Physics textbook (Wiley)

Spitzer Jr, L. 1948, ApJ, 107, 6

Terenin, A. & Vilessov, F. 1964, Advances in Photochemistry, 2, 385

Tielens, A. 2005, The Physics and Chemistry of the Interstellar Medium

Tielens, A. G. 2021, Molecular astrophysics (Cambridge University Press)

Tielens, A. G. G. M. 2008, Annual Review of A&A, 46, 289

Tschurl, M. & Boesl, U. 2006, International Journal of Mass Spectrometry, 249-250, 364, chava Lifshitz Memorial Issue

Verstraete, L., Leger, A., D'Hendecourt, L., Defourneau, D., & Dutuit, O. 1990, A&A, 237, 436

Visser, R., Geers, V., Dullemond, C., et al. 2007, A&A, 466, 229

Wahab, A., Pfuderer, L., Paenurk, E., & Gershoni-Poranne, R. 2022, Journal of Chemical Information and Modeling, 62, 3704, pMID: 35881922

Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296

Weingartner, J. C. & Draine, B. T. 2001, ApJS, 134, 263

Wenzel, G., Joblin, C., Giuliani, A., et al. 2020, A&A, 641, A98

Werner, M. W., Uchida, K. I., Sellgren, K., et al. 2004, ApJS, 154, 309

Zhen, J., Castillo, S. R., Joblin, C., et al. 2016, ApJ, 822, 113
# Appendix A: Derivation of $\gamma$ and $R_e$ in NGC 7023

## Appendix A.1: Derivation of $\gamma$ in the NGC 7023 NW PDR

In the following, we describe how we derived values for the physical conditions, that is $G_0$, $T$ and $n_e$, necessary to determine $\gamma$ in NGC 7023, for 8 angular distances ranging from 20 to 60" from the star. To derive the radiation field intensity, we used the procedure described in Pilleri et al. (2012) and adopted in Joblin et al. (2018). Briefly, this consists in using the sum of two identical synthetic spectra with effective temperatures $1.5 \times 10^4$ K from the Kurucz (1993) library and a radius of 10 $R_\odot$ for the HD 200775 stars (Aleczian et al. 2008). An additional extinction of $A_V = 1.5$ was applied to account for the presence of dust in the immediate surrounding of the stars. The radiation field intensity as a function of distance, in units of the Habing field, is derived assuming that the star-to-PDR inclination is of 63 degrees, which is necessary to reconcile the sky projected distance to the modeled star-to-PDR distance of 0.143 pc (Joblin et al. 2018). For angular distances from the star above 45", the effect of dust absorption in the PDR is taken into account using a factor $e^{-\tau(\lambda)}$ where $\tau(\lambda) = \sigma_{ext}(\lambda)N_H$, with $N_H$ the total hydrogen column density along the axis between the PDR and the star. With a constant density in the PDR of $n_H = 2 \times 10^4 \text{ cm}^{-3}$ (Pilleri et al. 2012), the estimated values for $N_H$ at positions situated 50", 55" et 60" from the star are $10^{21}$, $2 \times 10^{21}$, and $3.2 \times 10^{21} \text{ cm}^{-2}$, respectively. This corresponds to extinction values of $A_V = 0.5$, 1, and 1.6 mag, respectively. The resulting radiation field intensities are computed by integrating the specific intensities in the 5.17 to 13.6 eV range, and are reported in Table 1. The gas temperature values for the positions up to 45" from HD 200775 are taken from the study of Boulais (2013) who modeled the spatial profile of the [C II] 158 µm line emission using an LTE radiative transfer model to determine a stable temperature of ~ 750 K. The values adopted inside the PDR, that is above 45", are from the models of Montillaud et al. (2013). Gas densities $n_H$ in the cavity region up to 25" from the star are from Berné et al. (2015). The authors provided a value of $n_H \sim 1.5 \times 10^4 \text{ cm}^{-3}$ at ~ 25", which we assume to be constant up to the PDR dissociation front from which we use a value of $n_H \sim 2 \times 10^4 \text{ cm}^{-3}$ (Pilleri et al. 2012). The gas temperature and density for all studied positions are reported in Table 1. The electron density is derived from the gas density by considering that all electrons originate from the ionization of carbon, that is $n_e = 1.6 \times 10^{-4} n_H \text{ cm}^{-3}$. The gas temperature and the electron density can then be used to derive the ionization parameter $\gamma = \frac{G_0 \sqrt{T}}{n_e}$ (see values in Table 1). Andrews et al. (2016) used $\gamma \sim 2 \times 10^4 \text{ K}^{1/2} \text{cm}^3$ at the border of the PDR, (~ 42" from HD 200775), which is consistent with the values derived in the present study.

## Appendix A.2: Derivation of $R_i$ in the NGC 7023 NW PDR

We limit the study to the field of view (FoV) presented in Fig. 1, which excludes the parts of the nebula that are closer than 20" from the star, and farther than 60". The former have to be excluded because other features attributed to species such as $C_{60}$ and $C^+_{60}$ (Berné et al. 2013), and a strong thermal continuum due to hot dust grains (Berné et al. 2007) are present in the spectra. This contamination makes a clean analysis of the PAH emission with PAHTAT difficult. Similarly, at distances above 60", the signal-to-noise ratio in the spectra becomes too low to use PAHTAT. In addition, regions too far from the star become optically thick in the mid IR (see Fig. 4 in Pilleri et al. 2012), which is also a challenge when using PAHTAT. We used the dataset

obtained by Werner et al. (2004) with the IRS instrument aboard Spitzer. The data were previously analyzed by Berné et al. (2007) and Boersma et al. (2013), notably. The cube was obtained at a spatial resolution of 3.6" with a sampling of pixels of 1.8" and a spectral resolution of ~ 80, over a spectral range ranging from 5.5 to 15 µm. We applied the PAHTAT tool (Pilleri et al. 2012), which adjusts spectral templates of PAHs (notably cationic and neutral PAHs) to observations, to all spectra in the mid-IR cube to derive the integrated intensities of the $PAH^+$ and $PAH^0$ components (respectively $I_{PAH^+}$ and $I_{PAH^0}$) as well as their sum $I_{PAH} = I_{PAH^+} + I_{PAH^0}$. The $R_i$ map is then obtained using $I_{PAH^+}$ and $I_{PAH^0}$ following Eq. 1.

## Appendix A.3: Derivation of $R_e$ in the NGC 7023 NW PDR

Gas cooling in the neutral ISM of galaxies is dominated by emission in fine-structure lines of [C II] at 158 µm, and of [O I] at 63 µm and 145 µm. In addition, emission in the rotational lines of H$_2$, in the mid-IR, and emission in CO lines and other molecular species in the far IR can contribute. Bernard-Salas et al. (2015) provided a detailed study of gas cooling in the far IR in NGC 7023 with Herschel-PACS and SPIRE. We use their data for the [C II] line at 158 µm, and the [O I] lines at 63 and 145 µm (Fig. A.1). We complement them with the 0-0 S(1), S(2), and S(3) rotational lines of H$_2$ at 17.0, 12.3, and 9.7 µm derived from Spitzer observations (Fig. A.1). The emission maps for the S(2) and S(3) H$_2$ lines were derived using PAHTAT (Pilleri et al. 2012) on the Spitzer-IRS SL data cube of the NGC 7023 NW PDR obtained by Werner et al. (2004). The map of the S(1) line integrated intensity was extracted using the CUBISM software (Smith et al. 2007a) on the LL data cube obtained by Werner et al. (2004). Overall, the gas line emission is largely dominated by the [O I] and [C II] lines. H$_2$ emission is important only at the dissociation front (up to 30% maximum of the total line emission), but overall it represents less than 10% of the gas cooling. In addition, cooling by H$_2$ emission occurs in regions where gas heating is not dominated by PE heating, but rather by the formation of H$_2$, and the collisional de-excitation of vibrationally pumped H$_2$. Hence, we do not include H$_2$ emission in the gas cooling budget, $I_{gas}$ being then the sum of [O I] and [C II] lines emission. As for the derivation of the map of $R_i$, we limit the study to the FoV presented in Fig. 1, which excludes the parts of the nebula that are closer than 20" from the star, and farther than 60". The $R_e$ map is then obtained from $I_{PAH}$ (Sect. A.2) and $I_{gas}$ following Eq. 2.

# Appendix B: Derivation of $R_e$ and $\gamma$ in a sample of sources.

## Appendix B.1: Galactic photodissociation regions

We include data for the Orion Veil, which is the extended nebular region in front and south of the M42 region of Orion A, illuminated by the O6 star $\Theta^{1}$ Ori-C. It is a diffuse nebula where cooling is largely dominated by [CII] emission (Pabst et al. 2021). All the data used here were kindly provided by C. Pabst. It consists in the table of values of $\gamma$ obtained for this region totaling over 11 000 pixels, as well as the values of the [CII] line intensities $I_{[CII]}$ derived from SOFIA observations. In addition, we use the Spitzer-IRAC 8 µm maps of the region provided by Megeath et al. (2015). Emission in this IRAC filter ($I_{8\mu m}$) is largely dominated by PAH emission, hence, in the absence of spectroscopic measurements, it can be used to evaluate $I_{PAH}$. We have used the data of NGC 7023, where $I_{PAH}$ is obtained from spectroscopic
<figure><img src="image_14.png" /><figcaption>Fig. A.1: Emission maps of gas cooling. Maps on the upper row show the [C II] line at 158 µm, as well as the [O I] lines at 145 µm and 63 µm, which were observed with Herschel-PACS at a resolution of 11", 8.8", and 4.5" respectively (Bernard-Salas et al. 2015). Maps on the lower row show the emission from the H₂ rotational S(1), S(2) et S(3) lines at 17.035, 12.278, and 9.7 µm, respectively. These maps were extracted from the Spitzer archival data using CUBISM Smith et al. (2007a). The S(2) and S(3) maps are from the SL module of IRS (resolution of 3.6"), and the S(1) map is from the IRS LL module (resolution of 10"). Emission from the H₂ v=1-0 S(1) line at 2.121 µm is presented in contours as a spatial reference for the dissociation front, but it is not included in the cooling budget, since this line is pumped by UV photons and hence does not cool the gas.</figcaption></figure>

analysis (Sect. A.2) as well as the $I_{8\mu m}$ IRAC maps, and found $I_{PAH} \approx 1.15 \times 10^{-7} I_{8\mu m}$, with $I_{PAH}$ in W m$^{-2}$ sr$^{-1}$ and $I_{8\mu m}$ in MJy sr$^{-1}$. This is used to derive $I_{PAH}$ from the IRAC data, and $R_e$ is then given by $\frac{I_{[CII]}}{I_{PAH}+I_{[CII]}}$. The values presented in Fig. 3 and summarized in Table 2 are obtained by determining the median and first and 99th percentiles for $R_e$, and median and 25th and 75th percentiles for $\gamma$. The used interval for $\gamma$ is smaller than for $R_e$ because of the presence of a larger number of outlier values of $\gamma$.

We include data for the Orion Bar (also illuminated by $\Theta^1$ Ori-C) from Bernard-Salas et al. (2012) for the cooling lines ([OII], [CII]) at position ($\alpha = 5:35:19.778, \delta = -5:25:30.65$). We derive a total intensity for the cooling lines (dominated by [OI] at 63 µm) of $6.6 \times 10^{-5}$ W m$^{-2}$ sr$^{-1}$. We extract the value of $I_{8\mu m}$ at this position both from the IRAC data and available IRS spectroscopy, which respectively yields $I_{8\mu m} = 5500$ and $I_{8\mu m} = 7850$ MJy sr$^{-1}$, which gives $I_{PAH} = [6.32-9.02] \times 10^{-4}$ W m$^{-2}$ sr$^{-1}$. $R_e$ is found to fall in the range [0.07-0.11]. The range of $\gamma$ values is derived from model results for the Orion Bar presented in Joblin et al. (2018), and corresponding to the range of $\Lambda_V$, for which the calculated [OI] and [CII] emissions are the dominant cooling lines. Values of $R_e$ and $\gamma$ for the Orion Bar are summarized in Table 2.

We include the IC 63 nebula, which is illuminated by the B0.5 star $\gamma$-Cas, and for this source $I_{PAH}$ is taken from Table 3 in Pilleri et al. (2012). The values for the cooling lines are taken from Andrews et al. (2018). The error on absolute intensities of

the cooling lines is ~30%, including uncertainty on beam filling factor and instrumental error. The values of $\gamma$ are from Pabst et al. (2021).

Ced 201 is a nebula formed by a runaway B9.5 star penetrating a molecular cloud in Cepheus. For this source, $I_{PAH}$ is also taken from Pilleri et al. (2012). The values for the cooling lines intensities are directly taken from the observations available in the Herschel Science Archive and extracted on the central spaxel of the PACS instrument. The same error as for IC 63 on the intensity, that is 30%, is considered. The values of $\gamma$ are from Pabst et al. (2021).

The Horsehead is illuminated by the $\sigma$-Ori AB stars of spectral type O9.5V and B0.5V. The values for $I_{PAH}$ are also taken from Pilleri et al. (2012), and the values for the cooling lines are also from the observations available in the Herschel Science Archive and extracted in the central spaxel of PACS (error = 30%). To estimate the values of $\gamma$ we use a FUV radiation field of $G_0 = 100$, a temperature of $T = 100-500$ K, and a gas density of $n_H = 1-20 \times 10^4$ cm$^{-3}$ (Habart et al. 2005). This yields $\gamma$ values in the range [30-1400].

## Appendix B.2: The diffuse interstellar medium

For the diffuse interstellar medium, we use the IR emission intensity for PAHs deduced by Compiègne et al. (2011), that is $I_{PAH} = 1.45 \times 10^{-31} \text{W H}^{-1}$ and the cooling measured for the [CII] line from Boulanger et al. (1996), that is $\Lambda_{\text{gas}} = 3 \pm$
0.4 × 10⁻³³W H⁻¹, which yields Rₑ = 0.021 ± 0.002. We neglect the cooling from [OI] which is marginal for this low density medium. The average gas density ranges between 50 and 60 cm⁻³ (Boulanger 1999), and the gas temperature is T ~ 80 K. The radiation field in the diffuse ISM is G₀ ~ 1.7. With nₑ = 1.6 × 10⁻⁴n_H, this yields γ ∈ [1580 – 1900].

## Appendix B.3: The Large Magellanic Cloud

For the Large Magellanic Cloud (LMC), we rely on the work of Lebouteiller et al. (2012) who provide a value of Rₑ = 0.07, and quote a factor of 2 uncertainty. We thus adopt Rₑ ∈ [0.05 – 0.01]. From Fig. 10 of Rubin et al. (2009), we derive G₀/nₑ ∈ [630 – 1260], which for a temperature of 75 K (Rubin et al. 2009) yields γ ∈ [5460 – 10 900].

## Appendix B.4: Star-forming galaxies

For star-forming galaxies in the local Universe, we use the ratio of the [CII] emission to the integrated flux in the 5 to 10 µm range provided in Helou et al. (2001). Following McKinney et al. (2021), we apply a factor 0.4 to this ratio to obtain Rₑ, to account for the fact that PAH emission is responsible for only about 40% of the mid-IR emission of galaxies (Helou et al. 2001). The range of values of γ for the galaxies are taken from McKinney et al. (2021) using the same approach as for ULIRGs, which is based on the relationship of Díaz-Santos et al. (2017) and presented in the following section.

## Appendix B.5: Ultra luminous infrared galaxies

McKinney et al. (2021) present measurements of the ratio between the PAH to [CII] emissions, which we adopt to derive Rₑ. These authors also provide the values of G₀/n_H in a large sample of ultra luminous infrared galaxies (ULIRGs). Using these data, kindly provided in electronic format by the authors, we derive the median values of Rₑ and γ and 25th and 75th confidence intervals (see values in Table 2). We compute γ from G₀/n_H considering a temperature of 100 K, in agreement with estimated values for T in these galaxies (Díaz-Santos et al. 2017), and using nₑ = 1.6 × 10⁻⁴ n_H.

## Appendix B.6: High-redshift galaxy GS IRS20

McKinney et al. (2020) provide a measurement of the ratio of the [CII] cooling line to the intensity of the PAH 6.2 µm band in the GS IRS20 galaxy at a redshift z ≈ 1.9. This value and associated error bar can be converted to Rₑ using the calibration factor of 0.11 obtained by Smith et al. (2007b) on a large sample of nearby galaxies for the ratio of the 6.2 µm band to total PAH emission. Additional uncertainty stems from the lack of precise measurement of the 11.2 µm feature which implies a possible overestimation of the total PAH emission, of a factor of ~2 (McKinney priv. com.), which we include to obtain the maximum value of Rₑ. The range of values for γ is obtained by using the observed value of the IR surface density of this galaxy given in McKinney et al. (2020), that is ΣIR ≈ 5 × 10¹¹ L☉kpc⁻², converted to a value of G₀/n_H using the law calibrated on ULIRGs by Díaz-Santos et al. (2017), and finally to a value of γ using a temperature of 100 K and nₑ = 1.6 × 10⁻⁴ n_H.

## Appendix B.7: Protoplanetary Disks

The disks for which we derived Rₑ are Herbig Ae/Be objects, which show prominent PAH emission. Amongst the disks observed by ISO and presented in Acke & van den Ancker (2004), we selected those for which there are Herschel measurements of the gas cooling lines (Meeus et al. 2012). This represents a total of 10 objects, which are listed in Table B.1. For some objects the [OI] 145 µm and [CII] lines were not detected. However this is not an issue since the [OI] 63 µm line is by far the most important cooling line in these objects, where the gas is dense and highly irradiated. The point presented in Fig. 3 represents the median value of Rₑ from Table B.1 and the 25th and 75th percentiles. The range of values for γ can only be obtained indirectly. We used the results of the detailed model for PAHs in Herbig Ae/Be stars developed by Visser et al. (2007). In particular, these authors show that over 80% of the PAH emission arises from regions beyond the inner 100 AU of the disk, and their Fig. 5 provides the values for G₀/nₑ, which range between 10³ and 10⁴ in this region of the disk. For a gas temperature of ~1000 K, this implies values for γ in the range [3 × 10⁴ – 3 × 10⁶]. The used temperature value is a rough estimate, however this choice is not critical since γ depends only on √T.

## Appendix C: Derivation of the ionization potentials, ionization yields and partition coefficient.

### Appendix C.1: Ionization potentials

Following Wenzel et al. (2020), the values of the ionization potentials as a function of the charge Z, for Z ≥ 0, are given by:

$$IP^{(Z+1)+} = 3.9 \text{ eV} + \frac{1}{4\pi\epsilon_0} \left[ \left(Z + \frac{1}{2}\right) \frac{e^2}{a} + \left(Z + 2\right) \frac{e^2}{a} \frac{0.03 \text{ nm}}{a} \right] \frac{1C}{e} \text{ eV}, \quad (C.1)$$

where ε₀ is the vacuum permittivity, e the elementary charge, and a the diameter of the molecule. The change of variable between a and the number of carbon atoms is given by (WD01):

$$a = \sqrt[3]{\frac{N_C}{468}} \text{ nm}. \quad (C.2)$$

For Z = -1, we use IP = 0 (but see Appendix C.3).

### Appendix C.2: Ionization yields

Several studies provide measurements of ionization yields for neutral PAHs (Jochims et al. 1996; Verstraete et al. 1990) and cationic PAHs (Zhen et al. 2016; Wenzel et al. 2020). From these studies, empirical laws have been derived to model the ionization yields. For neutral species (Z=0), we use the law derived by Jochims et al. (1996) (energy values in eV):

$$Y(E) = \begin{cases} \frac{E - IP^+}{9.2} & \text{for } \frac{E \le IP^+ + 9.2}{E > IP^+ + 9.2} \\ 1 & \text{for } E \le IP^+ + 9.2 \end{cases} \quad (C.3)$$
<table><caption>Table B.1: Intensities of gas cooling lines and PAH emission in Herbig Ae/Be protoplanetary disks.</caption><thead><tr><th>Object Name</th><th>[OI] 63 µm</th><th>[OI] 145 µm</th><th>[CII] 158 µm</th><th>PAHs</th><th>R<sub>e</sub></th></tr></thead><tbody><tr><td>AB AUR</td><td>851.2</td><td>44.60</td><td>51.00</td><td>17.9</td><td>0.0053</td></tr><tr><td>HD 97048</td><td>1592.0</td><td>65.6</td><td>106.8</td><td>55.9</td><td>0.0029</td></tr><tr><td>HD 100453</td><td>66.10</td><td>-</td><td>-</td><td>7.10</td><td>0.0009</td></tr><tr><td>HD 100546</td><td>6043.0</td><td>194.7</td><td>203.8</td><td>92.2</td><td>0.0070</td></tr><tr><td>HD 135344</td><td>47.90</td><td>-</td><td>-</td><td>1.09</td><td>0.0044</td></tr><tr><td>HD 141569</td><td>245.3</td><td>24.90</td><td>11.40</td><td>2.60</td><td>0.0108</td></tr><tr><td>HD 139614</td><td>44.50</td><td>-</td><td>-</td><td>10.4</td><td>0.0004</td></tr><tr><td>HD 142527</td><td>53.80</td><td>-</td><td>-</td><td>2.80</td><td>0.0019</td></tr><tr><td>HD 163296</td><td>208.4</td><td>-</td><td>-</td><td>4.10</td><td>0.0051</td></tr><tr><td>HD 169142</td><td>91.50</td><td>-</td><td>-</td><td>10.16</td><td>0.0009</td></tr></tbody></table>

**Notes.** Intensities for cooling lines are in $10^{-18}$ W m$^{-2}$ (from Meeus et al. 2012) and for PAHs in $10^{-14}$ W m$^{-2}$ (from Acke & van den Ancker 2004). $R_e$ values are calculated using Eq. 2.

For the second ionization ($Z = 1$), we use the law provided by Wenzel et al. (2020):

$$ Y^+(N_c, E) = \begin{cases} 0 & E < \text{IP}^{2+} \\ \frac{\alpha}{11.3-\text{IP}_{2+}}(E - \text{IP}^{2+}) & \text{IP}^{2+} \le E < 11.3 \\ \alpha & \text{for } 11.3 \le E < 12.9 \\ \frac{\beta(N_c)-\alpha}{2.1} (E - 12.9) + \alpha & 12.9 \le E < 15.0 \\ \beta(\bar{N}_c) & E \ge 15.0, \end{cases} \quad (C.4) $$

where $\alpha = 0.3$ and $\beta$ is given by:

$$ \beta(N_c) = \begin{cases} 0.59 + 8.1 \cdot 10^{-3} N_c & 32 \le N_c < 50 \\ 1 & N_c \ge 50. \end{cases} \quad (C.5) $$

Finally, for anions ($Z=-1$), we use $Y = 1$ following Visser et al. (2007).

### Appendix C.3: Threshold energy for electron detachment

In Sect. 3.1, we have assumed that $Y^-(E) = 1$ at all energies. This approximation can be justified by the low values of both the electronic affinity and the effective photoabsorption cross section at low energies (see Fig. 5). Recently, Iida et al. (2021) have shown that electron detachment from the pentacene anion, ($C_{22}H_{14}^-$), does not occur at the energy given by the electron affinity but is shifted to higher energies due to competition with radiative cooling. The authors estimated the survival rate of $C_{22}H_{14}^-$ to be 50% for an internal energy of 2.5 eV and 0% above 2.8 eV. We calculated the corresponding microcanonical temperature to be 775 and 824 K for 2.5 and 2.8 eV, respectively. The calculations use the density of states derived from a modified version of the Beyer-Swinehart algorithm (Mulas et al. 2006) and the harmonic vibrational frequencies from the *Theoretical PAH Spectral Database*. The electron affinity of circumcoronene is calculated to be similar to that of pentacene (see *Theoretical PAH Spectral Database*). We therefore assumed that efficient electron detachment would occur for temperatures higher than 800 K, which correspond to internal energies higher than 6 eV for $C_{54}H_{18}^-$. The calculations were performed with the same method as for $C_{22}H_{14}^-$ and using the harmonic frequencies from the *NASA Ames PAH IR Spectroscopic Database* (Bauschlicher et al. 2018; Ricca et al. 2012).

In Fig. C.1, we compare the model results for the PAH heating rate $\Gamma_{PAH}$ and the photoelectric heating efficiency of $\epsilon_{PAH}$ as a function of the ionization parameter $\gamma$, which were obtained using $Y^-(E) = 1$ for $E \ge 6$ eV (and $Y^-(E) = 0$ for $E < 6$ eV) compared to the case in which $Y^-(E) = 1$ for all energies. The values of the heating rate and efficiencies are lower by a factor ~1.5 at $\gamma = 10$, by a few percent at $\gamma = 100$, and are not affected at $\gamma = 100$ and above. This indicates that the choice of threshold for $Y^-(E)$ does not have a significant effect on the results of the model, especially for the large values of $\gamma$ that correspond to the objects studied in this paper (Fig. 9).

### Appendix C.4: Partition coefficient

A value of $\langle\gamma_e(E)\rangle = 0.5$ has been adopted in the literature (Verstraete et al. 1990; Tielens 2005; Bakes & Tielens 1994), based on the photoelectron kinetic spectra obtained for specific incident photon energies on benzene by Terenin & Vilesov (1964). We provide a revised estimate for this parameter based on the photoelectron spectroscopy measurements presented for the coronene molecule in Bréchignac et al. (2014). In particular, these authors provide (see their Fig. 2), for a given photon energy $E$ absorbed in the range from 7.3 eV (IP of coronene) to 10.5 eV, the probability density function of kinetic energy of the photoelectrons $E_K(E)$, from which $\langle\frac{E_K(E)}{E-IP}\rangle$ can be derived. Figure C.2 shows the values of $\gamma_e(E)$ for coronene, from which $\langle\gamma_e(E)\rangle$ is obtained by averaging over energies, yielding $\langle\gamma_e(E)\rangle = 0.46 \pm 0.06$ (68% confidence interval), in agreement with the classical value of 0.5.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Fig. C.1: Model derived heating rate Γ<sub>PAH</sub> (upper panel) and PAH photoelectric heating efficiency ε<sub>PAH</sub> (lower panel) obtained for IP<sup>0</sup> = 0 and 6 eV. In both cases, the model is computed with N<sub>C</sub> = 54, T = 500 K, and T<sub>Eff</sub> = 3 × 10<sup>4</sup> K. For the heating rate, we use f<sub>C</sub> = 0.1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. C.2: Fraction of the energy that goes into kinetic energy of the photoelectron, γ<sub>e</sub>(E), as a function of incident photon energy, derived from the photoelectron spectroscopic measurements of Bréchignac et al. (2014) for coronene, C<sub>24</sub>H<sub>12</sub>. The sharp features are induced by the presence of auto-ionization states.</figcaption></figure>