EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 664 (August 2022) A&A, 664 (2022) A67 Full HTML
Open Access
Issue
A&A
Volume 664, August 2022
Article Number A67
Number of page(s) 43
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202141854
Published online 09 August 2022

© M. Röllig and V. Ossenkopf-Okada 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

Numerical models of photodissociation regions (PDRs) are a valuable tool to model their chemical and thermodynamic structure and to simulate their line and continuum emission (Tielens & Hollenbach 1985; Hollenbach et al. 1991; Le Bourlot et al. 1993a; Sternberg & Dalgarno 1995; Kaufman et al. 1999; Le Petit et al. 2006; Röllig et al. 2007, 2013). This is particularly relevant since it is commonly assumed that all neutral atomic hydrogen gas and a large fraction of the molecular gas in galaxies lie in PDRs (Hollenbach & Tielens 1997).

The transition from the atomic to the molecular phase in the interstellar medium (ISM) is controlled by the balance between photo-dissociation by far ultraviolet (FUV) radiation and the recombination of atomic hydrogen on grain surfaces. Due to its nature as a line absorption process (Federman et al. 1979), photo-dissociation of H2 is subject to efficient shielding that leads to a sharp transition from H to H2. This is particularly true for dense gas (n > 105 cm−3). This transition has also been studied analytically and compared with corresponding numerical results (Sternberg 1988; Krumholz et al. 2008, 2009; McKee & Krumholz 2010; Sternberg et al. 2014, 2021; Bialy & Sternberg 2016; Bialy et al. 2017).

Numerical models of the physical and chemical conditions of the ISM compute the balance between the formation and destruction of the considered chemical species. The respective reaction rates depend on the local temperature which is the result of balancing local heating and cooling processes. Cooling is predominantly done by spectral line emission of energetically excited states of atomic and molecular species. In PDRs the heating is dominated by photo-electric heating, vibrational de-excitation of electronically excited H2, and cosmic ray heating. Both heating and cooling depend on the chemical abundances of the relevant species, which numerically asks for an iterative solution.

Astrochemical models of the ISM chemistry date back to the early 1950s (Bates & Spitzer 1951) and they have been used to explain the abundance of early observed molecules. They have been constantly extended since then and most modern astrochemical ISM models also simulate chemical kinetic processes on dust-grain surfaces with various degrees of complexity (Tielens & Hagen 1982; Hasegawa & Herbst 1993; Garrod 2013). The first species that has been recognized as almost exclusively formed in space on grain surfaces is molecular hydrogen (Gould & Salpeter 1963; Hollenbach & Salpeter 1971) and all PDR models have to include this reaction to yield physically reasonable results. For an extended review of grain surface models, we refer the reader to Cuppen et al. (2017).

The first PDR models were presented almost 50 yr ago to explain fine-structure emission observed with submillimeter instruments (Melnick et al. 1979; Russell et al. 1980, 1981; Stacey et al. 1983; Crawford et al. 1985). The earliest models (Glassgold & Langer 1974, 1975, Glassgold & Langer 1976; Langer 1976; Black & Dalgarno 1977; Clavel et al. 1978; de Jong et al. 1980) have been applied to low gas density n < 103 cm−3 illuminated by the ambient UV radiation. Higher gas density models with more intense radiation have been performed by Tielens & Hollenbach (1985), van Dishoeck & Black (1986, 1988), Black & van Dishoeck (1987), Sternberg & Dalgarno (1989, 1995), Le Bourlot et al. (1993a), and Draine & Bertoldi (1996). Some of these models have been further developed since then, reaching a high level of sophistication. Here, we can only name a few prominent ones: the PDR Toolbox1 (Kaufman et al. 1999, 2006; Pound & Wolfire 2008; Hollenbach et al. 2009) provides an extensive set of computed model grids and user-friendly tools to apply the model results to data. The Meudon PDR code2 (Le Petit et al. 2006; Goicoechea & Le Bourlot 2007; Le Bourlot et al. 2012; Bron et al. 2014) is a plane-parallel code with a particular focus on a detailed treatment of physical and chemical processes, namely radiative transfer and shielding processes. It computes the line intensities of tens of atomic and molecular species. The chemistry includes hundreds of species with reactions in gas phase and on surfaces. The code is publicly available together with an extensive online database of model grids and analysis tools. The KOSMA-τ group collaborates with both teams to integrate the model results into their online databases and tools.

The CLOUDY3 code is a spectral synthesis code designed to simulate conditions in interstellar matter under a broad range of conditions (Ferland et al. 2013, 2017). With an initial focus on the ionized gas, it has been extended to other regimes including PDRs (van Hoof et al. 2004; Shaw et al. 2005; Abel et al. 2005; Abel & Ferland 2006; Abel 2006). PDR model physics has also been included in a number of other numerical models. ProDiMo4, a model of protoplanetary disks, including gas-phase, X-ray, as well as UV photo-chemistry, gas heating and cooling balance, and disk structure and (dust and line) radiative transfer, is actively developed and extended (Woitke et al. 2009, 2016; Kamp et al. 2010, 2017; Thi et al. 2019, 2020). Other disk models include, for instance, DALI, which is a thermo-chemical model of protoplanetary disk atmospheres (Bruderer et al. 2009b,a, 2010, 2012; Bruderer 2013). A general model for the thermochemical structure of disks presented by Gorti & Hollenbach (2004) and Gorti & Hollenbach (2008) was recently updated to include the effects of gas-grain chemistry in protoplanetary disks using a three-phase approximation (Ruaud & Gorti 2019). Agündez et al. (2018) developed a generic model of a protoplanetary disk, with a focus on the photochemistry, and they applied it to a T Tauri and a Herbig Ae/Be disk.

UCL_PDR5 is a publicly available, 1D, and time-dependent PDR model (Bell et al. 2005, 2006; Priestley et al. 2017 that has been extended into the first 3D PDR code 3D–PDR6 (Bisbas et al. 2012; Gaches et al. 2019). The KOSMA-τ team used their clumpy PDR model as building blocks for KOSMA-τ-3D, to model arbitrary 3D geometries (Cubick et al. 2008; Andree-Labsch et al. 2017; Yanitski, in prep.).

Most numerical models assume stationary conditions and do not solve the time dependence; however, the importance of dynamical effects in the ISM lead to a number of new nonstationary models. A first approach was to couple (magneto) hydro-dynamical (MHD) models with PDR models to post-process their output (e.g., Levrier et al. 2012; Bisbas et al. 2015) and as subgrid models (Li et al. 2018). With the recent growth of computing power, MHD models started to numerically solve PDR physics and small chemical networks in parallel (Nelson & Langer 1997; Glover et al. 2010; Grassi et al. 2014; Walch et al. 2015; Seifried & Walch 2016; Girichidis et al. 2016; Valdivia et al. 2016; Seifried et al. 2017; Franeck et al. 2018; Inoue et al. 2020; Bisbas et al. 2021). Some PDR codes solve the PDR physics and chemistry self-consistently with the dynamics, but in a reduced dimensionality (Kirsanova et al. 2009, 2020; Akimkin et al. 2015, 2017; Bron et al. 2018).

Because of the high model complexity, the large variety in model assumptions, and numerical approximations, it remains difficult to compare results of different PDR model codes. The PDR benchmark was a first attempt to cross-calibrate different PDR model codes against each other (Röllig et al. 2007). The KOSMA-τ model is one of the models that participated in the benchmark study. It is the only PDR model employing a spherical model geometry and it has been applied to a large variety of studies in the last 20 yr.

PDR models are constantly improved and calibrated against observations of Galactic and extragalactic PDRs. Major constraints were obtained in recent years by a growing number of observations of the main PDR cooling lines, in particular the fine-structure lines of [C II] and [O I]. Modern instruments such as (up)GREAT (German REceiver for Astronomy at Terahertz Frequencies Heyminck et al. 2012; Risacher et al. 2018) on board the Stratospheric Observatory for infrared Astronomy (SOFiA, Young et al. 2012) allow one to observe Galactic PDRs with high spatial and very high spectral resolution showing complex line profiles due to source kinematics, optical thickness effects, and foreground absorption (e.g., Pabst et al. 2019; Kirsanova et al. 2019; Schneider et al. 2021; Kabanovic et al. 2022). Similar studies for the Large Magellanic Clouds (LMC; Okada et al. 2019a) together with the first detection of velocity resolved [13C II] emission in the LMC (Okada et al. 2019b) show the potential of these observations for extragalactic studies (other examples include Röllig et al. 2016; Cormier et al. 2019; Madden et al. 2020; Bigiel et al. 2020; Tarantino et al. 2021). Future PDR models may also be applied to Atacama Large Millimeter Array (ALMA) observations of protostars and hot cores and corinos which reveal a rich complex organic chemistry (COM) that current PDR models have yet to explain. Data from future instruments, such as the James Webb Space Telescope (JWST, Gardner et al. 2006) and FYST (Parshley et al. 2018), will challenge present PDR model predictions. The near- and mid-infrared (iR) instruments on board of JWST allow for very detailed studies of the hot gas and dust on PDR surfaces and the detailed knowledge of the H2 excitation and the PAH energetics (Okada et al. 2013) will help to further calibrate these models (Berné et al. 2022). Consequently, an active development of sophisticated PDR models is essential to deal with the wealth of data from current and near-future observations. The KOSMA-τ model will try to contribute with its particular strengths.

In this paper we report on recent major developments of the model and give an updated overview of the model properties. In particular, we present an upgrade of the model chemistry in KOSMA-τ to account for the full grain-surface chemistry in a quasi-three-phase model allowing us to study the interplay between ice formation on grains and PDR physics.

2 Models of PDRs – general properties

Figure 1 shows the general structure of a PDR. Attenuation of FUV photons produces a chemical stratification with the more stable chemical species more abundant closer to the surface of the PDR, while species that are easily dissociated can only survive if sufficiently shielded from the FUV field. The chemical stratification is accompanied by a strong temperature gradient from the surface to the shielded parts of the PDR because the main heating takes place via ejection of photo-electrons from dust grains and subsequent collisional energy transfer from the ejected hot electrons to the gas particles (Bakes & Tielens 1994; Weingartner & Draine 2001a). A decreased FUV flux reduces the photo-electric heating (PEH) significantly and results in lower temperatures. FUV attenuation similarly affects heating by vibrational de-excitation of H2 which dominates over PEH at higher densities (Röllig et al. 2006). in this paper we denote total FUV field strength with x in units of the Draine field χD = 2.6 × 10−3 erg s−1 cm−2 (integrated from 91.2 to 200 nm) (Draine 1978)7.

Figure 2 shows the numerical scheme any PDR model has to solve. it can be divided into four distinct problems (gray boxes). Three of them are local but the fourth one couples them nonlocally, asking for a multilayer iteration scheme. The local problems are: 1a) the local chemical problem, that is solving the chemical balance equations at a given position in the PDR. 1b) the local excitation, that is solving the energy level population of all species relevant for cooling and for comparison with observations at a given position in the PDR. 1c) the local energy balance, that is computing the gas and dust temperature at a given position in the PDR. Those problems are mutually coupled asking for an iterative local solution. As they have to be solved at every spatial point of the model, they are embedded in a spatial loop. The fourth problem is 2), the solution of the radiative transfer (RT) equation. This step provides quantities such as the photo-dissociation rates of important species, the FUV intensity across the model cloud, and others. Step 2) nonlocally couples all positions within a model PDR. Each individual point requires the solution of all other points. Hence a global iterative solution is needed. Problems 1a)–1c) are solved iteratively at every spatial location in the cloud. Using this set of solutions the radiative transfer problem 2) is solved over the whole cloud geometry. As a consequence, the local physical conditions computed in the previous iteration will be updated due to revised absorption and emission of IR and UV photons, which results in modified local cooling and heating capabilities leading to a new temperature solution and therefore a new chemical solution. Global iteration of the scheme (1a)–1c)–2)) is then performed until a predefined convergence criterion is met. An exception to this scheme are semi-infinite plane-parallel models. In these models, radiation can only enter and escape from one side of the cloud. As a consequence, the UV transfer is computed as the code advances along the spatial grid without the need for separate loops.

Any other model geometry with the possibility of radiation entering and escaping from different directions requires repeated iterations on the spatial grid. The global as well as the local numerical iterations are stopped after predefined convergence criteria are met.

Typical input quantities of PDR models are summarized in Table 1. Table 2 lists typical output quantities of PDR models.

A comparison with experimental data is usually limited to integrated model quantities, such as intensities or emission line ratios. The limited spatial resolution of IR and FIR telescopes and the generally long distance to massive star forming regions allows only in very few cases to directly compare the modeled PDR structure to real observations. In terms of model calibration this is a substantial complication.

In the following we present the KOSMA-τ code. Section 3 lays out the details of the code structure and the implemented physics. Section 4 discusses the details of the numerical solution of the chemical problem. Section 5 discusses the energy solution and the occurrence of multiple solutions. Section 6 presents the new model predictions of the surface chemistry and compares with predictions by other models and with observations. Section 7 finally demonstrates the practical value of the clumpy PDR modeling approach.

thumbnail Fig. 1

Simplified structure of a PDR. UV radiation (coming from the left) is progressively attenuated while penetrating a gas cloud. This leads to a chemical stratification and a strong temperature gradient.
thumbnail Fig. 2

General numerical scheme that needs to be solved in KOSMA-τ. Local iterations (red) are performed on every spatial grid point, spatial iterations (green) are performed over all spatial grid points in the model, and global iterations (blue) repeat the spatial loop if necessary until numerical convergence is reached.
Table 1

Typical PDR model input.

Table 2

Typical PDR model output.

3 The KOSMA-τ code

The KOSMA-τ PDR model is based on the plane-parallel PDR code by Sternberg (1988). The spherical KOSMA-τ model has been developed at the University of Cologne in collaboration with Tel Aviv University and the first results have been published by Gierens et al. (1992) and Störzer et al. (1996, 2000). Historically, the model geometry was driven by the finding that observations showed molecular clouds to be clumpy, porous or fractal (e.g., Stutzki et al. 1988). The fractal properties of many clouds can be described by a fractional Brownian motion structure and such a structure can be reproduced by an ensemble of clumps with a well defined power-law mass spectrum and mass-size relation (Stutzki et al. 1998; Zielinsky et al. 2000). Theory and numerical simulations of the ISM also show that turbulent flows naturally lead to a highly fragmented structure (e.g., Gammie et al. 2003). Gong & Ostriker (2011, 2015) showed that dense cores form from supersonic turbulent converging flows and Guszejnov et al. (2018) argue that the observed power-law scaling is the generic result of scale-free structure formation where the different scales are uncorrelated. For a review on clump formation see Ballesteros-Paredes et al. (2020). Consequently, it is necessary to include the effect of clumpiness and fragmented structures in any PDR modeling as far as possible.

The fundamental difference between the spherical model geometry and plane-parallel models is the higher ratio of surface to mass (or volume). The infinite extent of plane-parallel models parallel to the surface is not real but a necessity of the 1-D setup. Two-sided plane-parallel models slightly improve on that aspect but have to introduce an additional model parameter, the cut-off AV. The simplest possible model with a finite configuration is the spherical clump, but as for the two-sided models, it comes at the cost of an additional model parameter, the clump mass. The spherical model is also a strong simplification since molecular clouds are never spherical but it allows for a better description of the clumpy structure of the ISM with a large fraction of internal surfaces. By matching the parameters to the observed clump-mass spectra we can approximate the fractal structure in PDRs, something that is impossible with plane-parallel configurations.

The model geometry naturally explains some observations that are hard to explain otherwise. Bolatto et al. (1999) and Röllig et al. (2006) showed that the observed metallicity influence on the emission ratio of [C II]/12CO and [C I]/12CO can be explained by a spherical (i.e., finite mass) cloud model. Izumi et al. (2021) observed unusually high ratios of [C I] 3P1 → 3P0 to 12CO(1–0) emission at high extinction that cannot be explained by standard plane-parallel PDR models and suggests highly clumpy gas. Schneider et al. (2021) performed a detailed PDR multiline model analysis of a globule in Cygnus X and showed that the observed line emission can only be explained with a two-component PDR: a clumpy internal and a nonclumpy external PDR. Other studies that used KOSMA-τ model predictions to successfully analyze and interpret PDR observations include Mookerjea et al. (2006), Schulz et al. (2007), Kramer et al. (2008), Pineda et al. (2008), Cubick et al. (2008), Sun et al. (2008), Röllig et al. (2011, 2012, 2016), Schneider et al. (2016, 2018), Garcia et al. (2021), Nayak et al. (2021), and Mookerjea et al. (2021). Andree-Labsch et al. (2017) investigated the spatial variations of PDR emission lines across the Orion Bar and showed that the observed spatial profiles can be explained by dense PDR clumps embedded in a thinner inter-clump medium.

The spherical geometry requires more computational power due to the need of angular averaging and the additional mass parameter. The clump ensemble approach requires the computation of large parameter grids with the consequence of a high computational effort. To compensate for this, physical complexity was reduced in certain calculations, for instance the full ro-vib structure of H2 was approximated by 15 virtual vibrational levels only. As computing power becomes more readily available, we are gradually increasing the model complexity again.

The general model iteration is preceded by the model setup and the precomputation of the FUV continuum radiative transfer. Using the multicomponent dust radiative transfer code MCDRT (Szczerba et al. 1997) we compute the spectrally resolved FUV radiation field in the model clump together with the dust temperature (fully resolved for all dust components and dust sizes) (Röllig et al. 2013). In a post-processing step following the global model convergence, we then compute the detailed line & continuum emission, spatially resolved as well as model clump averaged (Gierens et al. 1992).

3.1 Geometry

KOSMA-τ is a 1D model with spherical geometry using the cloud depth z = Rtot − r as spatial coordinate, with the total radius Rtot and the radius r from 0 to Rtot. We usually represent the spatial coordinate by the optical extinction AV along this coordinate8. The total gas density n is described by (1)

where n0 is the total gas density at the surface (r = Rtot) in cm−3. The standard parameters are: α = 1.5, Rcore = 0.2, approximating the structure of Bonnor-Ebert spheres with a density contrast from edge to center of 11.2. We note that similar density contrasts are found in isobaric PDR models (Joblin et al. 2018; Wu et al. 2018; Bron et al. 2018). In general we use the clump mass instead of the clump radius as the model parameter characterizing the size of the clumps.

The spatial loop starts at the surface z = 0 and proceeds into the model cloud until the cloud center is reached. A 1D setup in spherical coordinates automatically implies an isotropic FUV radiation field providing radiation from all directions. This assumption is approximately true for the average ambient FUV field in the Galaxy and for nearby FUV sources if the PDR is embedded in a diffuse medium because of the strong scattering of the interstellar dust in the FUV range. With albedos between 0.3 and 0.4 (Weingartner & Draine 2001b) and the FUV extinction being about three times stronger than at visible wavelengths, a dust column providing an AV ≲ 1 is sufficient to convert a directed FUV field into an isotropic field. Observations and corresponding models seem to support this assumption. Choi et al. (2015) showed that the FUV emission of the H II region around ζ Oph is dominated by dust grain scattering. Similar findings have been reported for the Spica nebula (Choi et al. 2013), the Orion-Eridanus Superbubble (Jo et al. 2011, 2012), and the Taurus-Perseus-Auriga Complex (Lim et al. 2013). A similar scattering effect has been proposed for diffuse Hα emission outside of bright H II regions (Seon & Witt 2012). Therefore the assumption of an isotropic FUV illumination is a good approximation for low to intermediate mass (size) clumps. It becomes questionable for very large model clumps where the radiation would have to be redistributed over parsec scales. However, the very massive clumps are actually a good approximation to the common plane-parallel setup. In Röllig et al. (2007) we could show that for these cases our spherical model agrees well with predictions from plane-parallel PDR codes with a perpendicular FUV irradiation producing comparable PDR structures. This is due to the fact that the emission of very massive structures typically arises from a thin surface layer due to optical thickness effects. This layer itself is dominated by the FUV radiation falling perpendicular onto it. Radiation from other directions and the backside is quickly absorbed in the large columns of material. Even if the approximation of the isotropic illumination breaks down in those cases it has no measurable effect because only the perpendicular contribution affects the large-clump model.

To assure a sufficient numerical resolution of the dissociation front and the transition from atomic carbon to CO we do not assume a fixed spatial model grid. Instead we rely on an adaptive spatial gridding using the Bulirsch-Stoer method (Press et al. 2007, Sect. 16.4) to determine the next step width during each spatial iteration.

3.1.1 Radiative transfer

Solving the radiative transfer (RT) problem is one of the key aspects in any PDR model. The local FUV photon density determines the heating efficiency as well as the strength of the local photo-chemistry. The problem is complicated by the fact that some of these processes are affected primarily by the FUV continuum dust shielding, while others depend on the combination of dust and spectral line shielding (e.g., Weingartner & Draine 2001c) and Heays et al. (2017). Solving the RT fully self-consistently is very time-consuming and many PDR codes compromise on certain aspects of the RT to reduce the computational efforts. The RT in KOSMA-τ is currently divided into three parts.

  • 1)

    The dust RT is done in a preprocessing step using the MCDRT-code (details described in Röllig et al. 2013). MCDRT computes the continuum RT within the spherical clump for a given dust composition based on the respective optical properties of the dust components and an assumed dust size distribution. The result of this computation is the internal FUV field, the emitted continuum spectrum, and the dust temperature distribution in the clump. These are used as input quantities for the KOSMA-τ calculation. MCDRT does not know the gas structure of the clump yet and cannot account for FUV absorption by H2 or CO.

  • 2)

    During the KOSMA-τ iterations the line RT to compute the local photo-processes (dissociation and ionization) are computed during each iteration via a ray tracing scheme that is described below. The line cooling is computed per cooling line based on the local energy level excitation using a spherical escape probability formalism.

  • 3)

    The final clump line emission is computed in a postprocessing step based on the final chemical and physical structure of the clump using the spherical RT code ONION fully accounting for non-LTE (local thermal equilibrium) effects (Gierens et al. 1992). These computations are performed per line ignoring any line overlap or line-pumping between different molecules. The error due to the inconsistent treatment of line cooling and final line emission is small and can be ignored (example [C II] 158µm: median ΔTex/Tex = 3.1% over full parameter grid).

MCDRT uses the Mie theory (Bohren & Huffman 1983) to compute the extinction efficiency Qext and albedo ω for each dust component, assuming the scattering properties of spherical grains. For silicates and graphite we use the dielectric constants from Draine (2003), while for very small grains we followed the approach given by Li & Draine (2001). Scattering by dust particles is currently assumed to be isotropic. More details are given in Röllig et al. (2013).

Computing the continuum radiative transfer only once in a preprocessing step is a reasonable approximating due to the weak thermodynamic coupling between gas and dust. The same approach is not possible for the FUV-line radiative transfer because the chemical structure changes during the iterations. Therefore, KOSMA-τ calculates the line radiative transfer for all relevant line transitions. This includes all absorption processes where spectral line absorption is important, that is the photodissociation of H2 (van Dishoeck 1987; Sternberg & Dalgarno 1989), atomic carbon C, and all isotopologues of CO (van Dishoeck & Black 1988; Visser et al. 2009). Line shielding might be important for other species as well, for example for N2 (Heays et al. 2017), but is not yet implemented into KOSMA-τ.

For all photo-processes, the relevant physical quantity is the local photon flux above critical energy available for the respective process. The flux is determined by the initial unattenuated FUV field and the column of absorbing material which the photons penetrated. Given the spherical geometry and the isotropic radiation field the column density of absorbing material is a function of local position and direction of the incoming photons. The final local flux is the result of the angular average over the full solid angle. To compute the attenuation along all directions we use the method described in Gierens et al. (1992).

Sternberg et al. (2014) showed that the absorption of UV photons in Werner and Lyman H2 lines can dominate the continuum absorption under certain circumstances. Treating these cases fully self-consistently requires a re-computation of the dust continuum radiative transfer and all dust properties after each iteration. During each iteration we compute the pumping of H2 levels and its dissociation as well as the dissociation of carbon monoxide. We use the local dust column densities as function of angle and the effective UV continuum absorption cross section σD (Sternberg & Dalgarno 1989) to compute the continuum absorption in all directions and integrate over all angles (see Sect. 3.1.2). σD depends on the current dust size distribution and is computed by MCDRT.

KOSMA-τ computes the self-shielding of all CO isotopologues using the shielding function provided by Visser et al. (2009) accounting for self-shielding, mutual shielding (line-overlap) and shielding by atomic and molecular hydrogen. The shielding factors fsh(N(X)) are parameterized as function of H2 and CO column densities. The CO shielding varies with the assumed Doppler line width and Visser et al. (2009) report a 26% increase of the CO photodissociation rate for an increase in the Doppler broadening bCO from 0.3 to 3 km s−1. We use CO shielding functions computed for bCO = 3 km s−1, Tex(CO) = 20 K, and a 12C/13C ratio of 69 and neglect variations with bCO. The H2 self-shielding for the radiative pumping by FUV photons is computed line by line based on the shielding prescription by Federman et al. (1979) accounting for the Voigt profile absorption lines. Draine & Bertoldi (1996) provided H2 shielding functions that account for line overlap but do not allow for a line-by-line treatment. KOSMA-τ offers the option to use their shielding prescription for the computation of the H2 photodissociation rate instead.

In this study, we do not resolve the rotational energy levels of H2 but only consider transitions between vibrational levels in the ground state (15 levels) and the electronically excited Lyman B1Σu (24 levels) and Werner C1Πu (10 levels) bands. The vib-level population in the ground state is computed from the balance between collisions, spontaneous decay and UV pumping (Allison & Dalgarno 1969; Dalgarno & Stephens 1970; Stephens & Dalgarno 1972, 1973; Abgrall et al. 1992). The photodissociation rate is computed from decay of Lyman and Werner states, excited by UV pumping, to an unbound state (Sternberg 1988).

Continuum processes that depend on the local FUV intensity, such as the photoelectric heating, are computed based on the intensity results of MCDRT and neglect the effect of line absorption on the local mean FUV intensity. We do not account for line-overlap and treat the radiative transfer line by line.

thumbnail Fig. 3

Structure of the shell ray tracing setup. The shells are arranged equidistantly for clarity. Real spatial grids are not equidistant. The angle index t is the second index given at the blue labeled points.

3.1.2 Ray tracing scheme

A ray tracing approach is applied to determine the local photon flux inside the model cloud in order to compute photo-processes such as photo-dissociation and photo-ionization. it is also applied to compute the final emission characteristic of the model cloud. KOSMA-τ is a 1-D code therefore the angular gridding has to be derived from the existing spatial grid.

The model consists of a series of concentric, nonequidistant shells with radius ri. The local conditions are constant across a particular shell, for instance all positions indicated by the blue points in Fig. 3 share the same physical and chemical conditions. The spatial computation iterates over all shells, that is computes all red points in the figure. The cloud is intersected by a series of parallel rays qi, one across the center and the others tangential to the shells. The crossing points between the shell boundaries and the intersecting rays (black in Fig. 3) and the angles of incidents at these points constitute the spatial and angular grid that is used to perform the radiative transfer computations. An angle index t = 0 corresponds to the tangential point between ray and shell. The number of spatial grid points and therefore shells and rays in real calculations is in the order of a few hundred depending on the cloud parameters.

The grid points are labeled with their shell index s and their angular index − s ≤ t ≤ + s as shown in Fig. 3. Intensities and relevant quantities such as column densities are computed point by point along a ray and for all rays q. We use the radiative transfer equation: (2)

where s, t and s′, t′ are adjacent points on the same ray and κ and S are the opacity and source function in the shell between the points. l is the distance between the two points. It is convenient to reorder the grid points s, t → q, t where q = s − |t| + 1 is the corresponding ray index. Along each ray q the angle index runs from to . Using the transfer equation we can recursively compute the intensities along a given ray q (3)

with and s, s − 1 the corresponding shell indices. Ibg is the background intensity. Accordingly, the intensity at left surface points in Fig. 3 is Iq,t = Ibg = Bv(Tbg) with Planck’s radiation law Bv at the background temperature Tbg = 2.73 K, while for the rest of the points we have a step-wise integration of the radiative transfer equation. The detailed energy level population of a species is obtained from the corresponding set of non-LTE rate equations which are solved with a generalized Newton–Raphson technique. We derive the line source function S1 from the corresponding excitation temperatures. S = S1 + Sc where Sc is the continuum contribution from the local dust temperature also contributing to the pumping of the quantum levels. Similarly, we compute the total opacity from their line and continuum contributions κ = κ1 + κc. Contributions from external dust continuum are not considered. We use 16 frequency points to resolve the (Gaussian) line profile and apply a Gauss–Hermite quadrature for the frequency integral. Collision rates, A-values and line frequencies are taken from the Leiden Atomic and Molecular Database (Schöier et al. 2005; van der Tak et al. 2020).

The final result is the emergent intensity for each ray and the intensity inside the cloud for all shell-ray intersections. For any given cloud depth, that is shell, this allows us to compute the angular average over intensities coming from all directions: (4) (5)

where I(rs, θ) are the intensities across the shell s (e.g., at all blue points in Fig. 3) and θs,t are the angles at the points s, t, θs,t = arccos(rs−|t|/rs]).

A similar ray-tracing scheme is also used to compute for instance the H2 and CO self-shielding of the UV absorption lines and the corresponding H2 pumping along each ray (Störzer et al. 1996). In our model geometry column densities are functions of the angle θs,t. We compute N(X)s,t the column density of species X at point s, t using a simple trapezoidal rule. Equivalent to Eq. (2), we compute the attenuation of all photo-reactions in the chemistry by using the angle-dependent columns in the attenuation factor exp(−γAV) from Eq. (D.2) where we expressed column density by AV in order to compute the reaction-specific local photo-rate. To compute the clump averaged emission we compute: (6)

For more details on the radiative transfer scheme see Gierens et al. (1992). We note that Eqs. (2) and (4) have to be computed for all relevant radiative transitions. For species with many transitions this introduces a significant numerical complexity. In particular the solution to the full rotational-vibrational structure of the H2 molecule, including the UV transitions, involves almost 30 000 radiative transitions: Thus, the described scheme has to be solved for 30 000 transitions along several hundred rays at each spatial step. In the past this was a prohibitive effort, but we recently extended the code to include the full H2 problem. This will be described in detail in a subsequent paper (Röllig & Ossenkopf-Okada, in prep.).

3.2 Gas-phase chemistry

The chemistry in KOSMA-τ is modular. Chemical species can easily by added or removed from the chemistry and the code selects all participating reactions from the chemical database in use. The standard database is based on the 2012 edition of the UMIST Database for Astrochemistry (UDfA12; McElroy et al. 2013)9 including a number of modification and updates: isotopologues including 13C and 18O have been added (Röllig & Ossenkopf 2013); new branching ratios from Chabot et al. (2013) are used; addition of l-type isomers (UMIST only contains c-type); updated fractionation reaction rates (Mladenovic & Roueff 2014); refitted low-temperature rates (Röllig 2011).

From the list of reactions KOSMA-τ constructs a set of ordinary differential equations (ODE): (7)

The first two terms sum over all two-body processes forming species i and all photo-processes and/or cosmic-ray processes leading to the formation of species i, respectively. The final two terms are the equivalent sums over all processes leading to the destruction of species i by two-body reactions (index m) and photo-processes (index m′). The prefactor k is called the rate coefficient and tabulated in parameterized form in chemical databases. We use lower-case k and upper-case K to distinguish gas-phase and surface reaction rate coefficients, respectively. We list the implemented gas-phase reaction types in the Appendix D.

The set of N rate equations (Eq. (7)), one for each of the N chemical species is complemented by a set of conservation equations per chemical element M: (8)

where is the number of atoms of element M in species i. Analogously we also find a charge conservation equation: (9)

where is the charge of species i in units of the elemental charge where can be negative in contrast to . Charge exchange between grains and gas is currently not considered.

3.3 Surface chemistry

A large number of observations of star formation regions show that observed gas-phase abundances of species other than H2 such as NH3, N2H+ and CH3OH cannot be explained with pure gas phase chemistry (e.g., Geppert et al. 2005; Garrod & Herbst 2006; Herbst & van Dishoeck 2009; Bottinelli et al. 2010; Oberg et al. 2011; Boogert et al. 2015). The freeze-out of gas-phase species in the dark region of molecular clouds also removes important tracer species from the observable content of the ISM and modifies the IR/FIR line emission of the clouds.

Historically, the only grain surface chemistry included in PDR codes was the formation of H2 followed by immediate desorption to the gas-phase. Based on measurements in the diffuse medium a mean H2 formation rate of 3 × 10−17 cm3 s−1 was found (Jura 1974). This formation rate or comparable temperature dependent parametrizations (e.g., Tielens & Hollenbach 1985; Sternberg & Dalgarno 1995) was implemented in astrochemical models as gas phase reaction simulating the surface process. Later, other prescriptions have been suggested to simulate the formation of H2 more precisely, including considering Langmuir–Hinshelwood and Eley–Rideal processes (e.g., Le Bourlot et al. 2012), physisorption and chemisorption of H2 (Cazaux & Tielens 2002, 2004), as well as stochastic treatment of dust temperatures and grain populations (e.g., Barzel & Biham 2007; Bron et al. 2014). KOSMA-τ implements the H2 formation formalism described by Cazaux & Tielens (2002) and Cazaux & Tielens (2004). Details are described in Röllig et al. (2013).

It is still unclear whether surface chemistry takes place on very small particles such as polycyclic aromatic hydrocarbons (PAHs) (Boschman et al. 2015; Andrews et al. 2016; Foley et al. 2018). To estimate the chemically active grain surface, we used the H2 formation rate observed in the diffuse medium. Cazaux et al. (2016) derived the grain surface per H-atom that is available to H2 formation as 4 × 10−21 cm2/H-atom. We used this value to compute the minimum grain radius for surface chemistry, integrating over the dust size distributions from Weingartner & Draine (2001b) that apply to diffuse gas conditions. Assuming that the limiting factor is the thermal stability of the dust grains, carbonaceous dust can be approximately 20% larger than corresponding silicate grains (Li & Draine 2001). From these two conditions we find minimum radii of 24 Å and 29 Å for silicate and carbon dust by integrating over dust size distributions suitable for diffuse gas. This is comparable to a threshold of 20 Å given by Hollenbach et al. (2009). Using this as lower integration limit we find total surface areas Asil = 2.6 × 10−21 cm2/H-atom and Acarb = 1.4 × 10−21 cm2/H-atom for silicate and carbonaceous grains, respectively.

Other surface reactions have been introduced to numerical PDR models with a primary focus on formation routes of chemical species observed in PDRs, such as H2O and H2CO (for example Hollenbach et al. 2009; Guzmán et al. 2011; Le Bourlot et al. 2012; Esplugues et al. 2016; Putaud et al. 2019). The updated chemistry in KOSMA-τ now includes all relevant surface processes in a quasi-three-phase model (gas + surface + inert ice bulk). For the surface chemistry we follow the rate equation approach as described in Hasegawa et al. (1992) and Hasegawa & Herbst (1993) including the competition between different processes as described by Chang et al. (2007) and Garrod & Pauly (2011). In our quasi-three-phase model we limit the mobility of surface species to the top surface layer following (Cuppen et al. 2017) and desorption to the top two ice mantle monolayers (MLs) (see Eq. (D.17)). For surface reactions with an activation energy barrier we account for competing processes such as diffusion and desorption following (Garrod & Pauly 2011). Including surface chemistry increases the number of chemical species in the network because gas-phase species and surface species must be considered as two different chemical species for all mathematical aspects involved. In the following, the symbol ns denotes chemical species on grain surfaces. The rate equations of the gas-phase chemistry Eq. (7) is extended and modified. We get an additional set of equations governing the surface processes: (10)

where, K1 and Kjk denote the surface reaction rate coefficients for one-body and two-body reactions, respectively. We also have to extend the gas-phase chemistry Eq. (7) by the additional accretion and desorption terms: (11)

where (…) corresponds to the right hand side of Eq. (7). Here Kdes describes all thermal and non-thermal desorption processes, that is a conversion of a surface species to a gas-phase species while kacc is the rate coefficient for accretion (freeze-out) which converts gas-phase species into their surface equivalent.

We assume a surface density of binding sites on the grains nsite = 1.5 × 1015 cm−2 (Tielens & Allamandola 1987). For tunneling through activation energy barriers we assume an energy barrier width of a = 2 Å (Garrod & Pauly 2011). For the desorption (or binding) energies Ed, we use values provided by KIDA (Wakelam et al. 2017b) with the modifications given in Table 3. Typically, Ed is given for species bound to a H2O ice surface as the most abundant ice component. Desorption energies for other surfaces may differ significantly (Esplugues et al. 2016) and we include different desorption energies for species bound directly to a carbonaceous surface if they are known (see Table 4). For a more thorough discussion on desorption energy choice and the resulting model sensitivity we refer the reader for example to (Penteado et al. 2017) and (Kamp et al. 2017). We use the notation J(X) to distinguish surface species from their gas-phase counterpart X.

In Appendix D.2 we summarize all surface reactions that are currently implemented in KOSMA-τ. Here, we only report the new chemical desorption framework that we included in the surface chemistry.

H2 ice. A significant population of H2 ice in the ISM has been proposed already almost 30 yr ago based on early astro-chemical models (Sandford & Allamandola 1993). However, no clear observational evidence could be detected to date. In the context of numerical models, freeze-out of gas-phase H2 is potentially problematic because it can lead to a catastrophic freeze-out of all molecular gas onto dust grains. Several physical processes have been discussed to prevent nonphysical H2 ice populations and we include encounter desorption for J(H) and J(H2) as suggested by (Hincelin et al. 2015). The underlying idea is based on the significantly lower desorption energy for hydrogen on an H2 substrate compared to other substrates (e.g., Vidali et al. 1991; Cuppen & Herbst 2007; Das et al. 2021). For J(H) or J(H2) located next to a J(H2) molecule we decrease Ed from its canonical value to Ed = 23 K (Cuppen & Herbst 2007) greatly enhancing the possibility for a desorption of J(H) and J(H2). As a consequence, the further accretion of H2 ice is strongly suppressed after the build-up of the first ML.

Dust model. The dust properties are computed by MCDRT assuming a dust composition and dust size distribution, for example as given by (Mathis et al. 1977) or (Weingartner & Draine 2001b). The dust radiative transfer and dust properties are computed per dust size bin together with corresponding mean values, for instance for dust temperature, dust surface area, and dust density. Details are described in (Röllig et al. 2013). Presently, KOSMA-τ uses a single surface weighted average grain temperature Td for all types of grains independent of size and grain composition10. For the KOSMA-τ computations presented in this paper we assume the dust model #7 from (Weingartner & Draine 2001b) consisting of four dust types: carbonaceous grains, silicates, PAHs, and ionized PAHs. In this model the average radii for silicates and carbon grains (including VSGs and PAHs) are 〈〈a〉〉sil = 4.5 µm and 〈〈a〉〉carb = 4.0 µm, when 〈〈 〉〉 indicates a surface weighted average. For the photoelectric heating and grain recombination cooling we use the prescription from (Weingartner & Draine 2001a) corresponding to the assumed dust size distribution. This implicitly includes the computation of the grain charge distribution. We neglect any feedback of line UV absorption on the dust temperature.

Cosmic-ray induced desorption. Cosmic rays (CR) hitting dust grains deposit an energy of ECR ≈ 0.4 MeV into dust particles heating an average grain of 0.1 m radius (Leger et al. 1985) to approximately TCR,max = 70 K (Hasegawa & Herbst 1993), sufficiently hot to induce significant thermal desorption. Cooling of the grain takes place via radiative cooling and via evaporation cooling due to the sublimation of bound ice particles. The cooling time scale of a grain in s is given by (12)

where the cooling rate due to sublimation is given by: (13)

where the sum includes all ice species, 〈ndust〉 is the average number density of the dust particles, and fdes gives the fraction of surface species that are candidates for desorption so that Idesns(i)/〈ndust〉 gives the number of desorbable molecules of species i per dust grain. v(i) ≈ 3 × 1012 s−1 is the desorption attempt frequency given by the surface vibrations (Hasegawa & Herbst 1993). Radiative cooling is computed by with the Stefan-Boltzmann constant σ and qabs = 0.13 K−2 cm−1 for silicate grains. Cooling is always fast compared to the frequency of CR hits for any dust grain.

Assuming that most desorption occurs at about TCR,max = 70 K we can approximate the CR-induced desorption rate coefficient as: (14)

where is the fraction of time spent by grains at the temperature , given by the cooling time relative to the frequency of CR hits. Kevap,i(TCR,max) is the thermal desorption rate coefficient of species i at TCR,max. Assuming s (Hasegawa &Herbst 1993) find fcr = 3.16 × 10−19 assuming ζCR = 5 × 10−17 s−1. To use this formalism we can distribute the total available dust surface in our model across virtual average grains of radius 〈aD〉 = 0.1 µm and find an average grain density of 〈ndust〉 = 3.17 × 10−12 per H-atom. Using the same description of 0.1 µm grains (Sipilä et al. 2021) published an updated estimate based on time-dependent computations of the heating and cooling due to CRs. Their heating depends on ζCR and the local AV and can be approximated in the range ζcr = [1.3 × 10−17, 10−16] s, and AV = [0, 100] mag by: (15)

where z = log10ζCR. Energy deposition by CRs directly into the ice mantle (Wakelam et al. 2021) is not considered but will be implemented in a future code update.

Chemical desorption. (Minissale et al. 2016) and (Cazaux et al. 2016) presented an analytical expression for an additional desorption mechanism using the released binding energy from exothermic reactions in the ice-phase to desorb species i to the gas-phase first presented by (Dulieu et al. 2013). They fitted a set of experimentally measured desorption rates through a probability that the reaction product has an energy higher than the desorption energy ED,i: (16)

where ΔHR is the exothermicity of the formation of the product(s) and df specifies the degrees of freedom that the energy is distributed among. A good fit was obtained when using a description that is equivalent to an elastic collision with a grain mass of M = 120 amu where the energy fraction obtained by the product is ϵi = (M − mi)2/(M + mi)2. This approach, however, is not applicable if multiple reaction products are generated that can desorb. In this case the released energy ΔHR must be distributed over all reaction products. A modification of Eq. (16) is needed which accounts for the energy distribution between multiple reaction products: (17)

where ηi specifies the fraction of ΔHR that is available to species i. Unfortunately, no laboratory data seem to be available for the desorption rate in case of more than one reaction product so that ηi cannot be determined from a fit to the experimental data. (Wakelam et al. 2021) used an equal split between all atoms within the reaction products, assigning equal energy fractions to all atoms within the molecules. We propose a different approach based on the nature of the chemical reaction that produces the exothermicity. We assume that the fraction ηi of the binding energy that is transferred to the individual products is proportional to the number of open shell electrons ei that contribute to the formation of each product (ignoring closed shells). Taking for example the reaction J(O) + J(HCO) → CO2 + H this gives the following ei: eH = 1, eC = 2, eO = 4. (18)

where the sum is over all atoms involved in the reaction. For each product this quantifies the probability for desorption and for remaining on the surface . In the general case of two reaction products we can compute the branching ratios for all four possible reaction branches: (19)

where R1,2 and P1,2 stands for two reactants and products respectively, and the subscript s denotes species on the surface or in the ice mantle. Take for example the reaction J(O) + J(HCO) → CO2 + H: we find ηH = 1/11, and . Following Eq. (19) gives the branching ratios (values in parenthesis would result from the application of Eq. (16) to all products, this means a value of ηi = 1 for all products.): (20)

The final reaction rate coefficient for the chemical desorption of i resulting from the reaction between species j and k is: (20)

with Kjk from Eq. (D.8). For the computation of the BR we always use the ED,i for a H2O substrate. Tests have shown that the chemistry is not very sensitive to variations in the assumed desorption energies. The list of all chemical desorption reactions and their branching ratios is given in Appendix D.2, Table D.4.

Table 3

Modified desorption energies used in the surface chemistry.

Table 4

Desorption energies for binding on bare grain surfaces.

3.4 Thermal balance

Thermal balance is the local balance between all cooling Λk and heating Γh processes: (21)

where the efficiency of any individual process depends on the local conditions , such as the chemical vector n = (n1, …, nN), the column density vector N(θ) which is a function of direction θ, gas and dust temperature Tg/d, FUV intensity χFUV (in units of the Draine field Draine 1978), the cosmic ray ionization rate ζCR, and possibly others. The superscript i is the global iteration counter, the subscript x enumerates the positions.

Equation (21) is coupled to the problems of chemical balance, local excitation and radiative transfer (see Fig. 2). The importance of various heating and cooling mechanisms as a function of position in the PDR has been widely discussed in the past (e.g., Tielens & Hollenbach 1985; Hollenbach et al. 1991; Bakes & Tielens 1994; Sternberg & Dalgarno 1995; Röllig et al. 2006; Woitke et al. 2009). At low AV, photo-electric heating and H2 vibrational de-excitation are the dominant heating, at high AV, cosmic ray heating (CR) and gas-grain collisions are the most important heating processes. Cooling at the surface of the PDR is provided by [C II] and [O I] fine-structure line cooling, by CO rotational line cooling at high AV, and by dust (Ossenkopf et al. 2015). The photo-electric heating depends on the grain properties. For an assumed dust distribution given by (Weingartner & Draine 2001b) we implemented the formalism given by (Weingartner & Draine 2001c, their Eqs. (44), (45) and Tables 2 and 3). Details of our implementation are described in (Röllig et al. 2013). If the user selects a MRN dust distribution (Mathis et al. 1977) we use the photoelectric heating given by (Bakes & Tielens 1994, their Eqs. (41), (43), (44)). Table 5 lists all heating and cooling processes implemented in KOSMA-τ and Details are discussed in Appendix E.

3.5 Local convergence

Solving the local problem is equivalent to finding a local gas temperature such that the local chemistry and all the local thermodynamics are in equilibrium. To solve the chemistry as well as the heating and cooling we need to compute the local radiation field, the dust temperature and the local photon escape probability (see Local iteration loop in Fig. 2). In plane-parallel PDR models, this can usually be achieved by assuming exponential attenuation along the line of sight. In spherical PDRs, this is complicated by the fact that we assume an isotropic FUV irradiation and that absorbing columns of gas and dust depend on the direction. To compute the local intensity in a spherical cloud we then need to average over all directions: (22)

Along any given line-of-sight the attenuation of radiation at a given frequency can be described in terms of exp(−γvAV). γv relates the attenuation at a frequency v to the visual extinction. Using the relation Av = 1.086τv it is common practice to use the visual extinction AV as “optical” coordinate. A quantity such as the local optical depth τv = σvN (Sternberg & Dalgarno 1989), with frequency v dependent cross section σv and the related column density of the absorber N, has to be considered for all directions θ: (23)

The flow of the local iteration is summarized in Algorithm: Local iteration 1.

Table 5

Heating and cooling processes implemented in KOSMA-→.

Algorithm: Local iteration 1

Initial values:
z → znew, Tg → Tg,old
Compute Td(z)
angular average of radiative quantities
Compute →fuv(z)
reaction specific dust attenuation
self shielding CO
H2 UV shielding and pumping
Compute Tg(z)
while Etot → 0 do
  T → Tg
  solve level population: H2,CO,O,C, …
  solve chemistry
  compute local heating & cooling →(T), →(T)
   Etot → →(T) → →(T)
  choose root finding algorithm
  
end while
Tg(z) → T

3.6 Spatial loop

The physical and chemical conditions in a PDR are subject to nonlinear effects such as the photo-dissociation of H2 and CO and their respective shielding. The photo-dissociation rates show a sudden drop once the absorption lines become optically thick. As a consequence, we find a steep density gradient of H2 and CO in these transition regions covering many orders of magnitude in volume density. To avoid numerical problems, it is important to use a sufficiently dense spatial grid across these transition regions and the method described below performs well across the full parameter range.

During each spatial iteration KOSMA-τ performs an adaptive spatial gridding using the Bulirsch-Stoer method (Press et al. 2007, Sect. 16.4) to determine the next step width such that the numerical solution of the set of coupled first-order differential equations for the functions , having the general form (24)

is sufficiently accurate. In the context of solving the PDR structure, we identify yi(z) = Ni(z) as the perpendicular column density of species i from the surface to depth z in the cloud so that dyi(z)/dz = ni(z). This extends to the column densities of energy level populations of species of interest, in particular CO, C+, C, O, and H2. From a given position and column density vector, the Bulirsch-Stoer stepper routine returns the next step-width Δh together with the respective values for ni(z + Δh) and Ni(z + Δh) as well as an estimate of which step-width to attempt next. This successively propagates the solution into the cloud until the cloud center is reached.

3.6.1 Global convergence

KOSMA-τ uses the final column density vector, that is the column density of all species measured from the edge of the cloud to the cloud center, N to test for convergence. Global convergence is met if (25)

The standard value applied is ϵiter = 0.01. The minimum number of iterations is 2. We also specify a maximum number of iterations (default of 60) before stopping further iterations (Most models converge within less than 30 iterations). If more than 40 iterations are needed we incrementally relax the 0.01 criterion in steps of 0.01 to 0.1. We try another ten iterations with ϵiter = 0.1. If no convergence is found in 60 iterations the result of the last iteration is stored because for certain applications the results might still be sufficient. It is for example in some cases numerically difficult to find a stable solution for massive dense clouds. There, the low J transitions of CO are the only cooling options for the gas in the central parts because of the low gas temperatures. The respective optical depths however are so high that Eq. (21) becomes strongly nonlinear and small changes in Tg lead to large differences in the CO population numbers that can produce oscillating solutions for the population column density vector. However, this happens for such high optical depths that the effect is not observable and has no impact on the total CO emission of the model cloud.

In general, we find that for a given number of chemical species the code converges faster for low densities, low masses (i.e., low optical depths) and not too high FUV fields. The total computation time for one cloud model, that is one specific set of model parameters range between few minutes and 10–20 h on an Intel Core i7 CPU.

3.6.2 Chemical time scales

To assess the limits of the model assumption of stationary chemistry we have to compare the involved chemical time scales with other relevant time scales, such as the dynamical time scale of the clump or its total lifetime. The model looses predictive power if those times are comparable or larger than the time the chemistry needs to reach steady-state. We derive chemical times from the total formation or destruction rate of a species divided by its total number density. Taking the inverse is the relevant formation/destruction time per particle and it will depend on the local conditions. Generally, the chemistry is very fast close to the surface with significant local photon densities and high temperatures and slows down with AV. It typically peaks once AV > 1–3. The chemistry is slowest for low density and FUV illumination (e.g., n = 103 cm−3 and χ = 1) with gas-phase time scales ≲ 106 yr and fastest for high density and FUV illumination (e.g., n = 106 cm−3 and χ = 106) with gas-phase time scales ≲104 yr. The surface chemistry is considerably slower once AV > 10 with values of 108 yr and longer. Our model showed the time scales for J(CO) and J(CH4) to be about at least 1–2 orders of magnitude faster compared to J(H2O) and J(CO2).

Typical dynamic time scales are on the order of 105–106 yr which is larger than the gas-phase time scales under typical PDR conditions and of the same order for low density and FUV conditions. However, considering that PDRs typically evolve from cold molecular gas it is reasonable to assume a total life time of up to ~107 yr which is sufficient to establish chemical steady-state independent of AV. For the surface chemistry even this time may be too short to reach a stationary state for example for water ice deep inside the clump (AV ≳ 10). We note that our models sometimes show a significant deviation between total formation and destruction rates of some ice species, due to our small network of surface reactions. We conclude that the assumption of steady-state is justified for the gas-phase chemistry under most PDR conditions and may become problematic in low n0 and low χ environments. The ice chemistry can reasonably well be approximated with a stationary chemistry for AV ≲ 5–10 and even deeper in the cloud for some ice species, such as J(CO; see also Hollenbach et al. 2009). On the other hand, time-dependent computations by (Esplugues et al. 2019) show that J(H2O) abundance at strong FUV illumination and n = 105 cm−3 reaches a plateau at AV > 7 after 106–107 yr. We note that previous studies showed the existence of multiple chemical solutions (Le Bourlot et al. 1993b, 1995; Lee et al. 1998; Charnley & Markwick 2003; Boger & Sternberg 2006; Dufour & Charnley 2019) as well as the existence of sustained chemical oscillations (Roueff & Le Bourlot 2020) in gas-phase models of the ISM.

4 Chemical solution

4.1 Solving the chemical equations

KOSMA-τ computes an equilibrium solution to Eqs. (7) and 11, that is the solution to the system (26)

using the Newton–Raphson method with the Jacobi matrix . This is a locally convergent method that is guaranteed to converge for a starting point sufficiently close to the root. The elements of are: (27)

For the unknown densities ni of species built up from elements we find an over-determined set of equations. The default method in KOSMA-τ is to replace equations from the set of chemical rate Eqs. (7) and (11) with the corresponding elemental and charge conservation Eqs. (8) and (9). This can always be justified because the conservation equations should be strictly fulfilled while all rate equations are in any way only approximately known. Moreover, it improves the numerical stability of the iterative solution.

The Newton–Raphson method approaches the solution by using an (old) approximate solution (or initial guess) to compute an (new) improved approximate solution in the next step: (28)

This is a linear(ized) system of the form .To solve Eq. (28), we need to invert in order to determine the new density vector nnew. In standard Newton–Raphson algorithms, the new solution-step is computed as n(i)new = n(i)old + δni and used as updated n(i)old value in the subsequent iteration step. The steps can be expressed as relative changes: n = δni/n(i)old. KOSMA-τ solves Eq. (28) by LU decomposition but has alternative solvers implemented.

4.2 Newtonian stepping strategies

The Newton–Raphson algorithm is guaranteed to provide fast convergence if one is close to solution, but depending on the starting values and the overall topology of the Newtonian vector field δn, that is the solution to the system , the Newton–Raphson algorithm may also diverge or become trapped in a local minimum leading to infinite oscillations. Adaptive stepping strategies may prevent this.

We can express the step in the Newton–Raphson algorithm, n(i)new = n(i)old + δni, as a relative change, so that n(i)new = n(i)old · fstep with (29)

For large steps, this approach can be problematic because it does not prohibit negative solutions, that is negative densities. In general, Newton stepping through the negative domain is not problematic as long as we can guarantee that the steps will finally provide positive densities. However, this is problematic because the algorithm may converge on a local minimum involving negative ni. One way to avoid this is to prevent the Newton steps from producing sign switches in n.

The tanh() allows for a convenient construction of such a general stepping function whose symmetry and limits can be controlled with few numerical parameters: (30)

ω+ > 1 and ω− > 1. For small η Eq. (30) approaches Eq. (29) and levels off for large arguments preventing too large steps. We implement an adaptive choice of ω+, ω−, and η, with details described in the Appendix C.1.

5 Energy solution

The numerical results presented in this section have been computed with WL-PDR, a simplified, plane-parallel PDR model written in Mathematica (Wolfram Research Inc. 2020). WL-PDR is designed to act as numerical testing environment of PDR modeling aspects. A brief explanation of the code is given in Appendix G.

5.1 Numerical iteration scheme

The key part in Fig. 2 is the local solution module. A local solution in this context is the local balance between all cooling Ak and heating Γh processes (Eq. (21)), where the efficiency of any individual process depends on the local conditions . Equation (21) is coupled to the problems of chemical balance, local excitation and radiative transfer (see Fig. 2). Figure 4 shows an idealized form of Eq. (21). At any given position in the PDR, every point on the curve Etot(T) corresponds to a different chemical vector and a different temperature. At the solution Etot(Tsol) = Esol = 0 heating and cooling are balanced. If Etot > 0 (gray rectangle left of Tsol) we have heating excess Γ > Λ. In the opposite case, Etot < 0 (gray rectangle right of Tsol) we have a cooling excess Γ < Λ11. Physics ensures that limT→0 E(T) > 0 because all cooling process become inefficient at very low temperatures and heating is less dependent on T. Conversely, it is also ensured, that . Here Tmax corresponds to the upper temperature that is considered reasonable in a PDR and which depends on the parameter range that the model needs to cover. Typical values are Tmax ≈ 10 000–20 000 K. Even though many coolants are inefficient at very high temperatures – CO and H2 for instance are chemically destroyed – other cooling processes become strong. Typical candidates are Lyman α, O I 6300 Å cooling (Bakes & Tielens 1994; Spitzer 1978) and gas-grain collisional cooling. These processes ensure that E(Thigh) < 0. If Etot is a continuous function the intermediate value theorem states that at least one root must lie in the interval [0, Tmax] and a physical solution of the problem exists.

The solution Tsol is stable if ∂Etot/∂T < 0. Imagine the gas at Tsol is slightly perturbed toward higher temperatures, for example toward Eb. In this regime, the cooling exceeds the heating and drives the temperature back to Tsol. The same happens for perturbation toward Ea. The stronger heating drives the temperature back into the equilibrium.

thumbnail Fig. 4

Idealized form of Eq. (21) with a single root. Numerically, a temperature solution Tsol is determined by finding the root of the energy balance equation Esol = Etot(Tsol) = 0.

5.2 Multiple temperature solutions in PDR models

The energy brackets E(Tlow) > 0 and E(Thigh) < 0 guarantee that at least one root exits, but it does not prohibit the existence of multiple roots. For numerical root finding algorithms this poses a serious problem, because they usually find only one root. Which one depends on the choice of the algorithm, the initial starting temperature and the temperature range where the root is searched.

The form of Etot(T) can change when moving through the PDR. In Fig. 5, we sketch such a scenario. At low values of AV we find a high temperature solution Thigh that is the root of the upper, gray curve in Fig. 5. We denote the high and low temperature solutions as hot atomic medium (HAM) and warm molecular medium (WMM), respectively, as we see in Fig. 6b that the high temperature solution is usually also associated with significantly less H2 than the low temperature solution. At higher values of AV the local conditions allow for an additional cooling process to become efficient in an intermediate temperature range producing a local minimum of the thick line in the figure. We compute the PDR structure from the outside to the inside. After determining Thigh at the low AV value we move to a deeper cloud position and search the temperature root there. To speed up root search the new temperature root is usually searched in the neighborhood of the previous solution. Assuming that the spatial gridding is fine enough, we expect the temperature to change slowly between subsequent numerical steps. However, this strategy is doomed to always find the solution T3 and will not reach T1 (or T2).

If we continue the computation steps at even higher AV we have to find a significantly lower root at Tlow. In this three-step picture, we have a large temperature jump from T3 to Tlow. In real PDR model computations this is a sudden transition from temperatures of a few thousand K to few hundred K across a small spatial range, especially in models with high gas density and high FUV illumination. Depending on the search range of the applied temperature finding algorithm this sudden temperature jump can lead to severe numerical instabilities and prohibit global convergence.

thumbnail Fig. 5

Idealized form of Eq. (21) with multiple roots. Filled points indicate a stable solution, open points are unstable against perturbations. The energy balance for two different positions in the PDR are schematically drawn in gray together with their equilibrium temperatures Thigh and Tlow.

5.3 Physics of multiple PDR temperatures

The described multiple temperature solutions were first reported by (Burton et al. 1990) for n = 106 cm−3 and G0 = 103. They described two temperature solutions with different heating and cooling balances. We find similar effects for even higher FUV fields. It is important to include this behavior in the PDR code to evaluate knock-on effects for the chemical structure and the observable line intensities. Fig. 6 summarizes the model results (computed with WL-PDR) for the benchmark case V4 (n = 105.5 cm−3, χ = 105). In panel (a), we show the temperature and chemical abundance profiles. The dashed curves corresponds to the HAM solution, the solid curves to the WMM. In panel (b), we plot the total energy balance vs. the gas temperature. At low AV values we find a strong contribution by vibrational H2 de-excitation heating (pump heating) and H2 line cooling. They inherit the nonmonotonous temperature dependence from the H2 abundance with a weak maximum around 5000 K. H2 formation heating is also strong and even remains so at much higher temperatures while the pump heating and line cooling drops off beyond 7000 and 8000 K.

Moving even deeper into the PDR, H2 cooling becomes stronger and starts to suppress the roots at temperatures higher than the WMM solution Tsol = 910 K. In the case AV = 0.75 (only shown in panel b) as red curve) we find that Etot(T > 910 K) < 0. Looking at Fig. 6b shows that at 6000 K the energy is almost balanced with Etot = −5.5 × 10−20 erg s−1 cm−3.12

In (Röllig et al. 2006), we discussed how the dominant heating and cooling at the surface of PDRs changes with gas density n and with FUV intensity χ. The energetic behavior described above is associated with the cooling (and heating) capabilities of molecular hydrogen H2. However, we that this is not due to increased cooling in the high density case. The important point is that the cooling needs to be balanced by a sufficiently strong heating in order to produce a temperature solution (Etot = 0). In case of low FUV or low density, we do not have an efficient heating process in this regime and accordingly we do not find multiple solutions. Adding additional or stronger heating to the model can possibly alter this behavior.

Such phase transitions result from physical processes that become effectively inefficient once they cross an energetic threshold. The main difficulty is that even though such transitions are physically possible, their prescription in numerical models suffers from inherent uncertainties: (1) numerical approximations of complex physics, (2) unknown or inaccurate material constants (e.g., collision rates, spectroscopic constants, binding energies), (3) unavoidable model simplification (e.g., geometry, numerical resolution), (4) numerical inaccuracies. Any of these can either lead to a phase transition or prevent it.

In the model V4 at AV = 0.75 we noted in Fig. 6 that Etot(T) is smaller but close to zero between 5000 and 7000 K. This is because the total cooling is only slightly stronger in that temperature range than the total heating. In Fig. 7, we show how changing a single heating process can change the situation. The orange line shows Etot in case of a 20% enhanced H2 formation heating. This would result in two additional temperature solutions, that is a possible temperature solution about seven times higher than in the case of standard H2 formation heating.

Predicting the effect of the two different temperature solutions on observable line intensities is not straight forward because excitation conditions and column density effects are mixed. We found up to 20% higher CO emissivities for the upper lines and up to 70–80% higher CH+ line intensities for higher transitions when selecting the HAM solution instead of the WMM solution.

The final choice between WMM and HAM depends on the modeled physical scenario as already described by (Burton et al. 1990). An initially hot and fully atomic (or ionized) medium that becomes subject to significant lower FUV illumination will settle down at lower temperatures coming from a much higher temperature. This clearly prefers the HAM solution. In the opposite case, where the gas starts from cold and molecular conditions at the onset of FUV illumination, the WMM solution is favored because the gas is gradually transitioning from lower to higher temperatures. A scenario for the second case is the onset of massive star formation in the vicinity of the model cloud while the first scenario could occur in case of some shielding due to cloud motions or sudden changes in the illuminating source due to supernova events.

thumbnail Fig. 6

Two-temperature solution for the V4 benchmark model (n = 1055 cm−3, χ = 105). Left: Tg/d (top) and ni (bottom) profile. Solid and dashed lines correspond to the low and high temperature solutions (WMM and HAM. i.e., T1 and T3 from Fig. 5). Right: Etot(T) at selected positions in the PDR (marked by black points in Fig. 6a). The top panel magnifies the range around Etot = 0.
thumbnail Fig. 7

Energy balance for the V4 model at AV = 0.75. The blue line shows the standard computation. The orange line was computed using a 20% enhanced H2 formation heating.

6 New model results

The most important surface reaction in the ISM is the formation of molecular hydrogen (Cazaux & Tielens 2002, 2004; Le Bourlot et al. 2012; Bron et al. 2014; Wakelam et al. 2017a; Thi et al. 2020) and all PDR models already account for it with various degrees of complexity. Other surface reactions were typically not included until sufficient computing power became available. (Woitke et al. 2009) presented the disk PDR model ProDiMo with surface chemistry and (Hollenbach et al. 2009) published a plane-parallel PDR model with a small surface network to study the chemistry of H2O and O2 and (Guzmán et al. 2011) studied the chemistry in the Horsehead using an updated version of the Meudon code (Le Petit et al. 2006). Recently, (Esplugues et al. 2016, 2019) presented PDR model results including (time-dependent) surface chemistry. In a more recent work, they expanded their chemical desorption scheme to include partial desorption of the products (Esplugues et al. 2019). We compare our model computations with their results and in addition we investigate how the revised description of CR induced desorption affects the structure of the PDR. As shown by (Sipilä et al. 2021) Eqs. (12) and (15) result in higher values of fCR compared to 3.16 × 10−19 as suggested by (Hasegawa & Herbst 1993) making the CR-induced desorption more efficient. We use the label CR1 and CR2 to indicate CR desorption by (Hasegawa & Herbst 1993) and (Sipilä et al. 2021), respectively. Chemical species and model parameters for computations in Sect. 6 are summmarized in Tables 6 and 7, respectively.

Table 6

Chemical species in the surface test models.

Table 7

Model parameter for models discussed in Sect. 6.

6.1 Selective freeze-out

The ice composition in a molecular cloud varies with extinction because it is mostly governed by thyhe dust temperature, which decrease with AV. The gas temperature does not always steadily decrease with cloud depth because cooling radiation may become trapped in the cloud for large optical depths, leading to increasing gas temperature for high values of AV. Table 8 lists the average condensations (dust) temperatures for the most abundant ice species, derived from the balance between accretion and desorption rates. It is important to note, that major carbon-bearing ice species condense below 50 K (with the exception of methanol). If Tdust is sufficiently high, this may selective favor oxygen-bearing ice components or lock up elements in certain ice species that are favored under current dust temperature conditions.

We compare how the structure of a strongly illuminated PDR (n = 104 cm−3, M = 103 M⊙, χ = 104) depends on the surface chemistry. Figure 8 shows the gas and dust temperatures of a PDR model cloud as function of AV. The dust temperature does not alter with adding surface chemistry. The gas temperature reaches a minimum at AV ≈ 5. At this depth, the FUV radiation is sufficiently attenuated and does not dominate the gas heating any more. Other processes, such as heating by gas-grain collisions start to dominate. Deeper in the cloud it becomes increasingly more difficult for cooling line photons to escape the model cloud due to optical thickness. This leads to an increase in gas temperature until it approaches Tdust. This effect is enhanced by any process that further limits the cooling capacity of the gas, for instance freeze-out of relevant line cooling species. Figure 8 shows how the gas temperature increase is stronger if surface chemistry is included and reaches Tgas ≈ 40 K, almost 20 K warmer compared to the model without freeze-out. This effect also occurs in case of the more efficient CR induced desorption model CR2.

Figure 9 shows the corresponding chemical structure excluding (left panel) and including (two right panels) surface chemistry. The incident FUV field in this particular case is strong leading to the typical transition from C+ → C → CO at AV ≈ 3. In the center of the model cloud without surface chemistry, the majority of the carbon is bound in CO with a temperature slightly above 20 K. The atomic carbon density peaks at AV = 2–5. In the models with surface chemistry, we find the same behavior up to an AV ≈ 7. At higher visual extinction, the CO abundance decreases while other gas phase species have higher abundances, for example C and CH.

We note that the dust is too warm to host a significant J(CO) population. Under these conditions the ice consists mainly of J(CO2) which locks up most of the available C and O atoms in the ice mantle. This renders the main destruction reactions of atomic carbon C, for example by collisions with O2 and SO2, inefficient. These reactions would otherwise replenish the CO population after destruction by H3+ and He+. Instead, the carbon remains locked in atomic form and the oxygen is converted to ice species. That affects related species like HCO+ and OH. The higher C abundance leads to an enhanced abundance of light carbon hydrides CHn. As a consequence CO is not available as coolant any more and the gas temperature rises by 15–20 K compared to the pure gas-phase model as shown in Fig. 8. This effect is more pronounced in models with lower FUV illumination where the additional warm C core is visible through 50–100% stronger [C I] emission at 610 and 370 µm.

For χ = 104, the relatively high Tdust prevents most ices to remain on the surface. Comparing the frost temperatures from Table 8 with Tdust from Fig. 10 shows that J(CO2) is the only viable carbon reservoir in the ice mantle that could survive temperatures above 40 K. CH3OH condensation temperatures would be sufficiently high, but the formation in the ice is prevented because the relevant precursor species do not survive long enough in the solid phase. There are no efficient chemical formation routes for CO2 → H2O available, which prevents the formation of significant amounts of H2O ice.

The influence of surface chemistry on the CO abundance has already been described by (Hollenbach et al. 2009). They modeled a significantly simpler chemistry and are mainly focused on the predictions for the H2O and O2 gas-phase abundance but we find a similar formation and destruction behavior for CO when we compare with our χ = 103 results. For the higher FUV field shown in Fig. 9, we find that for AV ≳ 5 CO is mostly destroyed in gas-phase models by H3+ + CO → HCO+ + H2 and He+ + CO → O + C+ + He. Adding surface chemistry does not introduce fundamental different channels with the exception of freeze-out.

The increasing C density at high visual extinction can also be understood by looking at the respective formation/destruction rates. Surface chemistry suppresses the dominant destruction channel of C at high AV: C + O2 → CO + O because most of the oxygen is locked up in J(CO2) and J(SO2) ice. This reaction can no longer replenish the CO population leading to a strong increase in C abundance and a decrease in CO in the dark cloud portion of the model cloud.

thumbnail Fig. 8

Influence of the surface chemistry on the thermal structure of a PDR. The model parameters are: n0 = 104 cm−3 (α = 1.5), M = 103 M⊙, χ = 106. CR1 and CR2 indicate the different models for the CR induced desorption.
Table 8

Condensation temperatures for some relevant species in K derived from the balance between accretion and desorption.

6.2 Illumination effects on the surface chemistry

For a more systematic analysis, we investigate how the additional surface chemistry affects individual species depending on the FUV field strength χ and whether the chemical desorption reactions discussed in Sect. 3.3 changes the overall chemistry.

The gas temperature remains mostly unaffected for χ ≫ 102 as shown in Fig. 11a. In Fig. 12 we show how different values of χ = 1, 102, 104, 106 affect the chemistry for C, HCO+, and O2. The high FUV case χ = 106 shows no significant effect of the surface chemistry. All relevant coolants remain abundant in the gas-phase due to the warm grain surfaces. For lower FUV field strengths we find that the gas temperature increases above the pure gas-phase case once the dominant coolant CO starts to deplete from gas-phase (see Fig. 11b). A dominant CO gasphase population survives until the dust temperature falls below 〈Tdust〉 ≲ 20 K (see Fig. 10) and desorption becomes too weak to replenish gas-phase CO. At lower χ = 1, 102 this occurs at: AV ≈ 0.2,1, respectively. At higher FUV this requires significantly higher values of AV. Any effect that desorbs J(CO) more efficiently will thus result in a lower gas temperature. This can be seen in Fig. 11b which compares the CO gas-phase abundance between the CR1 and CR2 models. In the CR2 model, the energy injected per grain as well as the grain cooling time are enhanced which effectively corresponds to a higher 〈Tdust〉. As a result we find that CO freeze-out sets deeper in the cloud if χ ≤ 102 and consequently we see a somewhat lower Tgas.

A general trend is that most carbon bearing species show enhanced abundances at high AV while most oxygen-bearing species have lower densities compared to the pure gas-phase chemistry. This is because ice composition is dominated by oxygen-bearing species, which locks-up a significant fraction of the available oxygen in the ice mantles. For models with significant CO freeze-out we find C to be the main carbon reservoir in the gas-phase. As a consequence, the CR induced ionization of atomic carbon can significantly contribute to the electron formation. Thus, we observe the same enhancement for the electron density. This strongly affects all species that form via dissociative recombination, for instance CH. We note similar inherited effects in Fig. 12, where HCO+ follows the abundance change of CO. The selective freeze-out of oxygen-bearing ice species also leads to a significantly decrease gas-phase abundance of O2. This is in agreement with observations having difficulties to confirm high molecular oxygen abundance predictions from pure gas-phase models.

Enhanced cosmic ray induced desorption in the CR2 models allows for a higher gas-phase abundance of certain species with lower binding energies, such as CO and CH4. The corresponding ice mantle composition is shifted toward tighter bound ice species, such as J(CO2) and J(H2O). A weaker dominance of overall depletion over gas-phase abundances seems to be in agreement with observations of for example S bearing species in PDRs that do not show signs of significant freeze out (Rivière-Marichalar et al. 2019). However, our current implementation of CR induced desorption (CR1 vs. CR2) remains a crude approximation with large uncertainties. Nevertheless, it might be feasible to use observations of low FUV sources as calibrators for more detailed models of CR induced desorption.

Chemical desorption does affect the formation and destruction of some species at various cloud depths. For instance, J(O) + J(O) → O2 takes over as the dominant O2 formation channel for 2 ≲ AV ≲ 5 so that we observe an enhanced O2 density before oxygen freeze-out starts to become important (Fig. 12c). Another example is the reaction J(O) + J(H) → OH which contributes approximately 10% to the total OH formation rate for AV > 30.

Figure 13 shows how the ice composition changes with cloud depth. Each panel in the figure corresponds to a different FUV illumination (assuming CR2). In the low FUV cases, χ ≤ 10 the ice forming closest to the cloud surface consist mainly of J(CO) (60–80%) up to a few AV. Deeper into the cloud the ice is converted to water ice (~50%) as well as methanol ice J(CH3OH) (~20–25%). Under increasing FUV conditions the ice composition shows some significant changes. For χ = 102, 103 we find J(H2O) ice closer to the cloud surface. J(CO) becomes less abundant (10–25%) and forms only deeper in the clump (AV > 2). For FUV strengths χ ≥ 103 the dust temperatures prohibit large amounts of J(CH4) and CO ice and the carbon ice consists to 50–90% of J(CO2), J(SO2) and water ice. For χ = 105 the ice is dominated by H2O ice and J(CO2) at AV > 50. The deeply embedded J(CO2) population vanishes at χ = 106 and water ice remains the dominant ice component.

The most obvious conclusion is that adding or removing some surface species from the chemical network may result in significantly different ice structures. The same is of course true for different sets of binding energies in use. Even so, another main effect of (any) existing surface chemistry is to open up new formation and destruction channels that become active below certain grain temperatures and will therefore significantly alter the gas-phase abundances. The freeze-out of CO is a good example for this effect. In terms of typical PDR tracers, it is important that they are removed from the gas-phase and not so much whether they end up forming J(CO) or J(H2O) ice.

The effective dust temperature is the most important factor in determining the ice structure. Figure 10 shows that assuming very low dust temperature for deeply embedded parts of the model PDR is not always justified. High FUV models may show enhanced dust temperatures and selectively prevent certain ice species to form. This modifies the overall ice composition on the grain surfaces also affecting the gas-phase abundances. As an example we showed how gas-phase CO is diminished in the dark cloud even though no explicit J(CO) ice is formed.

thumbnail Fig. 9

Influence of the surface chemistry on the chemical structure of a PDR. The results are for the same model as in Fig. 8. Solid (black) lines in the central (right) panel correspond to CR1 while dashed (red) lines are for CR2.
thumbnail Fig. 10

Dust temperature 〈Tdust〉 profile for n0 = 104 cm−3 and different FUV fields. Sublimation temperatures of selected ice species are indicated by red marks on the right.
thumbnail Fig. 11

Left: temperature profile changes with FUV strength (n = 104 cm−3). Solid lines show pure gas-phase results, dashed lines correspond to gas+surface chemistry. Right: influence of the surface chemistry on the CO gas-phase abundance for the same parameters as in Fig. 11. CR1 and CR2 indicate the different models for the CR induced desorption.
thumbnail Fig. 12

Chemical structure changes with FUV strength (n = 104 cm−3). Solid lines indicate pure gas-phase chemistry, dashed lines correspond to the gas+surface chemistry using the CR2 model.

6.3 Coupling between chemistry and line excitation

The previous section showed that the removal of CO from the gas-phase results in a reduced cooling capacity of the gas and an increase in gas temperature. In terms of observable line intensities the higher gas temperatures may partially compensate the reduced abundance of CO. Nevertheless, the general behavior is that 12CO and 13CO line intensities are lower in models with surface chemistry due to the freeze-out of CO. This is shown in Fig. 14 where we plot the CO line emission of the transition Jup → (Jup − 1) as function of the angular momentum quantum number Jup for pure gas-phase and surface chemistry models. We find the same behavior across the full density range.

Table 9 lists the clump averaged emission of the fine-structure lines. The selective freeze-out of O bearing species that leads to a dominant gas-phase population of C at large AV results in enhanced fine-structure emission for models with low to intermediate values of χ. This effect is ~2o% stronger in the CR2 models due to the slightly enhanced gas-phase abundance of atomic carbon. The enhancement of the [C I] intensities due to the surface chemistry occurs for the [C I] 3P1 → 3P0 and the [C I] 3P2 → 3P1 lines. The same effect does not occur for the [C II] line because of the higher energy of 92 K required to excite the [C II] 2P3/2 → 2P1/2 transition. The same holds for the [O i] lines at 63 and 145 μm.

Altogether, the surface chemistry in KOSMA-τ has a small effect on the total CO cooling budget, raising the gas temperatures by typically less than a few K (with some exceptions visible in Fig. 11a). The [C II] and [O I] cooling lines are rather insensitive to the modified chemical and temperature profile while the [C I] fine-structure lines may be affected depending on the effectiveness of overall desorption, which makes the atomic carbon an interesting tracer of surface chemistry in PDRs. We note that in diffuse and translucent clouds gasphase PDR models notoriously over-predict [C I] abundances (the so-called [C I] problem, (Gong et al. 2017). However, column densities are not a direct observable and therefore comparing model column densities to column densities derived from observations introduces additional uncertainties. Here, we are more concerned with denser clouds where PDR model prediction can over- and under-predict [C I] intensities depending on the detailed modeling approach. Standard chemical models are not able to reproduce the described [C I] enhancement because they do not solve the temperature self-consistently with the chemistry and the nonlocal radiative transfer through the model cloud.

The line intensities in Table 9 and Fig. 14 were computed for a model with n = 104 cm−3, M = 103 M⊙. Many observations of PDRs show significantly brighter emission lines. Examples are the Orion Bar and NGC 7023 (Joblin et al. 2018) and the Carina Nebula (Wu et al. 2018). They require significantly higher gas densities (n = 105−106 cm−3) to explain the observed intensities. At those high densities rotational lines of CO up to J > 20 can be excited; this is not possible for n = 104 cm−3. For comparison we provide plots of the CO SLED for higher densities in the appendix.

thumbnail Fig. 13

Percentile ice composition as function of AV for n = 104 cm−3, M = 103 M⊙ assuming CR2.
thumbnail Fig. 14

Effect of surface chemistry on the CO spectral line emission distribution (in K km s−1) as function of Jup. Each panel corresponds to a different FUV field strength χ. 12CO lines are in color, while 13CO lines are plotted in gray-level.
Table 9

Comparison of model fine-structure emission (in units of K km s−1).

6.4 Comparison with Other Models

Even though inter-model comparison tends to be difficult (compare for e.g., Röllig et al. 2007) we try to compare our results to other computations.

Hollenbach et al. (2009) added surface chemistry to their PDR model (Kaufman et al. 1999, 2006) to study H2O and O2 abundances observed in the gas-phase. They added only a small surface network to their chemistry and used significantly different desorption energies. In particular the values for C and OH are very different significantly affecting the ice composition. Nevertheless we find that H2O peaks at AV ≈ 5 with n(H2O) ≈ 10−7 cm−3 and drops to 10−9−10−10 cm−3 deeper in the clump, comparable to their result. We also see a corresponding increase of the water ice density at AV ≈ 3 − 5 locking-up most of the oxygen atoms in the ice mantle. On the other hand, they find a significantly different carbon ice structure, except for J(CO). For AV > 6 they find the carbon atoms locked-up in J(CH4) which does not occur in our models. They do not report details of their carbon surface network but the most likely reason for this is the assumed binding energies. Using comparable values we also find a dominant J(CH4) population. They discuss the threshold AV for ice formation (ML = 1) and we provide the corresponding numbers in Fig. 15. For conditions applicable to the Taurus cloud (n0 = 103 cm−3 and χ ~ 1) they find a value of 2 almost identical to our result of 2.1. For n0 = 104 cm−3 and χ ~ 100 they find approximately 3 where we have about 2.3 which can be explained by our increasing total gas density profile, which shifts the ML = 1 threshold closer to the surface.

Esplugues et al. (2016) presented a significantly improved version of the Meijerink PDR code (Meijerink & Spaans 2005) including an up-to-date surface chemistry. Unfortunately, they computed their models only up to AV = 10. More importantly, they used a simple approximation for their dust temperature that produces much too low Tdust (compared to detailed computations) and is not valid for AV ≳ 10. Therefore, their model is not able to selectively freeze-out particular ices only and thus could not produce the dark cloud C population we find. Their C density profile for n = 104cm−3, G0 = 104 at AV < 10 (Esplugues et al. 2016) was very similar to our results (see Fig. 12a). We note that typical plane-parallel PDR computations often do not cover deeply embedded regions with AV ≫ 10 and therefore miss out the strong peaks in species such as C and CH. The common argument that beyond a certain visual extinction no significant chemical variation takes place is not necessarily valid any more once we include surface reactions to the chemical network and account for higher dust temperatures due to diminishing cooling efficiencies. Looking at their Model 1 results (n = 104 cm−3, G0 = 104; Esplugues et al. 2016) we note that their ice composition is dominated by J(CO2) for AV > 3. Our models show J(CO2) as major ice component for slightly larger values of AV. They find water ice as second most abundant ice component for AV > 3. Our results indicate that water ice is dominating the ice composition at AV of a few with a J(SO2) peak around AV ≈ 5 which is not included in the model chemistry of (Esplugues et al. 2016). J(CO2) is dominating the ice composition at AV > 5. Figure 15 shows how the ice thickness (in MLs) in the center of our model clumps depends on the gas density and the FUV illumination. The number of MLs decreases with χ as well as with n0. The same is true for the depth where ML = 1. For n0 = 104 cm−3 and χ = 104 we find that the first ML builds at AV = 5.0 roughly consistent with (Esplugues et al. 2016) who find 0.6 ML at that extinction. However, our density increases with depth and therefore we have a slightly earlier onset of ice formation.

In a recent update, Esplugues et al. (2019) present computations with a modified dust temperature formulation (Hocuk et al. 2017) allowing for higher central dust temperatures. Comparing their dust temperatures with the detailed results from MCDRT shows a factor two lower values at the edge of the cloud but a factor ~2 higher central dust temperatures (except for χ = 106). At AV = 10 their Tdust is consistently ~ 10 K hotter except for χ = 1 where they have 8.8 K compared to our 3.6 K. As a result they also see a preferentially oxygen-bearing ice composition but with a different spatial behavior and a different ice composition because of their dust temperature exceeding the sublimation temperatures of CO and CH4. They present three models with different values for density and FUV strengths compared to Esplugues et al. (2016). A comparison with their results for the OH density profiles shows a similar behavior for AV < 5−6. The density peak shifts closer to the surface for increasing n and G0 and we observe the same trends. For n = 105 cm−3 they show a peak density of ni/nH ~ 10−8−10−7. We find similar densities for χ = 104. For the lower FUV illumination we find peak OH abundances of ni/nH ~ 10−6. Their model 3 (n = 106, G0 = 104) shows a much higher peak density for AV < 1 of a few 10−5 and we find similar values at AV ~ 0.2. However, all their models show a high OH density at their maximum AV = 10 of 2−3 × 10−8. Our models show a significantly lower central OH density of a few 10−10 in for n0 = 105 models. We find comparably high central densities only for χ = 106 which is consistent with our lower dust temperatures for χ < 106. Our result for their model 3 configurations shows similar densities at AV = 10 but differs significantly in shape and drops to 10−10 in the center of our cloud. For O2 we find significantly different results. At χ = 102, 104 we find peak relative densities of a few 10−6 at AV = 3−5. Our central O2 density drops with n0 to ni/nH < 10−12 for n0 = 105 cm−3 compared to values above 10−6−10−5 for their models 1 and 2. Such high values seem to be in conflict with observations (e.g., Goldsmith et al. 2011; Wirström et al. 2016, and references therein). Only for χ = 106 are we finding a strong central O2 population of larger than 106. A similar difference is also visible for H2O that has lower central and higher peak values in our results. (Nagy et al. 2017) gives relative H2O densities for the Orion Bar of ~1−5 × 10−12. The closest parameter set from Esplugues et al. (2019) is their model 3 with peak/center abundances of ~10−7. Within the same AV range we find factor ten lower values with central densities of 3 × 10−11 for n0 = 106 cm−3 and χ = 104. Our total H2O column densities at AV = 5 are 2 × 1013 cm2 for n0 = 105 cm−3 and 7 × 1013 cm2 for n0 = 106 cm−3, which is of the same order as the observed values of 2 × 1012 to 2 × 1013 cm2 (Nagy et al. 2017). Within their computations, changing the dust temperature prescription gave J(CH4) abundances at AV = 10 that differed by a factor 1016! Given their different dust temperature behavior and a significant different set of assumed ED they find a very different ice composition and we conclude that a detailed comparison with our results is difficult. Qualitatively, we find comparable amounts of J(CH4) and J(CO) ice for low values of χ but significantly lower abundances for higher FUV fields.

Guzmán et al. (2011) presented observations of H2CO emission from the Horsehead together with surface chemistry computations from the Meudon PDR code (Le Petit et al. 2006). H2CO is interesting because it can be efficiently formed in the gas-phase and on the surface of grains. They showed that adding the surface network leads to an increase of HCO densities by 1−2 orders and a pronounced density peak of H2CO at a few AV that is not produced in their gas-phase results. We find a similar qualitative behavior and can reproduce their H2CO peak densities (relative to ntot) of 10−8 for comparable PDR parameters (n0 = 105 cm−3, χ = 102).

thumbnail Fig. 15

Ice thickness in ML for the model clump center as a function of total gas density n0 and FUV strength χ, using CR2. The numbers give the threshold for the formation of one monolayer of ice in AV.
thumbnail Fig. 16

Ice abundances relative to the H2O ice column density. The median of all models with χ = 1−103 and M = 50, 1000 M⊙ is shown.

6.5 Comparison with Observations

Boogert et al. (2015) summarizes observed ice abundances based on column density derivations. The dust temperature variations along the observed lines of sight have a major effect on the ice but remain unknown which makes a direct comparison of the model ice composition with observations not trivial. Secondly, the limited number of chemical species included in the presented computations limits the predictive power for some of the ice species, for instance for J(NH3). Boogert et al. (2015) give ice abundances for different environments: massive young stellar objects (MYSOs), low mass young stellar objects (LYSOs) and BG stars. Naturally, they correspond to different local physical conditions, but we focus on general trends only. All abundances in this section will be relative to J(H2O) column densities unless stated otherwise.

The securely observed ices are in descending order: J(CO2) and J(CO) with 20–30% abundance each, J(CH3OH) with 6–10%, J(NH3) with ~6% and J(CH4) with ~5%. J(H2CO) is likely identified with a few % abundance. In Fig. 16 we show our model ice predictions where we average over a range of FUV strengths, χ = 1−103 and clump mass M = 50, 1000 M⊙ and plot as function of the gas density n0. J(H2O) is the most common ice across the whole density range and J(CO) varies between 20 and 40%. Both results are in agreement with observed numbers. We also find that J(CH3OH) is predicted with 3–20% relative abundance which is also consistent with observations. J(CH4) model abundances are higher than the observed ranges by a factor 2–3 and model J(CO2) is consistent with the observed 20–30% only for n0 ~ 103 cm−3. Our J(H2CO) predictions recover the observed values in the lower density range of Fig. 16. Given the crude averaging over the physical parameters and our small chemical network we find our model predictions to be in reasonable agreement with observations.

7 Clumpy Ensemble Model

Interstellar clouds are neither a plan-parallel slab nor of perfect spherical shape and results from these kind of models will always be a rough approximation to reality, a reality where the ISM is clumpy/fractal, turbulent, organized in filaments or fibers, and most importantly not in equilibrium. PDR models with more complex geometries have been designed to address this deficiency, but the higher complexity always comes with the price of much higher computation costs (Bisbas et al. 2012, 2021; Levrier et al. 2012; Grassi et al. 2014; Girichidis et al. 2016). The spherical setup of KOSMA-τ offers the attractive option of modeling clumpy clouds as a superposition of differently sized clumps following a well-defined clump-mass spectrum. This has been described in detail in Zielinsky et al. (2000), Cubick et al. (2008) and Andree-Labsch et al. (2017). For details see Appendix F.

Clumpiness has frequently been invoked to explain certain emission characteristics of the ISM in spatially unresolved observations. The main driver was always that clumpy gas has a higher surface to volume ratio and therefore shows an excess of emission primarily produced in the PDR surface regions of molecular clouds. A typical example is the [C II] emission, which is produced by strong FUV illumination and a good tracer of PDRs, this means a surface tracer. Conversely, rotational CO line emission is a good volume tracer because CO requires shielding from intense FUV illumination, which is effective only for extinctions AV > 1. Nonclumpy PDR models were not able to explain the observed excess in for instance the [C II]/CO(1−0) line ratio in active star forming regions (Stutzki et al. 1988; Spaans & van Dishoeck 1997; Dedes et al. 2010; Graf et al. 2012) but the observed ratios asked for models with a larger surface-to-volume ratio. A similar explanation has also been presented by Meixner & Tielens (1993), Hogerheijde et al. (1995), and Zielinsky et al. (2000) to explain variations in observed line ratios of several different PDR and molecular cloud tracers through clumpy ensembles of PDRs. Cubick et al. (2008) showed that the global far-infrared (FIR) emission of the Milky Ways can be explained in terms of clumpy PDR emission.

This has been supported by observations of clump mass spectra in molecular clouds. Based on molecular line observations Heithausen et al. (1998) measured the scaling relations for one source over several orders of magnitude and Kramer et al. (1998) measured the clump-mass distribution in various sources confirming a common power-law across all sources with power-law index 1.6 to 1.8, extending down to the resolution limit and clump masses as low as 10−3 M⊙ for at least two of the sources. These findings have been put at question by models that managed to explain the observed line ratios from high-pressure PDR models (e.g., Marconi et al. 1998). Joblin et al. (2018) showed that the pure line intensities in the Orion Bar can be explained from a plane-parallel PDR model. Their high-pressure models however, do not reproduce the observed spatial stratification of the different tracers but predict a very thin PDR layering. The key information on the clumpiness comes from the spatially resolved structures. Andree-Labsch et al. (2017) showed that the observed spatial stratification within the Orion Bar is in disagreement with any simple plane-parallel model but that some kind of clumpiness had to be invoked. Andree-Labsch et al. (2017) presented KOSMA-τ-3D, where individual volumetric elements (voxel) are populated with unresolved clumpy PDR ensembles. They modeled the 3D structure of the Orion Bar and found a good agreement of their results with a multiline data set from Herschel and Caltech Submillimetre Observatory (CSO) observations.

Velocity-resolved observations from Herschel, SOFIA and ALMA confirm the dynamical nature of PDRs (Goicoechea et al. 2016, 2017; Joblin et al. 2018; Wu et al. 2018; Luisi et al. 2021; Kabanovic et al. 2022). Photo-evaporation flows from globules and other dense clumps are ubiquitous (Mookerjea et al. 2012; Bron et al. 2018). Theoretical papers have predicted for a long time that dense clumps are carved out from their parental cloud by UV irradiation (Lefloch & Lazareff 1994; Henney et al. 2009; Bisbas et al. 2011). The inter-clump medium is then fed by the photo-evaporation flows from the PDR surfaces. The observations confirm the theoretical predictions that the inter-clump medium density is lower than the clump density by at least two orders of magnitude (Arkhipova et al. 2013; Mookerjea et al. 2019; Schneider et al. 2021). The photo-evaporation from the clumps provides a continuous low-density mass flow from the surfaces with velocities of 1−2 km s−1 (Makai 2015; Mookerjea et al. 2012, 2019; Sandell et al. 2015; Goicoechea et al. 2020). The inter-clump medium inherits the chemical properties of the PDR surfaces. Despite of the strong density difference their chemistry is rather similar because it is dominated by the FUV radiation rather than collisions. However, in principle this constitutes a nonstationary scenario due to the constant mass loss to the inter-clump medium (Bertoldi & Draine 1996; Störzer et al. 1997; Störzer & Hollenbach 1998, 1999; Maillard et al. 2021). It violates the assumption of a constant mass in the KOSMA-τ framework primarily affecting low-mass clumps (Mcl ≲ 10−2 M⊙ for n ≳ 106 cm−3; Decataldo et al. 2017, 2019). As a dynamic effect it is currently also not contained in isobaric PDR models where the low density gas only occurs as a very thin and hot surface layer.

The picture of the mass flow constitutes a recipe for us to represent the clump vs. inter-clump structure within the clumpy model framework. In KOSMA-τ we can model the dense PDR clumps by KOSMA-τ clump masses Mcl > 10−2 M⊙ while the inter-clump gas is represented by UV dominated conditions that are well modeled by low mass clumps with small AV.

With the continuously improving spatial resolution of modern observatories, such as ALMA and IRAM/NOEMA, it becomes possible to further test the clumpy ensemble picture by new observational data. However, even with today’s interferometric observations it remains difficult to resolve clumps with masses below 0.1 M⊙ in continuum observations. The completeness limit of the clump analysis is often reached at a few 0.1 M⊙ so that the slope of the mass distribution below one solar mass is uncertain and highly debated (see e.g., Pineda et al. 2009; Miville-Deschênes et al. 2017; Cheng et al. 2018; Kong 2019; Könyves et al. 2020). For an overview of how different clump finding strategies affect the derived size and mass distribution see Schneider & Brooks (2004) or Li et al. (2020). Moreover, the photo-evaporation modifies the clump distribution in molecular clouds under the irradiation forming PDRs. Hence, it is an important question to study whether the observational data provide further constraints on the slope in the clump size spectra below the scale resolved in the continuum observations.

7.1 Example: Orion Bar – Spatial Structure

High-resolution ALMA data of the Orion Bar directly resolve the larger clumps in our description. Goicoechea et al. (2016) presented ALMA observations of HCO+ 4−3 emission lines of the Orion Bar resolving the structure of the ionization and dissociation front with ~1″. Their data show fragmented clumps with sizes in the order of 2″ (see Fig. 17). This corresponds roughly to a clump mass of about 0.01 M⊙ at densities of n ≈ 4 × 106 cm−3 and ⊘ = 1.44″. KOSMA-τ computes Icl(HCO+4 − 3) = 57 K km s−1 clump averaged emission for such a clump assuming a FUV field of χ = 104. For a detailed discussion on suitable PDR model parameters for the Orion Bar we refer to Andree-Labsch et al. (2017). For a qualitative discussion we overlay in Fig. 17 contours for 60 K km s−1 to the HCO+ 4−3 data presented by Goicoechea et al. (2016). The comparison with the size of a 0.01 M⊙ clump in the lower left corner (radius 1.44″) shows a match of the typical structure size. Most of the emission at lower levels is not spatially resolved but forms an extended structure, mainly behind the Orion Bar. If we assume that this is due to smaller clumps, well below the beam size and not removed by photo-evaporation yet, we can also model them through KOSMA-τ. When following the original clump mass spectrum clumps with a mass of 0.001 M⊙ have a radius of 0.53″ and a clump averaged HCO+ 4−3 intensity of 39 K km s−1. For the comparison, we also include dashed contours for 40 K km s−1 in Fig. 17. Comparing the areas within the two contours we can count the required clumps to produce the observed emission and compare the ratio with our standard clump ensemble scaling.

Our standard clump ensemble setup approximately reproduces the properties of this HCO+ intensity map without any further fitting. Table 10 lists the parameters of the corresponding discrete clump ensemble. Here, we use just the two clump masses of m1 = 0.001 and mu = 0.01 M⊙ depicted in Fig. 17 and assume the standard scaling laws for the clump mass distribution (power law index α = 1.8) and the mass-size relation (power law index γ = 2.3; Heithausen et al. 1998). The flux of the 0.01 M⊙ clumps fully explains the 57 K km s−1 contour if we ignore higher intensities within the area. This is an obvious oversimplification since the map shows few smaller condensations with K km s−1. The low mass clumps explain the outer contours to about 50%. This is a remarkably good match given the fact that we did not perform any numerical fitting.

To fully explain the intensity area it would take about twice as many clumps with M = 0.001 M⊙ (numbers in parenthesis in Table 10). Consequently, this results in a steeper clump-mass distribution with a power-law index of α ≈ 2.1 or contributions from more even smaller clumps that were ignored in this simple picture. A steeper clump-mass distribution would agree with previous studies showing that denser regions, in particular in Orion A, tend to have steeper clump-mass indices a comparable to what we find (Bally et al. 1987; Maddalena et al. 1986; Nagahama et al. 1998; Schneider & Brooks 2004).

We can also compare the model column densities. For the model clumps from Table 10 we compute the mean column density of a species averaged over the projected clump area. The contour area filling factor of approximately unity allows to directly compare the column densities. Table 11 gives the observed values and our clumpy results. The observed column densities are mostly consistent with a clumpy PDR ensemble. The CO column density predicted by the model is higher than the value derived from 3−2 observations by Goicoechea et al. (2016). This is easily explained by the fact that the derivation used there is not sensitive to CO at temperatures below 120 K so that a large fraction of cool CO was not accounted for. The major discrepancy in the SO abundance may be attributed to some large uncertainties in the sulfur chemistry due to the unknown roles of vibrationally excited H2 (Goicoechea et al. 2021) and direct depletion (Fuente et al. 2016). The total column densities are also higher than the values given by Goicoechea et al. (2016) but closer to estimates by Hogerheijde et al. (1995) who derived at the peak of molecular emission.

The simple numerical experiment shows that (1) the HCO+ 4−3 intensity levels predicted by KOSMA-τ clumps are consistent with the observations of the Orion Bar. (2) The total flux predicted by a clumpy ensemble covers the observed flux values to a significant degree. (3) The column densities from a clumpy PDR ensemble are consistent with observed values. (4) The fact that we do not resolve the smallest ensemble clumps is consistent with their filling factor which leads to roughly homogeneous intensity distribution.

thumbnail Fig. 17

Integrated intensity map of the HCO+ 4−3 line emission in the Orion Bar (see Goicoechea et al. 2016, for details on the mapped area). The contours corresponds to integrated intensities of 40 K km s−1 (dashed) and 60 K km s−1 (solid). The colors indicate intensities between 0 and 100 K km s−1. The two circles in the lower left show the area of clumps with M = 0.01 M⊙ (white) and M = 0.001 M⊙ (gray).
Table 10

Properties of a simple discrete Orion Bar ensemble.

Table 11

Observed column densities vs. model results.

7.2 Nonstationary Clumpy PDR Evolution

The structure in Fig. 17 is the result of a high-pressure zone moving through the molecular cloud (Goicoechea et al. 2016) forming a fragmented and dynamically moving PDR surface. Compression by this wave and photo-evaporation lead to enhanced density contrasts with dense clumps, subject to erosion, and a thin interclump medium (Gorti & Hollenbach 2002). Goicoechea et al. (2016) find a dynamical crossing time for the wave front of a few 104 yr. 104 yr corresponds to the photo-evaporation destruction time for clumps with M < 0.01 M⊙ (Decataldo et al. 2017) and the chemical time-scale for HCO+ formation or destruction in the clumps from Table 10. Consequently, we do not expect the survival of any smaller clumps after the passage of the PDR zone. This is consistent with the fact that Fig. 17 shows almost no 40 K km s−1 emission on the side facing the ionization front (δx ≲ 15″). Clumps with M < 0.01 M⊙ did not survive the passage of the pressure front. The regions in Fig. 17 with δx ≳ 20″ shows volumes which have not yet been affected by the front and where the lower mass clumps are still surviving. They may already be subject to some evaporation but their embedded HCO+ population is not yet diminished. This picture also explains the emission peaks with levels above 60 K km s−1 at the PDR front in Fig. 17. Clumps with original masses above 0.01 M⊙ are brighter than 70 K km s−1 but may have been partially eroded so that we only see their larger fragments. We conclude, that the findings by Goicoechea et al. (2016) are consistent with a clumpy medium with a mass spectrum that changes from the original distribution by photo-evaporation during the passage of the PDR. The observed gas dynamics and the related life times of such clumps explain the spatial distribution of the observed emission.

8 Conclusion

We present the current status of the numerical PDR model code KOSMA-τ, which solves the coupled chemical and physical state of the ISM in a spherical model cloud under isotropic FUV illumination. KOSMA-T has been thoroughly tested in the last 20 yr and has been subject to a number of significant updates and improvements. In this paper, we provide a detailed description of the geometric and numerical setup of the code as a reference for comparison with other model codes. Details of adaptive stepping and numerical convergence are discussed and recommendations for application in other codes are derived. We discuss a series of numerical modifications of the code to improve the convergence stability and overall performance. The implementation of a time-dependent chemical solver allows for an efficient fall-back algorithm in case regular steady-state solvers do not converge. Finally we present a publicly available sandbox PDR model realized in the programming language Mathematica by Wolfram Research.

KOSMA-τ now includes a complete surface chemistry network fully coupled to the gas-phase chemistry via all relevant accretion and desorption processes. In particular, chemical desorption has been added to KOSMA-τ. We extended the prescription presented by Minissale et al. (2016) and included all possible reaction branches, also including partial desorption. We implemented the surface chemistry in a fully modular way analog to the gas-phase chemistry, allowing for simple addition or removal of chemical species. Adding surface chemistry to a PDR model can produce unexpected results. In particular, the detailed treatment of the dust temperature is important because PDRs can sustain large amounts of warm dust particles at very large extinctions strongly affecting the ice chemistry. We discuss how higher dust temperatures lead to a selective freeze-out of oxygen-bearing species compared to carbon-bearing species due to their higher condensation temperature. For high FUV models, the reduction in the CO column density due to selective freeze-out of oxygen-bearing species onto grain surfaces produces a surplus of atomic carbon, leading to a significantly enhanced [C i] line emission that could explain the difficulties of previous generations of PDR models in fitting observed levels of [C I] emission. The composition of the ice mantles changes in a complicated, which is not always intuitive way as a function of the FUV irradiation. This applies, for example, to the nonmonotonous production of methanol that is suppressed at elevated FUV fields. We discuss the ice composition under different FUV conditions and compare our results with observations and with predictions by other models. We used recent ALMA observations of the Orion Bar PDR to test the clumpy PDR model assumptions. Performing a simple numerical experiment shows that the PDR structure as well as the observed flux values can consistently describe the assumption of a clumpy ensemble of PDR clumps. Clumps with masses below 0.01 M⊙ are eroded in the dynamical evolution of the PDR, explaining the asymmetry of the observed profile around the high-pressure zone when taking the observed timescales into account.

Acknowledgements

We wish to thank J. R. Goicoechea for providing us with the ALMA data. We also would like to thank the anonymous referee for thoughtful suggestions and constructive feedback that helped to significantly improve the paper. The research presented here was supported by the Collaborative Research Centre 956, subproject C1, funded by the Deutsche Forschungsgemeinschaft (DFG), project ID 184018867.

Appendix A New model results

We show the spectral line energy distribution (SLED) of 12CO and 13CO lines as function of the FUV field strength χ for n0 = 103, 105,106 cm−3 in Figs. A.1, A.2, A.3, respectively. The effect of the surface chemistry on the line emission is most prominent for high density and higher-J transitions.13 We also show the density profile of selected species affected by the surface chemistry.

thumbnail Fig. A.1

Effect of surface chemistry on the CO spectral line emission distribution (in K km s−1) as function of Jup for n0 = 103 cm−3. Each panel corresponds to a different FUV field strength χ. 12CO lines are in color, while 13CO lines are plotted in gray-level.
thumbnail Fig. A.2

Same as Fig. A.1 for n0 = 105 cm−3.
thumbnail Fig. A.3

Same as Fig. A.1 for n0 = 106 cm−3.
thumbnail Fig. A.4

Chemical structure changes with FUV strength (n = 104 cm−3). Solid lines indicate pure gas-phase chemistry, dashed lines correspond to the gas+surface chemistry using the CR2 model.
thumbnail Fig. A.5

Chemical structure changes with FUV strength (n = 104 cm−3). Solid lines indicate pure gas-phase chemistry, dashed lines correspond to the gas+surface chemistry using the CR2 model.

Appendix B Numerical details

Sect. 3.6 described the spatial loops over the adaptive cloud depth grid. Input to the Bulirsch-Stoer stepper are the dependent variable vector and its derivative at the starting value of the independent variable z = 0. Also input are the step size to be attempted, htry (we start with 106 cm), the required accuracy, eps, and the vector yscal against which the error is scaled to give the desired accuracy ∆0 = eps × yscali. On output, N and z are replaced by their new values. hdid is the step size that was actually accomplished, and hnext is the estimated next step size. By default KOSMA-τ uses eps = 0.01 and yscali = ∣Ni∣ + ∣htry × ni∣ + 10−21 cm−3 × htry (see Press et al. 2007, for details).14 It is also possible to alter the value of htry depending on local conditions to influence the performance of the iteration in terms of accuracy and total computation time.

A disadvantage of the adaptive stepping is that testing for global convergence is not straight forward. In contrast to models with fixed spatial grids one cannot easily compare densities and population numbers between each iteration and check whether a given numerical tolerance is met because the spatial grid changes between iterations. One obvious solution is of course to define a spatial reference grid and interpolate the resulting structure of each iteration to this grid. However, this introduces the additional error of the interpolation to the numerical uncertainty of the solution. We nevertheless offer this option in a coming version of the code.

B.1 Chemical Solver Details

Steady-state chemistry. KOSMA-τ offers a variety of solution approaches to solve Eq. 28 including LU decomposition (Lapack routines DGESV and DGESVX, Anderson et al. (1999)), the minimum norm-solution to a real linear least squares problem: min using the singular value decomposition (SVD, Lapack routine DGELSD) as well as the linear least squares solver DGELS. We note that DGESV is the workhorse among these routines and the others provide a fallback option in case the standard approach fails.15 We summarize the chemical solution in Algorithm: Chemistry 1.

Time-dependent chemistry. KOSMA-τ is by default computing the chemical steady-state solution as explained above. In addition, we also implemented the possibility to approximate the local equilibrium solution by a time-dependent solution of the chemistry using a long time tequil. We use the DLSODES (double precision) solver from the ODEPACK package (Hindmarsh 1983) to solve the system of ODEs assuming the previous chemical solution as initial condition.16 when the steady-state solution can not be found by the code (see Algorithm: Chemistry 1). This approach makes the code more stable but reduces the overall performance in terms of computation time. Presently we do not store the individual time steps and only use the solution at the equilibrium time. We are currently implementing a fully time-dependent solution in KOSMA-τ.

Algorithm: Chemistry 1 Solution

Tg → Tg,new
n(i) → n(i)old
while iterationcounter < 600 do    → Newton-Raphson
iteration
   compute
   compute
   apply preconditioning
   solve for → use e.g., DGESV to invert →
   converged → TRUE
   foreach species i do
    if n(i) > 10→33 cm→3 then
     →i → →ni/n(i)
     if →i > 0.05 AND n(i) > 10→15 cm→3 then
      converged → FALSE
    end if
    n(i)new → n(i)/fstep(→i)       → see Eq. 29
    else
     n(i)new → 10→33 cm→3
    end if
    OPTIONAL: sanity checks on n(i)
   end foreach
   if converged then return → chemistry converged
   end if
end while
apply fallback strategy, e.g., call DLSODES    → not converged
       → not converged

B.1.1 Benchmark of the Chemical Solver

In this section we compare how the different approaches to solve the chemical network perform in terms of total CPU time and overall solution accuracy. We solve the eight benchmark problems from Röllig et al. (2007) and compare the individual results against each other and against the PDR benchmark results. We use the set of model parameters described in Röllig et al. (2007) and the same chemical database with the addition of state-to-state formation rates of CH+ and SH+ (Agündez et al. 2010; Nagy et al. 2013). Changes to the past KOSMA-τ setup are the updated UV radiative transfer and dust temperature computation (Röllig et al. 2013), H2 self-shielding according to Draine & Bertoldi (1996), CO self-shielding according to Visser et al. (2009), updated heating and cooling computations, and improved numerics. Surface chemistry is ignored and H2 formation is simplified as described in Röllig et al. (2007).

Table B.1 shows the total CPU time used by the different model setups for the eight benchmark problems.17 The numbers in parenthesis give the results for runs where the time-dependent solver DLSODES was used as backup of the stationary solver in case the local chemical solution did not converge.

For the models with fixed temperatures, the choice of the steady-state solver does not significantly affect the total CPU time; all models finish with in about half a minute. The fast model convergence results from removing the requirement of a self-consistent temperature solution for the problem. Fig. 2 shows that this simplifies the “Local iteration” to solving the chemistry and the energy level populations. As a result the global iteration reaches convergence after about four iterations. Computing the time-dependent approximation to the steady-state solution converges also after 4-6 iterations. Models F1, F3, and F4 required about 104 total chemical iterations while F2 required > 105 iterations.

The variable temperature models show a larger variation in CPU times. In general, the time to find a stable solution is 10-1000 times longer compared to the “F” models. For almost all problems, allowing the code to use a time-dependent solution as fallback in case of convergence difficulties decrease the total CPU time by a significant factor. This was the case for the V1 problem using DSGESV. The choice of linear system solver moderately affects the total CPU times. Using DGESVX yields CPU times up to a factor two longer compared to DGESV, except for V1. CPU time consumption of DSGESV is also higher than DGESV. Across all solver, the problem V1 required the most CPU time. The time to immediately compute a time-dependent solution (last column in Table B.1) instead of the steady-state solution is growing from V1 to V4. For the V4 models DLSODES is about 20 times slower than the other solvers. We can summarize that the approach of using the DLSODES as fallback solver in case the linear solvers fail is successful and leads to a stabilization of the system, in particular if one considers the additional nonlinearity introduced by the choice of the chemical network, as described in Sect. B.2. It is also surprising that it is difficult to predict computation times depending on the parameter regime especially across different numerical solvers.

B.2 Numerical Aspects of Solving the Chemical Problem

Few publications on numerical details of solving the chemical network in the astrochemical context are available in the scientific literature. Nejad (2005) gives an extensive overview over available numerical solvers and methods to improve stability and speed of the solution. An important quantifier is the computative complexity of the solving algorithm, that is the number of computational operations necessary to solve the problem. In terms of solving Eq. (28), we need to invert ℚ in order to determine the new density vector . Inverting a matrix using classical Gaussian elimination can be done in steps. We note that this is also the cost of classically multiplying two matrices. More advanced algorithms can perform matrix multiplications in fewer steps. For example, Strassen (1969) presented an algorithm that performs matrix multiplication with complexity . For large problems this is a non-negligible improvement.

For the LU decomposition applied in the routines DGESV and DGESVX we find an approximate complexity of . Please note, that the number of chemical reactions does not influence the time to solve the problem, but it does influence the time to set up the system Eq. (7). Since the rate coefficients involve expensive operations such as exponential functions and fractional powers, increasing the number of reactions might lead to increased CPU running times of the code. This can be partially compensated by precomputing and storing computationally expansive expression that do not change during the chemical iteration. In the following we discuss various strategies to reformulate the chemical problem in a mathematically equivalent but numerically more favorable way.

Table B.1

Total PDR benchmark CPU time (in units of s) with different chemical solvers.

B.2.1 Preconditioning of the Jacobi Matrix

Matrix preconditioning in the context of solving stiff systems of ODEs is a wide and technical topic. Nejad (2005) gives an extended discussion on preconditioning. Here, we would like to present a simple yet sometimes helpful preconditioning strategy that can be used in KOSMA-τ if numerical convergence is problematic, for example due to round-off errors. With preconditioning, we mean a transformation of the system from Eq. 28 to the following form (Viallet et al. 2016): (B.1)

where and and ℝ are the left and right preconditioning matrices. Various choices for and ℝ are possible and KOSMA-τ employs the following definition: where specifies a diagonal matrix with vector as diagonal elements and all other elements set to zero. In index notation we can write: (B.2)

with the row index j of the right hand side vector and Lj the nonzero element in row j of . The tilde denotes the preconditioned version of the original vector/matrix. The Jacobi matrix is computed with: (B.3)

where Rj the nonzero element in row j of ℝ. The new system Eq. (B.1) solves for and we have to multiply the solution with ℝ, that is with to get the final solution.

Another strategy discussed in Nejad (2005) is to perform a row-reorder in the system Eq. (28). The idea is to make use of the sparsity of the problem, which is common in astrochemical computations because most chemical species react with few reaction partners and because many reaction rates are effectively zero in many circumstances. It is numerically favorable if large submatrices of the Jacobi matrix are null matrices. Modern solvers of systems ODEs and linear problems are often optimized to invert sparse matrices efficiently.

If we reorder the rows for example by the descending count of zero elements per row in the Jacobi matrix we get a new matrix with large submatrices with zero elements which might improve solution speed. KOSMA-τ allows to perform a row reordering such that rows of ℚ with the most zero elements are at the top and the densest rows at the bottom. The inverse scheme would also be an option and we advice any interested modeler to experiment with either scheme.

Because of the high magnitude difference in chemical densities of the involved species we find a large variance in the magnitudes of the Jacobi matrix entries. This might reduce efficiency of the solution and re-scaling the system is an option to somewhat assist the algorithms. We offer the option to rescale every row in Eq. (26) such that the maximum entry in the respective Jacobi matrix row is unity.

In cases where KOSMA-τ tries to approximate a steady-state solution by solving the time-dependent system for a sufficiently long time tequil, we apply the additional heuristic scaling factor 1/(10−2tequil) to the row scaling. We mention this here as experimental approach and invite other groups to perform additional tests.

Results. A performance analysis of the introduced numerical tweaks in full detail is beyond the scope of this paper. Their performance and their applicability depends on detailed model parameters. We performed test computations against a reference model of n0 = 104 cm−3 with pure gas-phase chemistry solving for the steady-state without any preconditioning. The following scenarios have been tested:

  1. preconditioned Jacobian as described in Sect. B.2.1

  2. Jacobian and r.h.s vector rescaled to unity first AND precondition the rescaled Jacobian

  3. the (already) preconditioned Jacobian and r.h.s vector rescaled to unity

  4. resorted Jacobian with dense rows first

  5. rescaled AND resorted Jacobian with dense rows first

  6. rescaled AND resorted Jacobian with sparse rows first

  7. adaptive choice of elemental conservation rows in ℚ, replacing the most abundant species

  8. adaptive choice of elemental conservation rows in ℚ, replacing the least abundant species

In Table B.2 we summarize the number of global iterations for each of the described variations plus the total computation time. The first row shows the reference computation applying no preconditioning. Preconditioning the Jacobian with and ℝ (strategy 1) shows the strongest impact in terms of performance while all others lead to comparable computation times or even a much worse performance. This will also depend on the detailed chemical network used and the parameter regime. It demonstrates though that the computational effort to solve the chemical problem depends strongly on the detailed numerical approach even though they are all mathematically equivalent. From the CPU times shown in Table B.2 we conclude that the LR-preconditioning outlined in Sect. B.2.1 performs best. Currently, the standard setup in KOSMA-τ is no preconditioning until a more systematic parameter study for typical PDR model parameter ranges has been performed.

Table B.2

Performance comparison of various preconditioning strategies of the chemical problem

Appendix C Newton-Raphson stepping

Historically, KOSMA-τ used the following ArcTan stepping: (C.1)

with λ = 0.9. The scaling by (2/π) limits the second term to the range ± 1 and the prefactor of λ = 0.9 limits fstep,A to the range 0.1,…, 1.9 and dampens the slope around the origin to 0.9. The stepper in Eq. C.1 has several advantages over the linear stepper in Eq. 29: it prohibits steps into the negative domain and ensures positive densities; it is almost Newtonian for small relative steps η and therefore converges fast once it is sufficiently close to the root; it prohibits too wide and too small steps because 1 − λ ≤ fstep,A(η) ≤ 1 + λ ∀η (seeFig.C.1).

An alternative approach is a modification of Eq. 29 to only prevent Newton steps through the negative domain: (C.2)

By construction guarantees a positive stepping factor Eq C.2 converges to zero for large negative steps η ≪ −1. It is also symmetric in the log-domain, for instance for η = 1 we find n(i)new = 2 × n(i)old, and conversely η = −1 returns n(i)new = n(i)old/2. However, this stepper diverges for large positive steps, η ≫ 1 and the symmetry between negative and positive steps may lead to oscillations around the actual solution.

Inheriting from both approaches we construct a more general stepping function fstep,Tanh (see Eq. 30) whose symmetry and limits can be controlled with few numerical parameters. The min/max values of fstep,Tanh are given by (C.3)

Consequently, ω+ controls the maximum factor to increase n(i)old and ω− controls the minimum factor for decrease. Choosing slightly different values for ω+ and ω− avoids oscillations in symmetric problems.

The slope ∂fstep,Tanh/∂η for small η, that is close to the solution, is given by λ. A value close to unity guarantees the quadratic convergence behavior of the Newton-Raphson method for small values of η. However, far from the minimum of a function f, a full Newton-Raphson step not necessarily decreases the function. We only know that the stepping direction initially decreases f. Reducing the step width but a damping parameter λ < 1 then limits nonconvergent steps18.

thumbnail Fig. C.1

Comparison of various stepper functions. The inset shows the vertical axis in log-space. fstep,Tanh1 assumes (ω− = ω+ = 2, λ = 1), fstep,Tanh2 assumes (ω− = 20, ω+ = 5, λ = 0.5).

Figure C.1 compares the discussed fstep variants. The log-plot inset shows the symmetry of fstep,Tanh that can be controlled by with ω+ and ω−. We see how the Tanh remains closer to the linear Newton behavior for larger η compared to the ArcTan. We note that the ArcTan stepper fstep,A can be reasonably well approximated with ω− = 11, ω+ = 2, and λ = 0.9.

C.1 Choosing Efficient Stepping Strategies

Ideally, one would like to find a single optimal stepping strategy to reliably solve the Newton-Raphson scheme. In the numerical PDR context this turns out to be a demanding task. We tested a large variety of different steppers and compound stepping strategies. We computed a small grid of models with 107 species including surface species (see Table 6). The model density (constant), mass and FUV field were varied on the following regular grid: n = (102,104,106) cm−3, M = (50, 100) M⊙, and χ = (1, 103, 106) χDraine. We set ζCR = 5 × 10−17 s−1 and assume RV = 3.1. The maximum number of Newton-Raphson steps is set to 600 per solution attempt. Global convergence was determined as described in Sect. 3.6.1.

At first, we test the performance of different choices of (ω−, ω+, λ) of the Tanh stepper and compare it against the old stepper fstep,A. Table C.1 shows the total count of Newton-Raphson steps per computed model. Each row corresponds to a particular stepper choice. fstep,A is the standard ArcTan stepping. The last seven rows correspond to a compound strategy that is discussed below. Models that failed to converge and aborted are marked with a dash. Reasons for abortion are NaN values in the chemical solution or an array overflow in case of requesting more than 1900 spatial grid points because of too steep chemical gradients.

For some settings we also allow for a more aggressive Newton-Raphson approach where we keep the chemical vector after the final iteration even if the convergence criterion for the local chemistry was not met (indicated by the subscript ()new in Table C.1). This is a fragile approach because it may fail if the Newton-Raphson stepping ends far from the true solution, but it may speed up convergence if the steps approached the true solution. The standard approach is to discard the (not converged) chemical vector and use the previous solution. This ensures that a ”true” solution of the chemistry is used.

There is no superior stepper and their individual performance depends on many details of the model numerics and the model parameters. We found that asymmetric choices ω+ < ω− perform better than symmetric setups. We note that ω+ > ω− is generally not recommended, for example KOSMA-τ frequently failed to finish computation because of density overflow in the case of ω+ = 5, ω− = 2. The push toward larger densities leads to numerically unstable behavior. The general trend shown in Table C.1 is the number of steps growing with model density and with model FUV strength. There is no clear correlation with the model mass.

As λ determines the slope of the stepping function for small arguments one would assume that solutions are found more quickly for λ closer to unity. However, the values from Table C.1 do not generally confirm this assumptions, we find many counterexamples. At first glance the total number of steps shows no clear correlation with λ but we note that missing models distort the row total, for instance (10, 5, 0.5) and (10, 5, 0.6) did not converge for two high density models each. These models typically took significant number of steps to finish and are not included in the row total. When approximately accounting for the missing models there is a slight performance improvement with increasing λ, excluding the case of λ = 1. Generally, symmetric choices of ω− = ω+ tend to perform worse than runs with ω− > ω+.

A general conclusion from Table C.1 is that for low to intermediate densities the details of the Newton-Raphson stepping are hardly affecting the number of steps with the exception of a few cases that failed to converge. The high density model performance shows a much higher variance with stepping, with some models performing up to ten times worse than others and some even more extreme cases.

Computations where we did not discard the density vector after the last Newton-Raphson step even if convergence was not yet reached (subscript ()new) led to more nonconverging models overall but successfully reached convergence in some cases where the standard approach failed. In particular the choice of ω− = 11, ω+ = 2, and λ = 0.9 seems to perform better in the ()new case.

For a given set of PDR parameters it is possible to find an optimized set of stepping parameters (ω−, ω+, λ) that will minimize the computing time. Changing the PDR parameters may significantly change the topography of the Newtonian vector field and consequently may require a modified stepping strategy to succeed. To mitigate these complications we tested a compound strategy, named adaptive 1 in Table C.1. In this approach we change the stepping parameters along the Newton search to find a chemical solution. Our general approach is to use a relatively robust stepping strategy for the first 100 Newton-Raphson steps. Typically, 99.99% of all solution attempts require less than 30-40 steps. If the standard stepper is not successful after 100 steps, it will most likely not converge at all, for example because of oscillating or global nonconvergent behavior. To dampen oscillations we then progressively reduce step widths. First we reduce ω− for the next 100 steps, followed by stronger dampening by reducing λ. If the more careful stepping is unable to escape a local minimum we allow for a more aggressive stepping, for example by allowing for larger values of ω+,− or even ω− < ω+. Algorithm adaptive 1 shows the details of the adaptive strategy.

adaptive 1

start from
if # steps → 100 then
   use fstep,A
else if 101 → # steps → 200 then
          → start again from
   use fstep,Tanh with (→→ = 5, →+ = 2, → = 1.0)
else if 201 → # steps → 400 then
   use fstep,Tanh with (→→ = 5, →+ = 2, → = 0.5)
else
          → start again from
   use /step,Tanh with (→→ = 5, →+ = 5, → = 1.0)
end if
Table C.1

Total Newton-Raphson step count for different stepper choices.

In its final stepping attempt adaptive 1 applies a more aggressive stepper setup. In addition, it restarts the root search from the initial density after 100 and again after 400 steps.

The large number of Newton steps computed in some models highlights the importance of the stepper choice. Comparing the adaptive and the single stepper strategies shows better performance for the adaptive strategy across all PDR parameters. adaptive 1 performs best with 1/2 total Newton steps for the whole grid compared to the reference stepper fstep,A. From Table C.1 we conclude that adaptive 1 is the most successful algorithm.

Appendix D Chemical network

D.1 Gas-Phase Chemistry

The gas-phase chemistry in KOSMA-τ is based upon the UDfA2012 database with updates listed in Table D.1. It contains the following reaction types: For two-body reactions the Arrhenius form for the rate coefficient is (D.1)

with the parameters α, β, γ tabulated in the database. UDfA12 distinguished between the following two-body reactions: neutralneutral, ion-neutral, charge exchange, ion-ion neutralization, dissociative recombination, radiative recombination, associative electron detachment and radiative association. In addition it contains two direct photo-processes: photo-dissociation and photo-ionization both described by: (D.2)

with the photo-rate in the unshielded interstellarultra-violetradiation field P0, the unit-less (self-)shielding factor fss, and the extinction by interstellar dust at visible wavelengths AV.

Direct cosmic-ray ionization is described by (D.3)

and cosmic-ray induced photo-processes with a rate (D.4)

where ζH2 is the direct CR ionization rate per H2, γ describes the probability per CR ionization that a photo-process takes place and ω is the dust grain albedo in the FUV (McElroy et al. 2013). We assume ω = 0.5.

For reactions with isotopologues we take the same rate coefficients as for the standard isotopologue with the addition of the fractionation reactions listed in Table D.2. Reaction rates for the reactions c−C3H2 + HE+, c−C3H + C+, c−C3H + H+, and c−C3H + He+ have been replace with rate coefficients from KIDA (Wakelam et al. 2015). New reactions involving the linear isomers l−C2H3+, l−C3H2+, l−C3H2, l−C3H have also been taken from KIDA.

Table D.1

Changes and addition to the UDfA12 reaction network

Table D.2

Updated fractionation reactions

D.2 Surface Chemistry

We follow the rate equation approach as described in Hasegawa et al. (1992) and Hasegawa & Herbst (1993) and account for competing surface processes (Garrod & Pauly 2011) and chemically inactive bulk ice. The species on the surface need to be mobile in order to scan the surface for a suitable reaction partner. The thermal hopping time scale to move between adjacent binding sites on the surface is given by Hasegawa et al. (1992) as (D.5)

with the characteristic vibration frequency for the adsorbed species (D.6)

where nsite is the surface density of binding sites and m the mass of the particle and kB is the Boltzmann constant. We assume nsite = 1.5 × 1015cm−2 (Tielens & Allamandola 1987). ED,i is the desorption energy for physical desorption, is the energy barrier between binding sites. We note that in this section all energies ED, EB, Ea are given in units of K. Following Hasegawa et al. (1992) we assume EB,i ≈ 0.3ED,i. For H and H2 tunneling might also be a more rapid migration process compared to thermal hopping with a tunneling time: (D.7)

where a = 2Å is the energy barrier width assuming a rectangular barrier (Garrod & Pauly 2011). We note that Eq. D.7 is not valid for particles heavier than H2. The time for an adsorbed particle to sweep over a number of binding sites is tdiff,i and the diffusion rate Rdiff,i is defined as the inverse of the diffusion time tdiff,i = NS × max(thop,i, ttunnel,i) [s]. NS is the total number of surface sites per grain. The surface reaction rate coefficient can then be written as: (D.8)

where ngrain is the number density of dust grains. The denominator gives the total number of surface binding sites. The probability that a reaction occurs is given by feff,ij and is assumed to be unity for an exothermic reaction without activation energy.19 For an exothermic reaction between surface species i and j with activation energy Ea,ij the probability of reaction during a single collision between the reactants can be expressed as a Boltzmann factor exp [−Ea,ij/Td] or using the quantum mechanical tunneling probability. Later probability is only relevant for light species and we apply it only if one of the reactants is J(H) or J(H2). (D.9)

where μij = mimj/(mi + mj) is the reduced mass. During the collision between i and j overcoming the activation energy barrier is in competition with the migration of either reactant to a neighbor site and with desorption of either reactant. We use the reaction probability feff,ij as described by Chang et al. (2007) and Garrod & Pauly (2011) to account for the competing processes in case an activation energy barrier exists (D.10)

with v0 = max(v0,i, v0,j) and the species dependent diffusion rate coefficient is a placeholder for species i and j). Kdes,_ is the sum of all desorption probabilities.

Accretion. The adsorption or accretion rate coefficient depends on the thermal velocity of the incoming particle vth,i, the total effective cross section area of all target grains σdust(cm2), that is the cross section per grain times the dust number density nd, and the sticking coefficient S (T, Td) as: (D.11)

Table D.3

Photo-desorption yields

with the mean thermal velocity (D.12)

For the sticking of H we use the coefficient provided by Hollenbach & McKee (1979, their Eq. (3.7)), in case of H2 we assume S (T, Td) = 0.5 (Acharyya 2014) and for any other species we assume S (T, Td) = 1/3 (Willacy & Williams 1993).

Thermal desorption. Species on the grain surface are removed depending on the grain surface and their binding energy which is itself a function of the grain material and ice composition. The rate coefficient is given by (D.13)

where v0,i is from Eq. (D.6).

H2 formation induced desorption. The H2 binding energy of 4.48 eV is released during the formation process. It is commonly assumed that the energy is distributed evenly between internal energy and kinetic energy of the desorbed molecule and lattice energy of the surface. This local “hotspot”-heating of the surface may be sufficient to desorb additional species. Following Willacy et al. (1994) we apply this only for species with a binding energy of less than 1210 K (Roberts et al. 2007). The rate coefficient is: (D.14)

where is the H2 formation rate on grains (Röllig et al. 2013) and e measures the number of atoms and molecules that are desorbed by this process. The value of e is uncertain and we assume a low efficiency of ϵ = 0.01. Division by the total density of all surface species ns(tot) = ∑i ns(i) distributes the binding energy across all species on the grain. Depending on the value of ϵ this process can dominate other desorption processes and we found it to significantly contribute to the overall desorption. For a more detailed analysis see Roberts et al. (2007) and more recently von Procházka & Millar (2021).

Photo-desorption. The absorption of UV photons at the dust surface can sufficiently increase the internal energy of the surface species to induce desorption. Following Cuppen et al. (2017) we write the rate coefficient for photo-desorption as: (D.15)

where Yi is the photo-desorption yield per UV photon. Experimental data on Yi is only available for very few species and we assume a general value of Yi = 10−3 (molecules/photon) except for the species listed in Table D.3. For H2O we use the expressions given by Öberg et al. (2009) accounting for ice thickness and gas and dust temperature (their Eq. 10). χFDraine is a measure of the local UV energy density with photons cm−2 s−1. fss is the self-shielding factor which is unity for all species except CO where the shielding of the ice species is also provided by the column in the gas. 〈exp (−γiAV)〉 accounts for the dust attenuation of FUV radiation from the cloud surface to the local position. The 〈〉 indicates an average over the full solid angle.

Total desorption. The total desorption rate coefficient is the sum of all individual desorption rate coefficients: (D.16)

where fdes describes the fraction of all surface species that are considered candidates for desorption: (D.17)

where Σnsite gives the total number of all grain surface binding sites per volume summed over all grains (see also Aikawa et al 1996; Woitke etal 2009) InEq (D 17) we make the assumption, that surface species from the top two layers of the ice mantle can desorb to the gas phase (Aikawa et al 1996) For the chemical desorption we assume fsurf = fdes, where fsurf is the fraction of particles that can undergo surface reactions We note that the effects of ice opacity to FUV radiation and the fact that photo-desorption only occurs from the top few monolayers of the ice is implicitly included in the photo-desorption yields (Cuppen et al 2017)

Table D.4 list all chemical desorption reactions with their respective branching ratios. ∆HR is the exothermicity of the reaction, BR is the branching ratio and Ea denotes the activation barrier of the reactions. The BR are computed per reaction (left-hand side). In case of multiple product channels for reaction J(A) + J(B) they have to be rescaled accordingly. We assumed the following desorption energies ED,i for the reaction products in the Table: H(600K), H2(430K), CO(1300K), N2(1100K), O2(1200K), O3(1800K), OH(4600K), CO2(2600K), H2O(5600K), HCO(2400K), O2H(3650K), CH3O(4400K), H2CO(4500K), H2O2(6000K), CH3OH(5000K). See Table 3 and Wakelam et al. (2015) (https://kida.astrochem-tools.org/ for references.)

Table D.4

Chemical desorption reactions.

Appendix E Heating and cooling processes

Table 5 lists all cooling and heating processes implemented in KOSMA-τ with their references. For those that are not copied from the literature we specify their details here. We chose the notation such that heating and cooling processes are denoted with ΓID and ΛID, respectively with ID taken from Table 5. Processes which can act as heating or cooling processes are named according to their dominating behavior.

H2 heating and cooling. Table 5 lists four heating and cooling processes related to H2 which are all based on the transfer of internal molecular energy to kinetic gas energy by means of inelastic collision.

UV radiation pumps H2 molecules to excited electronic states (Lyman and Werner bands) which decay quickly back to ground electronic states (except for η ≈ 10% of the excited states which decay into the continuum as pointed out by Solomon (1965) (cited in Field et al. 1966)). The decay leads to a population of superthermal states and collisional de-excitation leads to heating of the gas with the rate: (E.1)

where the sums over i and j are over all considered states. is the level population of the i-th level, C(i → j) is the total collision rate from level i to j and Ei − Ej is the energy difference between both levels (see also Sternberg & Dalgarno 1989). Depending on the level population this process can also cool the gas. Various approximations to this computation have been published, for example by Tielens & Hollenbach (1985) and Röllig et al. (2006).

Photo-dissociation of H2 molecules produces two energetic atoms which transfer their energy to the gas. The heating rate depends on the FUV pumping rate P, the dissociation probability η and the average kinetic energy 〈Ej*〉 ≈ 0.4eV of the atoms produced by the spontaneous decay out of the excited states of the Lyman and Werner electronic bands (Stephens & Dalgarno 1973): (E.2)

where i indicates all considered electronic ground state levels of H2 and j* indicates all electronically excited levels.

The formation of H2 molecules releases the binding energy of erg. Under the assumption of energy equipartition between grain lattice energy, inner energy and kinetic energy, 1/3 of this energy is converted to kinetic energy. With the H2-formation rate Rf(T) we find for the heating rate: (E.3)

Details on the H2 formation treatment in KOSMA-τ are given in Röllig et al. (2013).

In hot gas, H–H2 and H2–H2 collisions are energetic enough to dissociate molecular hydrogen. This removes the amount of from the kinetic gas energy (Lepp & Shull 1983). The cooling rate can be written as: (E.4)

where the rate coefficients for the collisional dissociation reactions H + H2 → H + H + H and H2 + H2 → H2 + H + H are taken from UDfA12: (E.5) (E.6)

Cosmic ray heating. The interaction of energetic cosmic ray protons with H and H2 produces hot electrons that transfer their energy to the ISM in the following collisions. A detailed treatment as been presented by Glassgold et al. (2012). Generally, the heating rate of cosmic rays can be written as: (E.7)

where ζCR is the CR ionization rate per H nucleus in s−1 and Q is the average energy deposited as heat per ionization. Glassgold et al. (2012) provide values of Q for different interstellar environments as function of total density (their Table 6). Q can then be approximated by the following expression: (E.8)

with x = log10(n/cm−3). The prefactor converts from eV to erg.

Grain photo-electric heating. The heating by hot electrons released during photo-electric (PE) absorption of FUV photons by dust particles is the most important PDR heating process under most circumstances. The photo-electrons can also recombine with the grains and cool the gas with a rate Λ8, therefore ГPE = Г8 − Λ8.

KOSMA-τ offers a choice of how to treat PE heating ГРЕ. The user can decide which approximation of the PE heating (and electron recombination) to use. We implemented the description given by Bakes & Tielens (1994, (their Eqs. (41),(43),(44))) and by Weingartner & Draine (2001c, (their Eqs. (44),(45) and Table 2 and 3)). They mainly differ in the applied dust size distributions and the estimates for photoelectric thresholds, yields, and electron capture rates. See Röllig et al. (2013) for more details.

Gas-grain collisions. Collision between gas particles and dust grains can transfer energy between the two. If the gas is hotter than the grains the collisions cool the gas. If the gas is cooler the collisions act as gas heating. Burke & Hollenbach (1983) give the rate as (E.9)

We apply an approximate accommodation factor when using the PE heating description by Bakes & Tielens (1994). In case of the PE heating described by Weingartner & Draine (2001a) we employ a modified accommodation factor given by Groenewegen (1994): (E.10)

Carbon ionization. The ionization of atomic carbon releases hot electrons with an energy of about 1 eV = 1.602 × 10−12erg. We approximate the heating rate by (E.11)

Line cooling. Details on the computation of atomic and molecular line cooling can also be found elsewhere (e.g., Tielens & Hollenbach 1985; Sternberg & Dalgarno 1989; Röllig et al. 2006; Woitke et al. 2009; Papadopoulos et al. 2011). Following Sternberg (1988) the general expression for the cooling rate due to a transition from level i to level j of species x is: (E.12)

where Aij and Bij are the Einstein coefficients, ni is the level population of level i, Eij is the external radiation at the transition frequency vij and β(τij) is the escape probability at the optical depth τ. For all cooling lines we compute τ from the perpendicular column density of the respective species. We use the expression for a spherical cloud given by Stutzki & Winnewisser (1985): (E.13)

More details are also given in Störzer et al. (1996). The energy level population is computed by solving the system of rate equations balancing all populating and de-populating processes. Using the escape probability approximation these equilibrium equations take a simple form. For the level i we can write: (E.14)

where ∑j>i sums over all levels with Ej < Ei. The first term describes all populating collisions, the second term all depopulating collisions and emission events. For two and three-level systems the solution to Eq. (E.14) can be given analytically. Larger n-level systems are solved by LU decomposition. Presently, KOSMA-τ does not consider pumping by external radiation and we set Eij = 0 in Eq. (E.12).

We also added two high-temperature cooling processes: O I 6300 Å emission of the meta-stable 1D level of atomic oxygen, and H i Lyman-α emission, which become relevant for T > 5000 K. We use the following expressions (Sternberg 1988): (E.15) (E.16)

Table E.1

Fine-structure and rotational line cooling data

Appendix F The clumpy KOSMA-τ setup

Stutzki et al. (1998) showed that the observed clump-mass spectrum of the ISM can be described as power-law spectrum: (F.1)

giving the number of clumps dNcl in the mass bin dMcl. In addition the masses of the clumps are related to their radii Rcl by the mass-size relation (F.2)

The power-law indices α and γ can be derived from observations. Kramer et al. (1998) present clump mass spectra of seven different molecular clouds and find that α = 1.6…1.8 for all clouds from their sample. As a consequence low-mass clumps are much more numerous compared to massive ones. Schneider & Brooks (2004) find similar results with a slightly larger variation of α = 1.6…2.1 and with a tentative trend of higher values of α in regions with higher gas density. Other studies (Elmegreen & Falgarone 1996) find 〈γ〉 = 2.2 − 2.5 over a large sample of clouds. Heithausen et al. (1998) studied the Polaris Flare to derive the power law spectrum over a mass range of at least five orders of magnitude, down to masses less than 10−3 M⊙. They find α = 1.84 and γ = 2.31 and the default setup in KOSMA-τ assumes these values. Assuming that the clump masses in an ensemble are all within the range ml ≤ Mc1 ≤ mu we can write the total number of clumps in an ensemble Nens as (F.3)

and the total ensemble mass (F.4)

Taking Mens as model parameter of the clumpy Eq. (F.4) gives the constant A. The constant C can be written as: (F.5)

where ρens is the averaged ensemble density.

We note that Ad scales linearly with Mens and therefore all ensemble quantities that result from the linear superposition of the individual clumps also scale with the ensemble mass. Eq. (F.2) leads to smaller clumps being denser while Eq. (F.1) returns many more low-mass clumps compared to high-mass clumps. To illustrate the consequences for the ensemble we summarize some ensemble properties in Table F.1.

Table F.1 shows that the clump number distribution is heavily dominated by the smallest clumps ( 84%) and the mass-size relation states that smaller, that is less massive clumps, have a higher density. The most mass however is concentrated in the higher mass clumps, for instance 62% of the total ensemble mass is in the three most massive clumps of the ensemble. The volume distribution is even more skewed to the massive clumps. Over 90% of the total ensemble volume is filled be the three most massive clumps and 40% of the projected area.

Within a single ensemble we assume here that the clumpy model ensemble is spatially unresolved and ignore mutual shielding and absorption within the ensemble. The beam (ΩB) averaged emission of the clumpy model is then just the weighted superposition of the individual clumps. The weighting factor is the total solid angle of all clumps in the mass bin dMcl with being the solid angle of a single clump at a distance D, assuming D >> R. (F.6)

The integral is taken over all masses in the range [m, mu]. Full radiative transfer is implemented in KOSMA-τ-3D (Andree-Labsch et al. 2017).

Table F.1

Properties of a discrete example ensemble

Appendix G WL-PDR - A sandbox PDR model code written in Wolfram Mathematica

It is useful to reduce the numerical complexity of a model in order to better study numerical or physical aspects of the model. For some numerical tests in this paper we use the simple 1-D, semi-infinite plane-parallel PDR model code WL-PDR. The code is ported to Mathematica from PyPDR (Bruderer 2019). PyPDR is used as an educational tool to learn about the basic physics and some numerical internals of PDR codes. We ported PyPDR to Mathematica by Wolfram Research (Wolfram Research Inc. 2020) to seamlessly integrate the toy PDR model into our research infrastructure and to make use of the advanced analysis, numerics and visualization capabilities of Mathemat-ica.

The features of WL-PDR are summarized in Table G.1. The code can be downloaded at https://github.com/markusroellig/WL-PDR.

Table G.1

WL-PDR features

References

  1. Abel, N. P. 2006, MNRAS, 368, 1949 [Google Scholar]
  2. Abel, N. P., & Ferland, G. J. 2006, ApJ, 647, 367 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abel, N. P., Ferland, G. J., Shaw, G., & van Hoof, P. A. M. 2005, ApJS, 161, 65 [NASA ADS] [CrossRef] [Google Scholar]
  4. Abgrall, H., Le Bourlot, J., Pineau Des Forets, G., et al. 1992, A&A, 253, 525 [Google Scholar]
  5. Abrahamsson, E., Krems, R. V., & Dalgarno, A. 2007, ApJ, 654, 1171 [NASA ADS] [CrossRef] [Google Scholar]
  6. Acharyya, K. 2014, MNRAS, 443, 1301 [Google Scholar]
  7. Agündez, M., Goicoechea, J. R., Cernicharo, J., Faure, A., & Roueff, E. 2010, ApJ, 713, 662 [CrossRef] [Google Scholar]
  8. Agündez, M., Roueff, E., Le Petit, F., & Le Bourlot, J. 2018, A&A, 616, A19 [CrossRef] [EDP Sciences] [Google Scholar]
  9. Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684 [NASA ADS] [CrossRef] [Google Scholar]
  10. Akimkin, V. V., Kirsanova, M. S., Pavlyuchenkov, Y. N., & Wiebe, D. S. 2015, MNRAS, 449, 440 [NASA ADS] [CrossRef] [Google Scholar]
  11. Akimkin, V. V., Kirsanova, M. S., Pavlyuchenkov, Y. N., & Wiebe, D. S. 2017, MNRAS, 469, 630 [NASA ADS] [Google Scholar]
  12. Allison, A. C., & Dalgarno, A. 1969, Atomic Data, 1, 289 [NASA ADS] [CrossRef] [Google Scholar]
  13. Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users' Guide, 3rd edn. (Philadelphia, PA: Society for Industrial and Applied Mathematics) [CrossRef] [Google Scholar]
  14. Andersson, S., & van Dishoeck, E. F. 2008, A&A, 491, 907 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Andree-Labsch, S., Ossenkopf-Okada, V., & Röllig, M. 2017, A&A, 598, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Andrews, H., Candian, A., & Tielens, A. G. G. M. 2016, A&A, 595, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Arkhipova, V. P., Egorov, O. V., Esipov, V. F., et al. 2013, MNRAS, 432, 2273 [NASA ADS] [CrossRef] [Google Scholar]
  18. Atkinson, R., Baulch, D. L., Cox, R. A., et al. 2006, Atmos. Chem. Phys., 6, 3625 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef] [Google Scholar]
  20. Ballesteros-Paredes, J., André, P., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 76 [Google Scholar]
  21. Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ApJ, 312, L45 [Google Scholar]
  22. Barinovs, G., van Hemert, M. C., Krems, R., & Dalgarno, A. 2005, ApJ, 620, 537 [NASA ADS] [CrossRef] [Google Scholar]
  23. Barzel, B., & Biham, O. 2007, ApJ, 658, L37 [NASA ADS] [CrossRef] [Google Scholar]
  24. Bates, D. R., & Spitzer Lyman, J. 1951, ApJ, 113, 441 [NASA ADS] [CrossRef] [Google Scholar]
  25. Baulch, D. L., Bowman, C. T., Cobos, C. J., et al. 2005, J. Phys. Chem. Ref. Data, 34, 757 [Google Scholar]
  26. Bell, T. A., Viti, S., Williams, D. A., Crawford, I. A., & Price, R. J. 2005, MNRAS, 357, 961 [NASA ADS] [CrossRef] [Google Scholar]
  27. Bell, T. A., Hartquist, T. W., Viti, S., & Williams, D. A. 2006, A&A, 459, 805 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Berné, O., Habart, É., Peeters, E., et al. 2022, PASP, 134, 1035 [Google Scholar]
  29. Bertoldi, F., & Draine, B. T. 1996, ApJ, 458, 222 [NASA ADS] [CrossRef] [Google Scholar]
  30. Bialy, S., & Sternberg, A. 2016, ApJ, 822, 83 [NASA ADS] [CrossRef] [Google Scholar]
  31. Bialy, S., Burkhart, B., & Sternberg, A. 2017, ApJ, 843, 92 [NASA ADS] [CrossRef] [Google Scholar]
  32. Bigiel, F., de Looze, I., Krabbe, A., et al. 2020, ApJ, 903, 30 [CrossRef] [Google Scholar]
  33. Bisbas, T. G., Wünsch, R., Whitworth, A. P., Hubber, D. A., & Walch, S. 2011, ApJ, 736, 142 [Google Scholar]
  34. Bisbas, T. G., Bell, T. A., Viti, S., Yates, J., & Barlow, M. J. 2012, MNRAS, 427, 2100 [NASA ADS] [CrossRef] [Google Scholar]
  35. Bisbas, T. G., Haworth, T. J., Barlow, M. J., et al. 2015, MNRAS, 454, 2828 [NASA ADS] [CrossRef] [Google Scholar]
  36. Bisbas, T. G., Tan, J. C., & Tanaka, K. E. I. 2021, MNRAS, 502, 2701 [CrossRef] [Google Scholar]
  37. Black, J.H., & Dalgarno, A. 1977, ApJS, 34, 405 [NASA ADS] [CrossRef] [Google Scholar]
  38. Black, J. H., & van Dishoeck, E. F. 1987, ApJ, 322, 412 [Google Scholar]
  39. Boger, G. I., & Sternberg, A. 2006, ApJ, 645, 314 [NASA ADS] [CrossRef] [Google Scholar]
  40. Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles, eds. C. F. Bohren, & D. R. Huffman [Google Scholar]
  41. Bolatto, A. D., Jackson, J. M., & Ingalls, J. G. 1999, ApJ, 513, 275 [NASA ADS] [CrossRef] [Google Scholar]
  42. Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541 [Google Scholar]
  43. Boschman, L., Cazaux, S., Spaans, M., Hoekstra, R., & Schlathölter, T. 2015, A&A, 579, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Bottinelli, S., Boogert, A. C. A., Bouwman, J., et al. 2010, ApJ, 718, 1100 [NASA ADS] [CrossRef] [Google Scholar]
  45. Bron, E., Le Bourlot, J., & Le Petit, F. 2014, A&A, 569, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Bron, E., Agündez, M., Goicoechea, J. R., & Cernicharo, J. 2018, A&A submitted, [arXiv: 1881.81547] [Google Scholar]
  47. Bruderer, S. 2013, A&A, 559, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Bruderer, S. 2019, PyPDR: Python Photo Dissociation Regions Astrophysics Source Code Library, [record ascl:1985.827] [Google Scholar]
  49. Bruderer, S., Benz, A. O., Doty, S. D., van Dishoeck, E. F., & Bourke, T. L. 2009a, ApJ, 700, 872 [NASA ADS] [CrossRef] [Google Scholar]
  50. Bruderer, S., Doty, S. D., & Benz, A. O. 2009b, ApJS, 183, 179 [Google Scholar]
  51. Bruderer, S., Benz, A. O., Stäuber, P., & Doty, S. D. 2010, ApJ, 720, 1432 [NASA ADS] [CrossRef] [Google Scholar]
  52. Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223 [NASA ADS] [CrossRef] [Google Scholar]
  54. Burton, M. G., Hollenbach, D. J., & Tielens, A. G. G. M. 1990, ApJ, 365, 620 [NASA ADS] [CrossRef] [Google Scholar]
  55. Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29 [NASA ADS] [CrossRef] [Google Scholar]
  56. Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222 [NASA ADS] [CrossRef] [Google Scholar]
  57. Cazaux, S., Minissale, M., Dulieu, F., & Hocuk, S. 2016, A&A, 585, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Chabot, M., Béroff, K., Gratier, P., Jallat, A., & Wakelam, V. 2013, ApJ, 771, 90 [NASA ADS] [CrossRef] [Google Scholar]
  59. Chang, Q., Cuppen, H. M., & Herbst, E. 2007, A&A, 469, 973 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Charnley, S. B., & Markwick, A. J. 2003, A&A, 399, 583 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Cheng, Y., Tan, J. C., Liu, M., et al. 2018, ApJ, 853, 160 [NASA ADS] [CrossRef] [Google Scholar]
  62. Choi, Y.-J., Min, K.-W., Seon, K.-I., et al. 2013, ApJ, 774, 34 [NASA ADS] [CrossRef] [Google Scholar]
  63. Choi, Y.-J., Min, K.-W., & Seon, K.-I. 2015, ApJ, 800, 132 [NASA ADS] [CrossRef] [Google Scholar]
  64. Clavel, J., Viala, Y. P., & Bel, N. 1978, A&A, 65, 435 [NASA ADS] [Google Scholar]
  65. Cook, G. R., Metzger, P. H., & Ogawa, M. 1968, J. Opt. Soc. Am., 58, 129 [NASA ADS] [CrossRef] [Google Scholar]
  66. Cormier, D., Abel, N. P., Hony, S., et al. 2019, A&A, 626, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Crawford, M. K., Genzel, R., Townes, C. H., & Watson, D. M. 1985, ApJ, 291, 755 [Google Scholar]
  68. Cubick, M., Stutzki, J., Ossenkopf, V., Kramer, C., & Röllig, M. 2008, A&A, 488, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Cuppen, H. M., & Herbst, E. 2007, ApJ, 668, 294 [Google Scholar]
  70. Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1 [Google Scholar]
  71. Dalgarno, A., & Stephens, T. L. 1970, ApJ, 160, L107 [NASA ADS] [CrossRef] [Google Scholar]
  72. Das, A., Sil, M., Ghosh, R., et al. 2021, Front. Astron. Space Sci., 8, 78 [NASA ADS] [CrossRef] [Google Scholar]
  73. de Jong, T., Boland, W., & Dalgarno, A. 1980, A&A, 91, 68 [NASA ADS] [Google Scholar]
  74. Decataldo, D., Ferrara, A., Pallottini, A., Gallerani, S., & Vallini, L. 2017, MNRAS, 471, 4476 [NASA ADS] [CrossRef] [Google Scholar]
  75. Decataldo, D., Pallottini, A., Ferrara, A., Vallini, L., & Gallerani, S. 2019, MNRAS, 487, 3377 [Google Scholar]
  76. Dedes, C., Röllig, M., Mookerjea, B., et al. 2010, A&A, 521, L24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Dewangan, D. P., Flower, D. R., & Alexander, M. H. 1987, MNRAS, 226, 505 [NASA ADS] [Google Scholar]
  78. Draine, B. T. 1978, ApJS, 36, 595 [Google Scholar]
  79. Draine, B. T. 2003, ApJ, 598, 1017 [NASA ADS] [CrossRef] [Google Scholar]
  80. Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269 [Google Scholar]
  81. Dufour, G., & Charnley, S. B. 2019, ApJ, 887, 67 [NASA ADS] [CrossRef] [Google Scholar]
  82. Dulieu, F., Congiu, E., Noble, J., et al. 2013, Scientific Rep., 3, 1338 [Google Scholar]
  83. Elmegreen, B. G., & Falgarone, E. 1996, ApJ, 471, 816 [NASA ADS] [CrossRef] [Google Scholar]
  84. Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, J. Mol. Spectr., 327, 95 [NASA ADS] [CrossRef] [Google Scholar]
  85. Esplugues, G. B., Cazaux, S., Meijerink, R., Spaans, M., & Caselli, P. 2016, A&A, 591, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Esplugues, G., Cazaux, S., Caselli, P., Hocuk, S., & Spaans, M. 2019, MNRAS, 486, 1853 [NASA ADS] [Google Scholar]
  87. Fayolle, E. C., Bertin, M., Romanzin, C., et al. 2011, ApJ, 739, L36 [Google Scholar]
  88. Federman, S. R., Glassgold, A. E., & Kwan, J. 1979, ApJ, 227, 466 [NASA ADS] [CrossRef] [Google Scholar]
  89. Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137 [Google Scholar]
  90. Ferland, G. J., Chatzikos, M., Guzmân, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385 [Google Scholar]
  91. Field, G. B., Somerville, W. B., & Dressler, K. 1966, ARA&A, 4, 207 [NASA ADS] [CrossRef] [Google Scholar]
  92. Foley, N., Cazaux, S., Egorov, D., et al. 2018, MNRAS, 479, 649 [NASA ADS] [Google Scholar]
  93. Franeck, A., Walch, S., Seifried, D., et al. 2018, MNRAS, 481, 4277 [NASA ADS] [CrossRef] [Google Scholar]
  94. Fuente, A., Cernicharo, J., Roueff, E., et al. 2016, A&A, 593, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Gaches, B. A. L., Offner, S. S. R., & Bisbas, T. G. 2019, ApJ, 878, 105 [NASA ADS] [CrossRef] [Google Scholar]
  96. Gammie, C. F., Lin, Y.-T., Stone, J. M., & Ostriker, E. C. 2003, ApJ, 592, 203 [NASA ADS] [CrossRef] [Google Scholar]
  97. Garcia, R., Abel, N., Röllig, M., Simon, R., & Stutzki, J. 2021, A&A, 650, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485 [Google Scholar]
  99. Garrod, R. T. 2013, ApJ, 765, 60 [Google Scholar]
  100. Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15 [Google Scholar]
  102. Garrod, R. T., Weaver, S. L. W., & Herbst, E. 2008, ApJ, 682, 283 [NASA ADS] [CrossRef] [Google Scholar]
  103. Geppert, W. D., Hellberg, F., Österdahl, F., et al. 2005, in Astrochemistry: Recent Successes and Current Challenges, 231, eds. D. C. Lis, G. A. Blake, & E. Herbst, 117 [NASA ADS] [Google Scholar]
  104. Gierens, K. M., Stutzki, J., & Winnewisser, G. 1992, A&A, 259, 271 [NASA ADS] [Google Scholar]
  105. Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19 [NASA ADS] [CrossRef] [Google Scholar]
  106. Glassgold, A. E., & Langer, W. D. 1974, ApJ, 193, 73 [NASA ADS] [CrossRef] [Google Scholar]
  107. Glassgold, A. E., & Langer, W. D. 1975, ApJ, 197, 347 [NASA ADS] [CrossRef] [Google Scholar]
  108. Glassgold, A. E., & Langer, W. D. 1976, ApJ, 206, 85 [NASA ADS] [CrossRef] [Google Scholar]
  109. Glassgold, A. E., Galli, D., & Padovani, M. 2012, ApJ, 756, 157 [Google Scholar]
  110. Glover, S. C. O., Federrath, C., Mac Low, M. M., & Klessen, R. S. 2010, MNRAS, 404, 2 [Google Scholar]
  111. Goicoechea, J. R., & Le Bourlot, J. 2007, A&A, 467, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Goicoechea, J. R., Pety, J., Cuadrado, S., et al. 2016, Nature, 537, 207 [Google Scholar]
  113. Goicoechea, J. R., Cuadrado, S., Pety, J., et al. 2017, A&A, 601, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Goicoechea, J. R., Pabst, C. H. M., Kabanovic, S., et al. 2020, A&A, 639, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Goicoechea, J. R., Aguado, A., Cuadrado, S., et al. 2021, A&A, 647, A10 [EDP Sciences] [Google Scholar]
  116. Goldsmith, P. F., Liseau, R., Bell, T. A., et al. 2011, in The Molecular Universe, 280, eds. J. Cernicharo, & R. Bachiller, 33 [NASA ADS] [Google Scholar]
  117. Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13 [Google Scholar]
  118. Gong, H., & Ostriker, E. C. 2011, ApJ, 729, 120 [NASA ADS] [CrossRef] [Google Scholar]
  119. Gong, M., & Ostriker, E. C. 2015, ApJ, 806, 31 [NASA ADS] [CrossRef] [Google Scholar]
  120. Gong, M., Ostriker, E. C., & Wolfire, M. G. 2017, ApJ, 843, 38 [NASA ADS] [CrossRef] [Google Scholar]
  121. Gorti, U., & Hollenbach, D. 2002, ApJ, 573, 215 [Google Scholar]
  122. Gorti, U., & Hollenbach, D. 2004, ApJ, 613, 424 [NASA ADS] [CrossRef] [Google Scholar]
  123. Gorti, U., & Hollenbach, D. 2008, ApJ, 683, 287 [NASA ADS] [CrossRef] [Google Scholar]
  124. Gould, R. J., & Salpeter, E. E. 1963, ApJ, 138, 393 [NASA ADS] [CrossRef] [Google Scholar]
  125. Graf, U. U., Simon, R., Stutzki, J., et al. 2012, A&A, 542, A16 [Google Scholar]
  126. Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [Google Scholar]
  127. Groenewegen, M. A. T. 1994, A&A, 290, 531 [NASA ADS] [Google Scholar]
  128. Grosch, H., Fateev, A., & Clausen, S. 2015, J. Quant. Spectr. Radiat. Transf., 154, 28 [NASA ADS] [CrossRef] [Google Scholar]
  129. Guszejnov, D., Hopkins, P. F., & Grudic, M. Y. 2018, MNRAS, 477, 5139 [NASA ADS] [CrossRef] [Google Scholar]
  130. Guzmân, V., Pety, J., Goicoechea, J. R., Gerin, M., & Roueff, E. 2011, A&A, 534, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  131. Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421 [Google Scholar]
  132. Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] [CrossRef] [Google Scholar]
  133. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [Google Scholar]
  134. Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Hebrard, E., Dobrijevic, M., Pernot, P., et al. 2009, J. Phys. Chem. A, 113, 11227 [NASA ADS] [CrossRef] [Google Scholar]
  136. Heithausen, A., Bensch, F., Stutzki, J., Falgarone, E., & Panis, J. F. 1998, A&A, 331, L65 [Google Scholar]
  137. Henney, W. J., Arthur, S. J., de Colle, F., & Mellema, G. 2009, MNRAS, 398, 157 [NASA ADS] [CrossRef] [Google Scholar]
  138. Herbst, E., & van Dishoeck, E. F. 2009, ARA&A, 47, 427 [NASA ADS] [CrossRef] [Google Scholar]
  139. Heyminck, S., Graf, U. U., Güsten, R., et al. 2012, A&A, 542, A1 [Google Scholar]
  140. Hincelin, U., Chang, Q., & Herbst, E. 2015, A&A, 574, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Hindmarsh, A. C. 1983, Scientific Computing, 55 [Google Scholar]
  142. Hocuk, S., Szucs, L., Caselli, P., et al. 2017, A&A, 604, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  143. Hogerheijde, M. R., Jansen, D. J., & van Dishoeck, E. F. 1995, A&A, 294, 792 [NASA ADS] [Google Scholar]
  144. Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 [Google Scholar]
  145. Hollenbach, D., & Salpeter, E. E. 1971, ApJ, 163, 155 [NASA ADS] [CrossRef] [Google Scholar]
  146. Hollenbach, D. J., & Tielens, A. G. G. M. 1997, ARA&A, 35, 179 [NASA ADS] [CrossRef] [Google Scholar]
  147. Hollenbach, D. J., Takahashi, T., & Tielens, A. G. G. M. 1991, ApJ, 377, 192 [NASA ADS] [CrossRef] [Google Scholar]
  148. Hollenbach, D., Kaufman, M. J., Bergin, E. A., & Melnick, G. J. 2009, ApJ, 690, 1497 [NASA ADS] [CrossRef] [Google Scholar]
  149. Inoue, S., Yoshida, N., & Yajima, H. 2020, MNRAS, 498, 5960 [NASA ADS] [CrossRef] [Google Scholar]
  150. Izumi, N., Fukui, Y., Tachihara, K., et al. 2021, PASJ, 73, 174 [NASA ADS] [CrossRef] [Google Scholar]
  151. Jo, Y.-S., Min, K.-W., Seon, K.-I., Edelstein, J., & Han, W. 2011, ApJ, 738, 91 [NASA ADS] [CrossRef] [Google Scholar]
  152. Jo, Y.-S., Min, K.-W., Lim, T.-H., & Seon, K.-I. 2012, ApJ, 756, 38 [NASA ADS] [CrossRef] [Google Scholar]
  153. Joblin, C., Bron, E., Pinto, C., et al. 2018, A&A, 615, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  154. Jonkheid, B., Faas, F. G. A., van Zadelhoff, G.-J., & van Dishoeck, E. F. 2004, A&A, 428, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  155. Jura, M. 1974, ApJ, 191, 375 [NASA ADS] [CrossRef] [Google Scholar]
  156. Kabanovic, S., Schneider, N., Ossenkopf-Okada, V., et al. 2022, A&A, 659, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  157. Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  158. Kamp, I., Thi, W. F., Woitke, P., et al. 2017, A&A, 607, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  159. Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795 [Google Scholar]
  160. Kaufman, M. J., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 644, 283 [Google Scholar]
  161. Keller-Rudek, H., Moortgat, G. K., Sander, R., & Sörensen, R. 2013, Earth Syst. Sci. Data, 5, 365 [NASA ADS] [CrossRef] [Google Scholar]
  162. Kirsanova, M. S., Wiebe, D. S., & Sobolev, A. M. 2009, Astron. Rep., 53, 611 [NASA ADS] [CrossRef] [Google Scholar]
  163. Kirsanova, M. S., Pavlyuchenkov, Y. N., Wiebe, D. S., et al. 2019, MNRAS, 488, 5641 [Google Scholar]
  164. Kirsanova, M. S., Ossenkopf-Okada, V., Anderson, L. D., et al. 2020, MNRAS, 497, 2651 [NASA ADS] [CrossRef] [Google Scholar]
  165. Kong, S. 2019, ApJ, 873, 31 [NASA ADS] [CrossRef] [Google Scholar]
  166. Könyves, V., André, P., Arzoumanian, D., et al. 2020, A&A, 635, A34 [Google Scholar]
  167. Kramer, C., Stutzki, J., Rohrig, R., & Corneliussen, U. 1998, A&A, 329, 249 [NASA ADS] [Google Scholar]
  168. Kramer, C., Cubick, M., Röllig, M., et al. 2008, A&A, 477, 547 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  169. Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2008, ApJ, 689, 865 [NASA ADS] [CrossRef] [Google Scholar]
  170. Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 693, 216 [NASA ADS] [CrossRef] [Google Scholar]
  171. Langer, W. 1976, ApJ, 206, 699 [NASA ADS] [CrossRef] [Google Scholar]
  172. Launay, J. M., & Roueff, E. 1977, A&A, 56, 289 [NASA ADS] [Google Scholar]
  173. Le Bourlot, J., Pineau Des Forets, G., Roueff, E., & Flower, D. R. 1993a, A&A, 267, 233 [NASA ADS] [Google Scholar]
  174. Le Bourlot, J., Pineau des Forets, G., Roueff, E., & Schilke, P. 1993b, ApJ, 416, L87 [NASA ADS] [CrossRef] [Google Scholar]
  175. Le Bourlot, J., Pineau des Forets, G., Roueff, E., & Flower, D. R. 1995, A&A, 302, 870 [Google Scholar]
  176. Le Bourlot, J., Le Petit, F., Pinto, C., Roueff, E., & Roy, F. 2012, A&A, 541, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  177. Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef] [Google Scholar]
  178. Lee, H. H., Roueff, E., Pineau des Forets, G., et al. 1998, A&A, 334, 1047 [Google Scholar]
  179. Lefloch, B., & Lazareff, B. 1994, A&A, 289, 559 [NASA ADS] [Google Scholar]
  180. Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147 [Google Scholar]
  181. Lepp, S., & Shull, J. M. 1983, ApJ, 270, 578 [NASA ADS] [CrossRef] [Google Scholar]
  182. Levrier, F., Le Petit, F., Hennebelle, P., et al. 2012, A&A, 544, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  183. Li, A., & Draine, B. T. 2001, ApJ, 554, 778 [Google Scholar]
  184. Li, Q., Tan, J. C., Christie, D., Bisbas, T. G., & Wu, B. 2018, PASJ, 70, S56 [NASA ADS] [Google Scholar]
  185. Li, C., Wang, H.-C., Wu, Y.-W., Ma, Y.-H., & Lin, L.-H. 2020, Res. Astron. Astrophys., 20, 031 [CrossRef] [Google Scholar]
  186. Lim, T.-H., Min, K.-W., & Seon, K.-I. 2013, ApJ, 765, 107 [NASA ADS] [CrossRef] [Google Scholar]
  187. Lique, F., Werfelli, G., Halvick, P., et al. 2013, J. Chem. Phys., 138, 204314 [NASA ADS] [CrossRef] [Google Scholar]
  188. Luisi, M., Anderson, L. D., Schneider, N., et al. 2021, Sci. Adv., 7, eabe9511 [NASA ADS] [CrossRef] [Google Scholar]
  189. Maddalena, R. J., Morris, M., Moscowitz, J., & Thaddeus, P. 1986, ApJ, 303, 375 [NASA ADS] [CrossRef] [Google Scholar]
  190. Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  191. Maillard, V., Bron, E., & Le Petit, F. 2021, A&A, 656, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  192. Makai, Z. S. 2015, PhD thesis, BCGS, University of Cologne, Germany [Google Scholar]
  193. Marconi, A., Testi, L., Natta, A., & Walmsley, C. M. 1998, A&A, 330, 696 [NASA ADS] [Google Scholar]
  194. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  195. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  196. McKee, C. F., & Krumholz, M. R. 2010, ApJ, 709, 308 [NASA ADS] [CrossRef] [Google Scholar]
  197. Meijerink, R., & Spaans, M. 2005, A&A, 436, 397 [CrossRef] [EDP Sciences] [Google Scholar]
  198. Meixner, M., & Tielens, A. G. G. M. 1993, ApJ, 405, 216 [NASA ADS] [CrossRef] [Google Scholar]
  199. Melnick, G., Gull, G. E., & Harwit, M. 1979, ApJ, 227, L29 [NASA ADS] [CrossRef] [Google Scholar]
  200. Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  201. Miville-Deschênes, M. A., Salome, Q., Martin, P. G., et al. 2017, A&A, 599, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  202. Mladenovic, M., & Roueff, E. 2014, A&A, 566, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  203. Mookerjea, B., Kramer, C., Röllig, M., & Masur, M. 2006, A&A, 456, 235 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  204. Mookerjea, B., Ossenkopf, V., Ricken, O., et al. 2012, A&A, 542, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  205. Mookerjea, B., Sandell, G., Güsten, R., et al. 2019, A&A, 626, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  206. Mookerjea, B., Sandell, G., Veena, V. S., et al. 2021, A&A, 648, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  207. Nagahama, T., Mizuno, A., Ogawa, H., & Fukui, Y. 1998, AJ, 116, 336 [NASA ADS] [CrossRef] [Google Scholar]
  208. Nagy, Z., Van der Tak, F. F. S., Ossenkopf, V., et al. 2013, A&A, 550, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  209. Nagy, Z., Choi, Y., Ossenkopf-Okada, V., et al. 2017, A&A, 599, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  210. Nayak, O., Meixner, M., Okada, Y., et al. 2021, ApJ, 907, 106 [NASA ADS] [CrossRef] [Google Scholar]
  211. Nejad, L. 2005, Astrophys. Space Sci., 299, 1 [NASA ADS] [CrossRef] [Google Scholar]
  212. Nelson, R. P., & Langer, W. D. 1997, ApJ, 482, 796 [NASA ADS] [CrossRef] [Google Scholar]
  213. Neufeld, D. A., & Melnick, G. J. 1987, ApJ, 322, 266 [NASA ADS] [CrossRef] [Google Scholar]
  214. Neufeld, D. A., Wolfire, M. G., & Schilke, P. 2005, ApJ, 628, 260 [Google Scholar]
  215. Öberg, K. I., Linnartz, H., Visser, R., & van Dishoeck, E. F. 2009, ApJ, 693, 1209 [Google Scholar]
  216. Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2011, ApJ, 740, 109 [Google Scholar]
  217. Okada, Y., Pilleri, P., Berné, O., et al. 2013, A&A, 553, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  218. Okada, Y., Güsten, R., Requena-Torres, M. A., et al. 2019a, A&A, 621, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  219. Okada, Y., Higgins, R., Ossenkopf-Okada, V., et al. 2019b, A&A, 631, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  220. Ossenkopf, V., Koumpia, E., Okada, Y., et al. 2015, A&A, 580, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  221. Pabst, C., Higgins, R., Goicoechea, J. R., et al. 2019, Nature, 565, 618 [NASA ADS] [CrossRef] [Google Scholar]
  222. Papadopoulos, P. P., Thi, W.-F., Miniati, F., & Viti, S. 2011, MNRAS, 414, 1705 [NASA ADS] [CrossRef] [Google Scholar]
  223. Parshley, S. C., Kronshage, J., Blair, J., et al. 2018, SPIE Conf. Ser., 10700, 107005X [NASA ADS] [Google Scholar]
  224. Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71 [Google Scholar]
  225. Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998, J. Quant. Spec. Radiat. Transf., 60, 883 [Google Scholar]
  226. Pineda, J. L., Mizuno, N., Stutzki, J., et al. 2008, A&A, 482, 197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  227. Pineda, J. E., Rosolowsky, E. W., & Goodman, A. A. 2009, ApJ, 699, L134 [NASA ADS] [CrossRef] [Google Scholar]
  228. Pound, M. W., & Wolfire, M. G. 2008, in Astronomical Society of the Pacific Conference Series, 394, Astronomical Data Analysis Software and Systems XVII, eds. R. W. Argyle, P. S. Bunclark, & J. R. Lewis, 654 [NASA ADS] [Google Scholar]
  229. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (New York, NY, USA: Cambridge University Press) [Google Scholar]
  230. Priestley, F. D., Barlow, M. J., & Viti, S. 2017, MNRAS, 472, 4444 [NASA ADS] [CrossRef] [Google Scholar]
  231. Putaud, T., Michaut, X., Le Petit, F., Roueff, E., & Lis, D. C. 2019, A&A, 632, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  232. Risacher, C., Güsten, R., Stutzki, J., et al. 2018, J. Astron. Instrum., 7, 1840014 [NASA ADS] [CrossRef] [Google Scholar]
  233. Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16 [Google Scholar]
  234. Roberts, J. F., Rawlings, J. M. C., Viti, S., & Williams, D. A. 2007, MNRAS, 382, 733 [NASA ADS] [CrossRef] [Google Scholar]
  235. Röllig, M. 2011, A&A, 530, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  236. Röllig, M., & Ossenkopf, V. 2013, A&A, 550, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  237. Röllig, M., Ossenkopf, V., Jeyakumar, S., Stutzki, J., & Sternberg, A. 2006, A&A, 451, 917 [Google Scholar]
  238. Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 [Google Scholar]
  239. Röllig, M., Kramer, C., Rajbahak, C., et al. 2011, A&A, 525, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  240. Röllig, M., Simon, R., Güsten, R., et al. 2012, A&A, 542, L22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  241. Röllig, M., Szczerba, R., Ossenkopf, V., & Glück, C. 2013, A&A, 549, A85 [Google Scholar]
  242. Röllig, M., Simon, R., Güsten, R., et al. 2016, A&A, 591, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  243. Ruaud, M., & Gorti, U. 2019, ApJ, 885, 146 [Google Scholar]
  244. Roueff, E., & Le Bourlot, J. 2020, A&A, 643, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  245. Ruaud, M., Loison, J. C., Hickson, K. M., et al. 2015, MNRAS, 447, 4004 [Google Scholar]
  246. Russell, R. W., Melnick, G., Gull, G. E., & Harwit, M. 1980, ApJ, 240, L99 [NASA ADS] [CrossRef] [Google Scholar]
  247. Russell, R. W., Melnick, G., Smyers, S. D., et al. 1981, ApJ, 250, L35 [NASA ADS] [CrossRef] [Google Scholar]
  248. Sandell, G., Mookerjea, B., Güsten, R., et al. 2015, A&A, 578, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  249. Sandford, S. A., & Allamandola, L. J. 1993, ApJ, 409, L65 [NASA ADS] [CrossRef] [Google Scholar]
  250. Schneider, N., & Brooks, K. 2004, PASA, 21, 290 [NASA ADS] [CrossRef] [Google Scholar]
  251. Schneider, N., Bontemps, S., Motte, F., et al. 2016, A&A, 591, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  252. Schneider, N., Röllig, M., Simon, R., et al. 2018, A&A, 617, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  253. Schneider, N., Röllig, M., Polehampton, E. T., et al. 2021, A&A, 653, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  254. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [Google Scholar]
  255. Schroder, K., Staemmler, V., Smith, M. D., Flower, D. R., & Jaquet, R. 1991, J. Phys. B: At. Mol. Opt. Phys., 24, 2487 [CrossRef] [Google Scholar]
  256. Schulz, A., Henkel, C., Muders, D., et al. 2007, A&A, 466, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  257. Seifried, D., & Walch, S. 2016, MNRAS, 459, L11 [NASA ADS] [CrossRef] [Google Scholar]
  258. Seifried, D., Walch, S., Girichidis, P., et al. 2017, MNRAS, 472, 4797 [NASA ADS] [CrossRef] [Google Scholar]
  259. Seon, K.-I., & Witt, A. N. 2012, ApJ, 758, 109 [NASA ADS] [CrossRef] [Google Scholar]
  260. Shaw, G., Ferland, G. J., Abel, N. P., Stancil, P. C., & van Hoof, P. A. M. 2005, ApJ, 624, 794 [NASA ADS] [CrossRef] [Google Scholar]
  261. Shimonishi, T., Nakatani, N., Furuya, K., & Hama, T. 2018, ApJ, 855, 27 [Google Scholar]
  262. Simön-Diaz, S. & Stasmska, G. 2011, A&A, 526, A48 [CrossRef] [EDP Sciences] [Google Scholar]
  263. Sipilä, O., Silsbee, K., & Caselli, P. 2021, ApJ, 922, 126 [CrossRef] [Google Scholar]
  264. Spaans, M., & van Dishoeck, E. F. 1997, A&A, 323, 953 [NASA ADS] [Google Scholar]
  265. Spitzer, L. 1978, Physical processes in the interstellar medium (New York: Wiley Interscience) [Google Scholar]
  266. Stacey, G. J., Smyers, S. D., Kurtz, N. T., & Harwit, M. 1983, ApJ, 265, L7 [NASA ADS] [CrossRef] [Google Scholar]
  267. Sternberg, A. 1988, ApJ, 332, 400 [NASA ADS] [CrossRef] [Google Scholar]
  268. Stephens, T. L., & Dalgarno, A. 1972, J. Quant. Spec. Radiat. Transf., 12, 569 [NASA ADS] [CrossRef] [Google Scholar]
  269. Stephens, T. L., & Dalgarno, A. 1973, ApJ, 186, 165 [NASA ADS] [CrossRef] [Google Scholar]
  270. Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197 [Google Scholar]
  271. Sternberg, A., & Dalgarno, A. 1995, ApJS, 99, 565 [NASA ADS] [CrossRef] [Google Scholar]
  272. Sternberg, A., Le Petit, F., Roueff, E., & Le Bourlot, J. 2014, ApJ, 790, 10 [CrossRef] [Google Scholar]
  273. Sternberg, A., Gurman, A., & Bialy, S. 2021, ApJ, 920, 83 [NASA ADS] [CrossRef] [Google Scholar]
  274. Störzer, H., & Hollenbach, D. 1998, ApJ, 495, 853 [CrossRef] [Google Scholar]
  275. Störzer, H., & Hollenbach, D. 1999, ApJ, 515, 669 [CrossRef] [Google Scholar]
  276. Störzer, H., Stutzki, J., & Sternberg, A. 1996, A&A, 310, 592 [NASA ADS] [Google Scholar]
  277. Störzer, H., Stutzki, J., & Sternberg, A. 1997, A&A, 323, L13 [NASA ADS] [Google Scholar]
  278. Störzer, H., Zielinsky, M., Stutzki, J., & Sternberg, A. 2000, A&A, 358, 682 [NASA ADS] [Google Scholar]
  279. Strassen, V. 1969, Numer. Math., 13, 354 [CrossRef] [Google Scholar]
  280. Stutzki, J., & Winnewisser, G. 1985, A&A, 144, 13 [NASA ADS] [Google Scholar]
  281. Stutzki, J., Stacey, G. J., Genzel, R., et al. 1988, ApJ, 332, 379 [NASA ADS] [CrossRef] [Google Scholar]
  282. Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697 [NASA ADS] [Google Scholar]
  283. Sun, K., Ossenkopf, V., Kramer, C., et al. 2008, A&A, 489, 207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  284. Szczerba, R., Omont, A., Volk, K., Cox, P., & Kwok, S. 1997, A&A, 317, 859 [NASA ADS] [Google Scholar]
  285. Tarantino, E., Bolatto, A. D., Herrera-Camus, R., et al. 2021, ApJ, 915, 92 [NASA ADS] [CrossRef] [Google Scholar]
  286. Thi, W. F., Lesur, G., Woitke, P., et al. 2019, A&A, 632, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  287. Thi, W. F., Hocuk, S., Kamp, I., et al. 2020, A&A, 634, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  288. Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge University Press) [Google Scholar]
  289. Tielens, A. G. G. M., & Allamandola, L. J. 1987, Composition, Structure, and Chemistry of Interstellar Dust, eds. D. J. Hollenbach, & J. Thronson, Harley, A., 134, 397 [NASA ADS] [Google Scholar]
  290. Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245 [Google Scholar]
  291. Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722 [Google Scholar]
  292. Valdivia, V., Hennebelle, P., Gérin, M., & Lesaffre, P. 2016, A&A, 587, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  293. van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15 [Google Scholar]
  294. van Dishoeck, E. F. 1987, in IAU Symposium, 120, Astrochemistry, eds. M. S. Vardya, & S. P. Tarafdar, 51 [NASA ADS] [CrossRef] [Google Scholar]
  295. van Dishoeck, E. F., & Black, J. H. 1986, ApJS, 62, 109 [Google Scholar]
  296. van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771 [Google Scholar]
  297. van Hoof, P. A. M., Weingartner, J. C., Martin, P. G., Volk, K., & Ferland, G. J. 2004, MNRAS, 350, 1330 [NASA ADS] [CrossRef] [Google Scholar]
  298. Viallet, M., Goffrey, T., Baraffe, I., et al. 2016, A&A, 586, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  299. Vidali, G., Ihm, G., Kim, H.-Y., & Cole, M. W. 1991, Surf. Sci. Rep., 12, 135 [NASA ADS] [CrossRef] [Google Scholar]
  300. Vieira, D., & Krems, R. V. 2017, ApJ, 835, 255 [NASA ADS] [CrossRef] [Google Scholar]
  301. Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  302. von Prochâzka, A. A. & Millar, T. J. 2021, MNRAS, 501, 1228 [Google Scholar]
  303. Wakelam, V., Loison, J. C., Herbst, E., et al. 2015, ApJS, 217, 20 [NASA ADS] [CrossRef] [Google Scholar]
  304. Wakelam, V., Bron, E., Cazaux, S., et al. 2017a, Mol. Astrophys., 9, 1 [Google Scholar]
  305. Wakelam, V., Loison, J. C., Mereau, R., & Ruaud, M. 2017b, Mol. Astrophys., 6, 22 [NASA ADS] [CrossRef] [Google Scholar]
  306. Wakelam, V., Dartois, E., Chabot, M., et al. 2021, A&A, 652, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  307. Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238 [Google Scholar]
  308. Weingartner, J., & Draine, B. 2001a, ApJS, 134, 263 [CrossRef] [Google Scholar]
  309. Weingartner, J. C., & Draine, B. T. 2001b, ApJ, 548, 296 [Google Scholar]
  310. Weingartner, J. C., & Draine, B. T. 2001c, ApJS, 134, 263 [CrossRef] [Google Scholar]
  311. Wiesenfeld, L., & Goldsmith, P. F. 2014, ApJ, 780, 183 [Google Scholar]
  312. Willacy, K., & Williams, D. A. 1993, MNRAS, 260, 635 [NASA ADS] [Google Scholar]
  313. Willacy, K., Rawlings, J. M. C., & Williams, D. A. 1994, MNRAS, 269, 921 [NASA ADS] [CrossRef] [Google Scholar]
  314. Wirström, E. S., Charnley, S. B., Cordiner, M. A., & Ceccarelli, C. 2016, ApJ, 830, 102 [CrossRef] [Google Scholar]
  315. Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  316. Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  317. Wolfram Research Inc., 2020, Mathematica, Version 12.1, Champaign, IL [Google Scholar]
  318. Woodall, J., Agündez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  319. Wu, R., Bron, E., Onaka, T., et al. 2018, A&A, 618, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  320. Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [Google Scholar]
  321. Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJ, 749, L17 [NASA ADS] [CrossRef] [Google Scholar]
  322. Zielinsky, M., Stutzki, J., & Störzer, H. 2000, A&A, 358, 723 [NASA ADS] [Google Scholar]
7

Habing (1968) provided an alternative description which is taken as the reference with scaling factor G0 in some models. To convert between the two descriptions use χD = 1.7G0.

8

If not stated otherwise we always use the perpendicular AV along the central line of sight (in contrast to e.g., the effective AV Röllig et al. 2007).

10

In a future update we will split the surface chemistry computation into dust size bins with individual values of Td and surface area for each bin.

11

The assignment of positive values to heating is arbitrary and can be reversed.

12

Numeric root search algorithms may find a root here depending on their numerical parameters even though this is mathematically not a root.

13

Some low J transitions show signs of level inversion that are unlikely to appear in practice.

14

During the first global iteration and for every first spatial step in all subsequent iterations we use a relaxed value of eps = 0.05.

15

DGESV fails primarily when the matrix ℚ is rank deficient, which is usually an effect of a stiff system of ODEs together with problematic numerics from other parts of the computation.

16

We also implemented the option to use other numerical solvers, such as the DLSODA (Hindmarsh 1983) and the DVODPK from the Netlib library dvodpk (http://www.netlib.org/ode/) and MA28, the advance solver of sparse systems of linear equations from the HSL library (http://www.hsl.rl.ac.uk/).

17

Computed by the Intel FORTRAN routine cpu_time. The time returned is summed over all active threads. The result is the sum (in units of seconds) of the current process’s user time and the user and system time of all its child processes, if any.

18

We limit ourselves to nonadaptive dampening strategies here. However, adaptive dampening strategies such as line searches and backtracking may further improve global convergence performance.

19

Notation: We denote surface reaction rate coefficients with a capital Kij and gas phase reaction rate coefficients with a lower-case kij.

All Tables

Table 1

Typical PDR model input.

In the text
Table 2

Typical PDR model output.

In the text
Table 3

Modified desorption energies used in the surface chemistry.

In the text
Table 4

Desorption energies for binding on bare grain surfaces.

In the text
Table 5

Heating and cooling processes implemented in KOSMA-→.

In the text
Table 6

Chemical species in the surface test models.

In the text
Table 7

Model parameter for models discussed in Sect. 6.

In the text
Table 8

Condensation temperatures for some relevant species in K derived from the balance between accretion and desorption.

In the text
Table 9

Comparison of model fine-structure emission (in units of K km s−1).

In the text
Table 10

Properties of a simple discrete Orion Bar ensemble.

In the text
Table 11

Observed column densities vs. model results.

In the text
Table B.1

Total PDR benchmark CPU time (in units of s) with different chemical solvers.

In the text
Table B.2

Performance comparison of various preconditioning strategies of the chemical problem

In the text
Table C.1

Total Newton-Raphson step count for different stepper choices.

In the text
Table D.1

Changes and addition to the UDfA12 reaction network

In the text
Table D.2

Updated fractionation reactions

In the text
Table D.3

Photo-desorption yields

In the text
Table D.4

Chemical desorption reactions.

In the text
Table E.1

Fine-structure and rotational line cooling data

In the text
Table F.1

Properties of a discrete example ensemble

In the text
Table G.1

WL-PDR features

In the text

All Figures

thumbnail Fig. 1

Simplified structure of a PDR. UV radiation (coming from the left) is progressively attenuated while penetrating a gas cloud. This leads to a chemical stratification and a strong temperature gradient.
In the text
thumbnail Fig. 2

General numerical scheme that needs to be solved in KOSMA-τ. Local iterations (red) are performed on every spatial grid point, spatial iterations (green) are performed over all spatial grid points in the model, and global iterations (blue) repeat the spatial loop if necessary until numerical convergence is reached.
In the text
thumbnail Fig. 3

Structure of the shell ray tracing setup. The shells are arranged equidistantly for clarity. Real spatial grids are not equidistant. The angle index t is the second index given at the blue labeled points.
In the text
thumbnail Fig. 4

Idealized form of Eq. (21) with a single root. Numerically, a temperature solution Tsol is determined by finding the root of the energy balance equation Esol = Etot(Tsol) = 0.
In the text
thumbnail Fig. 5

Idealized form of Eq. (21) with multiple roots. Filled points indicate a stable solution, open points are unstable against perturbations. The energy balance for two different positions in the PDR are schematically drawn in gray together with their equilibrium temperatures Thigh and Tlow.
In the text
thumbnail Fig. 6

Two-temperature solution for the V4 benchmark model (n = 1055 cm−3, χ = 105). Left: Tg/d (top) and ni (bottom) profile. Solid and dashed lines correspond to the low and high temperature solutions (WMM and HAM. i.e., T1 and T3 from Fig. 5). Right: Etot(T) at selected positions in the PDR (marked by black points in Fig. 6a). The top panel magnifies the range around Etot = 0.
In the text
thumbnail Fig. 7

Energy balance for the V4 model at AV = 0.75. The blue line shows the standard computation. The orange line was computed using a 20% enhanced H2 formation heating.
In the text
thumbnail Fig. 8

Influence of the surface chemistry on the thermal structure of a PDR. The model parameters are: n0 = 104 cm−3 (α = 1.5), M = 103 M⊙, χ = 106. CR1 and CR2 indicate the different models for the CR induced desorption.
In the text
thumbnail Fig. 9

Influence of the surface chemistry on the chemical structure of a PDR. The results are for the same model as in Fig. 8. Solid (black) lines in the central (right) panel correspond to CR1 while dashed (red) lines are for CR2.
In the text
thumbnail Fig. 10

Dust temperature 〈Tdust〉 profile for n0 = 104 cm−3 and different FUV fields. Sublimation temperatures of selected ice species are indicated by red marks on the right.
In the text
thumbnail Fig. 11

Left: temperature profile changes with FUV strength (n = 104 cm−3). Solid lines show pure gas-phase results, dashed lines correspond to gas+surface chemistry. Right: influence of the surface chemistry on the CO gas-phase abundance for the same parameters as in Fig. 11. CR1 and CR2 indicate the different models for the CR induced desorption.
In the text
thumbnail Fig. 12

Chemical structure changes with FUV strength (n = 104 cm−3). Solid lines indicate pure gas-phase chemistry, dashed lines correspond to the gas+surface chemistry using the CR2 model.
In the text
thumbnail Fig. 13

Percentile ice composition as function of AV for n = 104 cm−3, M = 103 M⊙ assuming CR2.
In the text
thumbnail Fig. 14

Effect of surface chemistry on the CO spectral line emission distribution (in K km s−1) as function of Jup. Each panel corresponds to a different FUV field strength χ. 12CO lines are in color, while 13CO lines are plotted in gray-level.
In the text
thumbnail Fig. 15

Ice thickness in ML for the model clump center as a function of total gas density n0 and FUV strength χ, using CR2. The numbers give the threshold for the formation of one monolayer of ice in AV.
In the text
thumbnail Fig. 16

Ice abundances relative to the H2O ice column density. The median of all models with χ = 1−103 and M = 50, 1000 M⊙ is shown.
In the text
thumbnail Fig. 17

Integrated intensity map of the HCO+ 4−3 line emission in the Orion Bar (see Goicoechea et al. 2016, for details on the mapped area). The contours corresponds to integrated intensities of 40 K km s−1 (dashed) and 60 K km s−1 (solid). The colors indicate intensities between 0 and 100 K km s−1. The two circles in the lower left show the area of clumps with M = 0.01 M⊙ (white) and M = 0.001 M⊙ (gray).
In the text
thumbnail Fig. A.1

Effect of surface chemistry on the CO spectral line emission distribution (in K km s−1) as function of Jup for n0 = 103 cm−3. Each panel corresponds to a different FUV field strength χ. 12CO lines are in color, while 13CO lines are plotted in gray-level.
In the text
thumbnail Fig. A.2

Same as Fig. A.1 for n0 = 105 cm−3.
In the text
thumbnail Fig. A.3

Same as Fig. A.1 for n0 = 106 cm−3.
In the text
thumbnail Fig. A.4

Chemical structure changes with FUV strength (n = 104 cm−3). Solid lines indicate pure gas-phase chemistry, dashed lines correspond to the gas+surface chemistry using the CR2 model.
In the text
thumbnail Fig. A.5

Chemical structure changes with FUV strength (n = 104 cm−3). Solid lines indicate pure gas-phase chemistry, dashed lines correspond to the gas+surface chemistry using the CR2 model.
In the text
thumbnail Fig. C.1

Comparison of various stepper functions. The inset shows the vertical axis in log-space. fstep,Tanh1 assumes (ω− = ω+ = 2, λ = 1), fstep,Tanh2 assumes (ω− = 20, ω+ = 5, λ = 0.5).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.