DESPOTIC – A New Software Library to Derive the Energetics and SPectra of Optically Thick Interstellar Clouds

Mark R. Krumholz

Department of Astronomy & Astrophysics, University of California, Santa Cruz, CA 95064 USA

I describe DESPOTIC, a code to Derive the Energetics and SPectra of Optically Thick Interstellar Clouds. DESPOTIC represents such clouds using a one-zone model, and can calculate line luminosities, line cooling rates, and in restricted cases line profiles using an escape probability formalism. It also includes approximate treatments of the dominant heating, cooling, and chemical processes for the cold interstellar medium, including cosmic ray and X-ray heating, grain photoelectric heating, heating of the dust by infrared and ultraviolet radiation, thermal cooling of the dust, collisional energy exchange between dust and gas, and a simple network for carbon chemistry. Based on these heating, cooling, and chemical rates, DESPOTIC can calculate clouds’ equilibrium gas and dust temperatures, equilibrium carbon chemical state, and time-dependent thermal and chemical evolution. The software is intended to allow rapid and interactive calculation of clouds’ characteristic temperatures, identification of their dominant heating and cooling mechanisms, and prediction of their observable spectra across a wide range of interstellar environments. DESPOTIC is implemented as a Python package, and is released under the GNU General Public License.

galaxies: ISM — line: profiles — methods: numerical — ISM: clouds — ISM: molecules — radiative transfer

1 Introduction

With the advent of powerful radio telescopes such as the Atacama Large Millimeter Array (ALMA), it has become possible to study the cold interstellar medium (ISM) in unprecedented detail and at greater distances than ever before. Observations from these facilities have stimulated a great deal of theoretical interest in the properties of the cold ISM, both nearby and in environments far-removed from those found near the Sun. One of the goals of these theoretical investigations has been to study how the thermodynamics of gas, and thus the nature of the star formation process within it, varies as a function of environment. A second goal has been to predict the observable emission of gas in a variety of environments.

Theoretical investigations of this sort often benefit from approximate calculations using idealized geometries that can produce relatively fast results, while also including a wide range of microphysical processes in order to determine which ones are important. However, there are few publicly-available tools capable of performing these functions for the dense, optically thick phase of the interstellar medium. Traditional photodissociation region (PDR) codes (e.g. Meijerink & Spaans, 2005; Le Petit et al., 2006), or codes that can handle a variety of ISM phases such as cloudy (Ferland et al., 1998), perform calculations in this regime, but the complexity of the problem means that these codes are too computationally-costly for either broad surveys or quick, interactive scans of parameter space. A number of authors have released codes capable of performing fast calculations of molecular emission line spectra using the large velocity gradient or various other forms of the escape probability approximation (e.g. CASSIS111 and RADEX, van der Tak et al. 2007). While these are useful tools for the analysis of observations, they are only capable of predicting line emission given fixed physical conditions, and they do not calculate many quantities of interest for theoretical modeling, such as rates of heating and cooling, thermal equilibria, or time-dependent thermal behavior.

The need for codes that are capable of performing calculations of this sort is apparent from the wide variety of applications they have found in the recent literature. For example, Goldsmith (2001) and Lesaffre et al. (2005) investigate the temperature structure within protostellar cores. Krumholz & Thompson (2007) use an escape probability model to study the relationship between star formation rates and emission is a variety of molecular lines. Krumholz, Leroy & McKee (2011) use thermal equilibrium models of the ISM to explore the relationship between star formation and the chemical state of the gas. Narayanan et al. (2011, 2012), Shetty et al. (2011a, b), and Feldmann, Gnedin & Kravtsov (2012a, b) all investigate the conversion between observed CO luminosity and molecular mass using simulations of galaxies, coupled to post-processing to predict the observable line emission. Narayanan & Davé (2012a, b) perform calculations of interstellar medium temperatures as a way of estimating the Jeans mass in molecular clouds, and its possible implications for changes in the stellar IMF over cosmological times. Papadopoulos (2010) and Meijerink et al. (2011) consider star formation in extreme environments with X-ray and cosmic-ray fluxes far higher than are found in the Solar neighborhood, and in the process rely on calculations of the thermal behavior of gas under these conditions. Similarly, Muñoz & Furlanetto (2013) study the ISM in high-redshift galaxies where the metallicity is much lower and the cosmic microwave background is much hotter than in the present-day universe. With a few exceptions, all of these authors developed their own custom codes to model the thermodynamics and line emission of the cold ISM. However, this effort is largely duplicative, since these calculations all involve the same related set of problems. Moreover, the results of the calculations can be difficult to compare due to the differing assumptions and approximations made by the various authors in their modeling, not all of which are well-documented in the literature.

In order to support theoretical investigations facing problems of this sort, reduce duplication of effort, and encourage calculations with documented, open-source tools to allow easy comparison between authors, I have developed a software library to Derive the Energetics and SPectra of Optically Thick Interstellar Clouds (DESPOTIC). DESPOTIC uses an escape probability formalism to calculate line emission, and couples this to a calculation of either equilibrium or time-dependent gas and dust temperatures, including the dominant processes in a wide variety of environments: cosmic-ray and X-ray ionization heating, photoelectric heating, grain-gas energy exchange, and radiative heating and cooling of dust grains. The software is implemented as Python package, enabling easy, interactive calculation, and also easy integration with other software. It also provides an automated interface with the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005). DESPOTIC is publicly available from a dedicated web page, from google code, and from the Python Package Index, and is released under the GNU General Public License. It comes with extensive documentation, including a User’s Guide with a full listing of all routines and options.

In the remainder of this paper, I describe the model system that DESPOTIC uses and the equations it can solve (§ 2) and the numerical methods by which it solves those equations (§ 3). I then provide some example applications (§ 4), provide some warnings about the limitations of the code (§ 5), and summarize (§ 6).

2 Model System and Equations Solved

2.1 Physical Model

The basic physical system treated by DESPOTIC is a uniform spherical cloud (though other simple geometries are provided as options, as described below). Such a cloud is characterized by several physical and chemical properties, which are taken to be uniform unless stated otherwise. The physical properties are a volume density of hydrogen nuclei and a mean column density of hydrogen nuclei , the gas temperature , the dust temperature , the non-thermal velocity dispersion , and (optionally) a bulk radial velocity gradient . Note that DESPOTIC defines the column density as an average over the cloud, i.e. it is the total number of hydrogen atoms in the cloud divided by the cloud’s cross-sectional area.

The dust within a cloud is characterized by six quantities. Three of these describe the dust cross-section per H nucleus to thermal radiation at temperature T = 10 K, , to radiation in the range of eV that dominates photoelectron production, , and averaged over the diffuse interstellar radiation field (ISRF) . The fourth quantity is the total dust abundance normalized to the Milky Way value, . The remaining two quantities are the dust spectral index for thermal radiation, and the gas-grain collisional coupling coefficient . I define all of these terms in detail below.

DESPOTIC parameterizes the radiation field (including cosmic rays) around the cloud by the following quantities: gives the primary ionization rate per H nucleus due to hard x-ray photons and cosmic rays, describes the energy density, normalized to the Solar neighborhood value, of the non-thermal interstellar radiation field produced primarily by stars, gives the infrared radiation field seen by the dust, and is the cosmic microwave background temperature.

Finally, the chemical composition of the cloud is described by the abundances of bulk constituents and trace emitting species. The abundances of the bulk constituents in the DESPOTIC model are given by , , , , , and , which describe atomic hydrogen, para-H, ortho-H, helium, free electrons, and free protons, respectively.222Although DESPOTIC includes free protons and electrons, it is only intended for use in regions where the gas is predominantly neutral, i.e.  , and similarly for . It does not include many heating and cooling processes that are important in highly ionized regions. The abundances of emitting species (e.g. CO, HCN, HO, etc.) are characterized in the same way, with representing the abundance of the th emitting species.

Given the bulk composition, one can also compute a number of additional quantities, of which we will make use below. Three of these are the mean mass per H nucleus , the mean mass per free particle , and the isothermal sound speed , given by


where and are measured in units of the hydrogen mass . Note this this expression neglects the mass of electrons, and assumes that emitting species contribute negligibly to the mass. Two additional quantities are the gas specific heat at constant volume and at constant pressure , which for convenience we express per H nucleus rather than per unit mass or per unit volume. Thus and have units of energy over temperature, and can be converted to the usual values per unit mass simply by multiplying by a . Calculation of the specific heats requires some care when the chemical composition includes molecular hydrogen. I discuss this topic in detail, and derive DESPOTIC’s expressions for and , in Appendix A.

The final quantity one can compute is the clumping factor for the cloud, which represents an enhancement in the rates of all collisional processes due to non-uniformity of the gas. The quantity is the volume-averaged density over the cloud, but in a non-uniform cloud the density at any position may be higher or lower than this. Since the rate of collisions per unit volume at a given position varies as , the rate of collisions per H atom in a non-uniform cloud exceeds that in a uniform cloud by a factor


where the angle brackets indicate an average over the cloud volume; thus is simply the factor by which the mass-weighted mean density exceeds the volume-weighed mean density. For a supersonically turbulent medium, this factor is approximately (Ostriker, Stone & Gammie 2001; Padoan & Nordlund 2002; also see Lemaster & Stone 2008, Federrath, Klessen & Schmidt (2008), and Price, Federrath & Brunt 2011)


2.2 Heating and Cooling Processes

The gas heating and cooling processes included in DESPOTIC are ionization heating, heating by the grain photoelectric effect, gravitational compression heating, line cooling, and either heating or cooling by collisional energy exchange between dust and gas. The grain heating and cooling processes included in DESPOTIC are cooling by thermal radiation, heating by the interstellar radiation field, heating by an infrared radiation field, heating by the cosmic microwave background radiation, heating by absorption of line radiation, and collisional coupling to the gas.

Given this list of processes, the time rate of change of the gas energy per H nucleus as


where , , and are the rates of ionization, photoelectric, and gravitational heating per H nucleus, is the rate of line cooling per H nucleus, and is the rate of dust-gas energy exchange per H nucleus. I give explicit formulae for all these terms in Appendix B. The corresponding time rate of change of the temperature is


where is the gas specific heat per H nucleus at constant volume and is the specific heat per H nucleus at constant pressure (which are calculated in Appendix A). The parentheses indicate that one can use either or in the above equation, depending on whether one wishes to consider gas cooling isochorically or isobarically.

Similarly, for the dust grains the total rate of change of specific energy per H nucleus is


where , , , and are the rates of heating due to the interstellar radiation field, line radiation, the cosmic microwave background, and infrared radiation, and is the rate of dust cooling by thermal radiation. As with the gas heating and cooling processes, I give explicit formulae for all these terms in Appendix B. In principle one could consider time-dependent temperature evolution of the dust as well as of the gas, but since the specific heat of the dust is far less than that of the gas, and is a complex function of the properties of the grains, DESPOTIC does not treat this case. Instead, it assumes that the grain population is always in thermal equilibrium.

2.3 Chemical Processes

In addition to thermal processes, DESPOTIC can also calculate chemical processes that cause the abundances of various species to change with time. DESPOTIC allows users to define arbitrary chemical networks by specifying a set of species and a set of chemical reaction rate equations of the form


where is the vector of fractional abundances for the various species in the network, and the reaction rates on the right-hand side can be a function of these abundances, of the overall volume density, column density, gas temperature, ionization rate, radiation field, or any of the other quantities that DESPOTIC uses to describe a cloud. Once specified, the equations can be integrated over a specified time or until the chemical state reaches equilibrium. The repository version of DESPOTIC implements the reduced carbon-oxygen chemistry network of Nelson & Langer (1999), which models the processes leading to the transition from C- to CO-dominated composition in molecular clouds.

2.4 Line Shapes

Figure 1: Diagram of the geometry used by DESPOTIC when calculating line shapes. The circle shows the cloud, with radius , and the dashed line is the observer’s line of sight through it.

DESPOTIC’s final major capability is calculating the profiles of spectral lines. In general this is not a useful calculation in a one-zone escape probability model; since the level populations in such a model are assumed to be uniform, the result is necessarily rather uninteresting, and is simply given by the usual solution to the transfer equation for media with emission and absorption coefficients that are independent of position. However, one can relax the assumption of uniform level populations by making another one: that the species is in LTE, and that the temperature is a known function of position.333In principle one in fact needs to know only the excitation temperature for the two levels that produce the line, together with the number density the atoms or molecules that are in the lower state . However, in practice it is unlikely that one will simultaneously know and in any situation other than when the levels are in LTE, and thus I limit the discussion to this case. If one does in fact know and , it is trivial to perform a calculation for that case simply by setting the level populations to their LTE values at , and adjusting the overall density of the species so that has the desired value. Solving for the shapes of lines in this limit allows the code to compute pCygni and inverse pCygni profiles, among other applications. This computation is performed for a spherical cloud following DESPOTIC’s general model, and consider a line of sight passing through it at an offset distance from the cloud center (see Figure 1). Details of how this calculation is performed are given in Appendix C.

3 Code Architecture and Algorithms

In this section I describe the architecture of the DESPOTIC code and the algorithms it uses to solve the equations introduced in the previous section.

3.1 Overall Architecture

DESPOTIC is a library intended not only to be used for stand-alone calculations, but also to allow easy extensibility, easy integration with other codes, and to allow users to conduct interactive, exploratory calculations. To this end, DESPOTIC is implemented as a Python package, which allows a very high level of abstraction such that many useful computations can be performed with no more than a single line of code on the part of the user. To achieve high performance, DESPOTIC makes extensive use of the ability of the numPy and sciPy libraries to interface with the fast, optimized numerical libraries LAPACK444 (Anderson et al., 1999), MINPACK555 (Moré, Garbow & Hillstrom, 1980), and ODEPACK666 (Hindmarsh, 1983). It is hard to provide a quantitative estimate of code execution times for DESPOTIC routines, since as I discuss below the most computationally-intensive ones require iterative methods, and the time required for such a solution is a strong function of the quality of the starting guess. Nonetheless, I give a general idea of code execution times, as tested on a single processor of a modern workstation, for some example applications in § 4. Individual instances of DESPOTIC classes use internal private storage, and thus are thread-safe should a user desire to use threading to accelerate the calculation of large grids of models via the standard Python threading interface. Threading of internal DESPOTIC calculations for single clouds will be added in a future release.

3.2 Capabilities and Algorithms

3.2.1 Level Populations and Line Luminosities

The most basic capability of DESPOTIC is to compute level populations and line luminosities for an emitting species embedded in a cloud of specified physical properties (, , , abundances, etc.). The emitted intensity for any line is given by equation (65), and the numerical algorithm for calculating level populations and line luminosities is given in Appendix B.3. The computation can be performed either assuming the cloud is optically thin, or using the escape probability approximation for an optically thick cloud. Note that this is the same computation performed by codes like RADEX (van der Tak et al., 2007) and lineLum (Krumholz & Thompson, 2007), and the latter is the direct ancestor of the corresponding portion of DESPOTIC. Appendix D provides a direct comparison between DESPOTIC and RADEX.

3.2.2 Cooling Rates, Thermal Equilibria, and Time-Dependent Temperature Evolution

In addition to computing line luminosities and level populations, DESPOTIC can also compute the heating and cooling rates of gas and dust. It does so by evaluating all the terms in equations (6) and (8); since one of these terms is , this procedure entails solving for the level populations and escape probabilities.

DESPOTIC can also solve for equilibrium dust and gas temperatures. DESPOTIC obtains these values by setting in equation (6) and in equation (8). The user can also add arbitrary additional heating and cooling terms to either equation, to represent processes not modeled by DESPOTIC (e.g. endothermic or exothermic chemical reactions). At the discretion of the user, DESPOTIC can fix either or and solve for the other, or it can solve for both simultaneously. If either or is fixed, DESPOTIC solves the equations using the secant method. If neither is fixed, it solves for and simultaneously using the MINPACK routine hybrd1, which implements the Powell hybrid method.

Finally, DESPOTIC can compute the time-dependent thermal evolution of a cloud. Starting from an initial gas and dust temperature, DESPOTIC can integrate equation (7) for the gas temperature evolution. At the user’s discretion, the calculation can be done either isochorically or isobarically. When evaluating the heating and cooling terms that appear on the right-hand side of equation (7), DESPOTIC assumes that both the level populations and the dust temperature reaches equilibrium instantaneously; the former are computed via the procedure described in Appendix B.3, and the latter by the solution to equation (8) with . DESPOTIC also calculates the temperature-dependent specific heat or on the right-hand side using equation (19). It integrates equation (7) using the ODEPACK routine lsoda, which automatically evaluates the stiffness of the system, and solves using a predictor-corrector method for non-stiff problems and backward differentiation formula methods for stiff problems.

3.2.3 Chemical Evolution and Chemical Equilibria

DESPOTIC’s implementation of chemistry has two parts. First, DESPOTIC provides a series of routines that can integrate the chemical evolution equations (9), either for a specified time interval or until the rates of change of all abundances are zero to within some specified tolerance. Second, DESPOTIC provides a generic interface that can be used to implement arbitrary chemical networks. Once implemented, one can use the chemical evolution routines to integrate that network in an automated fashion. One basic network, that of Nelson & Langer (1999), is included in the code repository.

3.2.4 Line Profiles

DESPOTIC’s final major capability is calculating line profiles for species in LTE. When performing this calculation, it accepts user-specified profiles for the number density of the emitting species, the bulk velocity, the non-thermal velocity dispersion, and the temperature as a function of radius. From these inputs, plus the identity of the line whose profile is to be computed, it calculates all the dimensionless quantities given in equations (78) – (80), and then numerically integrates equation (74) at a range of user-specified frequencies or velocities. The integration is performed via a call to the ODEPACK routine lsoda. DESPOTIC then returns the CMB-subtracted intensity and brightness temperature as a function of frequency / velocity.

3.3 Atomic and Molecular Data

DESPOTIC obtains the chemical data required for its computations (e.g., Einstein coefficients, reaction rate coefficients) from the Leiden Atomic and Molecular Database (LAMBDA; Schöier et al. 2005). Access to the database is automated: DESPOTIC automatically fetches whatever data files are needed without explicit user intervention. DESPOTIC makes three approximations in situations where data from LAMDA is not available. First, for some species, LAMDA provides estimates only of collision rate coefficients for H, not for oH and pH separately, or it provides only oH or pH. In such cases, DESPOTIC assumes that the oH and pH collision rate coefficients are equal, and, if only generic H rates are given, it sets both of them equal to those.

Second, for some species collision rate coefficients for H are available, but collision rate coefficients for He are not. In this case DESPOTIC assumes that He collision rate coefficients are related to those for H by (Schöier et al., 2005)


where is the reduced mass of the species with H, and similarly for .

Third, by default DESPOTIC will not extrapolate collision rates outside the range of temperatures provided in the LAMDA tables. However, the user can override this default behavior, in which case DESPOTIC will extrapolate by assuming that the downward collision rate coefficient varies as a powerlaw in the gas kinetic temperature. For linear molecules, a more accurate extrapolation motivated by a quantum mechanical treatment of the collision is possible (see the Section 6 of Schöier et al. 2005), but no such treatment is available for non-linear molecules.

4 Sample Applications

In this section I provide some sample applications to demonstrate DESPOTIC’s capabilities. Each of these applications operates on one or more example clouds, whose properties are specified in Table 1. The values given in this Table are intended to be examples only, but input files corresponding to each of them are included with the DESPOTIC library to provide example templates that users can modify to set up their own clouds. The code to perform each of the example calculations listed below is also included with the DESPOTIC download.

Cloud Name MilkyWayGMC ULIRG ProtostellarCore PostShockSlab
Physical Properties
[km s] 2.0 80.0 0.1 0.5
[K] 8 45 8 250
[K] 8 60 8 8
0.0 0.0 0.0 0.0
0.1 0.1 0.1 0.1
0.4 0.4 0.4 0.4
0.1 0.1 0.1 0.1
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
Dust Properties
[erg cm K]
[cm H]
[cm H]
[cm H]
1.0 1.0 1.0 1.0
2.0 2.0 2.0 2.0
Radiation Field Properties
[K] 2.73 2.73 2.73 2.73
[K] 0.0 60.0 8.0 8.0
[s H]
1.0 1.0 1.0
Emitting Species Abundances
CO - -
C - -
O - -
CS - -
HCO - -
pNH - -
oNH - -
pHCO - -
oHCO - -
pHO - -
oHO - -

The table gives initial properties for the example cloud models used in § 4. For applications where and are fixed, the values given in the table are the values used; for applications where and are to be calculated, they are used as initial guesses. For the ProtostellarCore model, the density is given as a range because a range of models are run. The abundances in this model have been chosen to roughly match those recommended in Goldsmith (2001). For the PostShockSlab model, emitting species marked with asterisks indicate species from which line emission is computed, but that are ignored for the purposes of calculating the thermal evolution. The molecular data from LAMDA used in evaluating these models are taken from the following sources: CO, CO, and CO: Yang et al. (2010); C: Schroder et al. (1991) and Staemmler & Flower (1991); O: Jaquet et al. (1992); CS: Turner et al. (1992); HCO: Flower (1999); NH: Danby et al. (1988); HCO: Green (1991); HO: Daniel, Dubernet & Grosjean (2011).

Table 1: Sample clouds

4.1 CO Spectral Line Energy Distributions

Figure 2: Spectral line energy distribution for the first 8 rotational transitions of CO and CO, computed for the models MilkyWayGMC and ULIRG described in Table 1. The plot shows the velocity-integrated brightness temperature in each line normalized by cm. The contribution of the CMB has been subtracted off.

As a first example of DESPOTIC’s capabilities, Figure 2 shows a calculation of CO and CO spectral line energy distributions (SLEDs) for the MilkyWayGMC and ULIRG clouds described in Table 1. For this computation, the gas temperature is left fixed to the input value, and the level populations are computing using the escape probability formalism. As expected, all lines of the ULIRG are much brighter due to its higher gas kinetic temperature and velocity dispersion – to first order, the velocity-integrated brightness temperature of an optically thick line is simply the product of those two. In addition, the falloff in luminosity with is much slower for the ULIRG than for the Milky Way cloud. This is as a result of the much higher density and temperature of the ULIRG. The former allows its higher levels to be close to thermally populated, and the latter causes their thermal populations to be large. We also see that the CO(1-0) to CO(1-0) ratio is larger for the ULIRG than for the Milky Way model, reflecting the higher optical depth of the ULIRG. At higher , where the optical depth drops, the line ratios of the two isotopomers vary less between the two models.

Note that this computation for CO(1-0) is equivalent to calculating the CO “X-factor" that relates CO intensities to cloud masses and column densities. The values calculated by DESPOTIC are cm for MilkyWayGMC, and cm. This is in line with other theoretical and observational estimates for normal galaxies and ULIRGs, respectively (Bolatto, Wolfire & Leroy, 2013). One should be wary of reading too much into this result, since neither the chemical and thermal states of the clouds have been specified by hand. This computation should be done by combining a three-dimensional simulation with chemical post-processing to determine the chemical state of the clouds self-consistently, and then using DESPOTIC or a similar package to calculate the resulting gas temperature and line radiation (e.g. Narayanan et al., 2011, 2012; Shetty et al., 2011a, b; Feldmann, Gnedin & Kravtsov, 2012a, b).

4.2 Temperatures of Protostellar Cores

Figure 3: Equilibrium gas and dust temperatures versus density for the ProtostellarCore model described in § 4.2.
Figure 4: Values of the various heating and cooling terms for dust and gas in the ProtostellarCore models, calculated at the equilibrium temperatures shown in Figure 3. The panels show gas heating terms, gas cooling terms, dust heating terms, and dust cooling terms, as indicated in the legends. The terms shown are gas ionization heating, , dust IR heating, , dust ISRF heating , dust line heating , gas line cooling, , dust thermal cooling, , and dust-gas energy exchange, . The calculations also include gravitational and photoelectric heating, but these terms are below the plotted range.
Figure 5: Contributions to the overall line cooling rate for individual atomic molecular species in the ProtostellarCore model. All line cooling rates are computed at the equilibrium temperatures shown in Figure 3. Note that these rates are computed for constant abundances, and thus do not properly account for depletion at high densities. They are therefore likely to be overestimates at the high-density end, as discussed in § 4.2.

As a second example application, I use DESPOTIC to calculate the equilibrium gas and dust temperatures in protostellar cores as a function of density, using the algorithms outlines in § 3.2.2. In this calculation I include a large number of cooling species (see Table 1) in order to assess their density- and temperature-dependent contribution to cores’ thermal balance. For this calculation I use the ProtostellarCore model in Table 1. I then compute a grid of models with densities in the range cm in steps of 0.2 dex. For each model, I compute the equilibrium gas and dust temperatures, and, once the equilibrium has been calculated, I record the values of all the heating and cooling terms.

Figure 3 shows the equilibrium temperatures as a function of density, Figure 4 shows the contributions of the various heating and cooling processes, and Figure 5 further subdivides the line cooling into the contributions made by individual species. The plots illustrate a number of phenomena. First, the gas temperature is relatively high at low densities, and drops as the density increases. At densities below cm this drop is driven by increasingly effective line cooling. Between and cm, dust-gas collisions become competitive with line cooling, and lock the dust and gas temperature together, such that dust-gas energy exchange becomes dominant in setting the temperature. The dust in turn is always locked close to the infrared radiation field temperature, because the IR heating rate and thermal cooling rate both exceed all other sources and sinks of energy for the dust by orders of magnitude.

In terms of molecular line cooling, at low densities the dominant coolants are CO, CO, and C. As the density rises and the dust temperature drops, these become less important because dust coupling lowers the gas temperature. This makes it more difficult to excite the higher lines that have lower optical depths. At the same time, other species make an increasing contribution to the cooling as the density approaches their critical densities and begins to provide efficient collisional excitation. However, this example also illustrates one of DESPOTIC’s limitations. These calculations assume density- and temperature-independent abundances, and do not properly model the effects of freeze-out onto grain surfaces. Over the density range I have explored freeze-out is probably significant only for CS, since sulfur-bearing molecules begin to freeze out at densities above cm, but carbon- and nitrogen-bearing ones do not experience significant freeze-out until the density rises above cm and cm, respectively (e.g. Bergin & Langer, 1997). Once could include freeze-out effects by defining an appropriate chemistry network, but the simple Nelson & Langer (1999) network that DESPOTIC currently implements does not model these effects.

This is the most computationally-intensive of the example applications provided, due to the high optical depth and the large number of molecular coolants included. The majority of the computational effort involves iterating to obtain the level populations at high optical depth. Evaluating the entire grid of 21 models requires a bit under 5 minutes. However, since only a few chemical species are actually important to the thermal balance, one could obtain the results far more quickly simply by ignoring the large number of energetically-unimportant species when calculating the temperature, and only calculating their line luminosities once the temperatures have converged. DESPOTIC includes a capability to mark certain species as energetically-unimportant, allowing them to be treated in precisely this manner, and I demonstrate this capability in the next example.

4.3 Time-Dependent Cooling of Post-Shock Gas

A third example, which makes use of DESPOTIC’s ability to calculate time-dependent temperature evolution (§ 3.2.2), is to calculate the cooling of out-of-equilibrium gas. I consider a slab of gas whose properties are given by the PostShockSlab model in Table 1. At time , the gas has just been shock-heated to an out-of-equilibrium temperature of K, and I calculate the time evolution of its temperature and line emission thereafter, assuming that the gas is isobaric and using a slab geometry to compute escape probabilities. In calculating the thermal evolution I include only the energetically-dominant coolants CO, CO, and O, but I also periodically compute the line emission of a large number of other species as well (see Table 1 for the full list). By making this assumption, the total computer time required to evolve the model 40 kyr, including periodic calculation of emission from many lines, is minutes.

Figure 6: Gas temperature, dust temperature, and gas density versus time for isobaric cooling of the PostShockSlab model, as described in § 4.3.
Figure 7: Rates of cooling provided by CO lines, CO lines, O lines, and dust versus time, for the PostShockSlab model shown in Figure 6. The gray thick line shows the sum of all coolants, including all a number of lines that are not shown because they all below the range of cooling rates plotted.
Figure 8: CO line spectral line energy distribution for the PostShockSlab model shown in Figure 6. Each of the lines shows the relative contributions of the indicated rotational transitions of CO to the total cooling rate at the indicated times of 0 kyr, 20 kyr, and 40 kyr. Contributions are normalized so that the sum over all transitions is unity.

Figure 6 shows the gas temperature, dust temperature, and gas density versus time as computed by DESPOTIC for this initial condition. Figure 7 shows the contributions of various species to the cooling. As the plot shows, cooling is dominated by CO lines, with minor contributions from CO, O, and dust, and negligible contributions from all other sources. In Figure 8 I further examine the cooling, by showing how the CO spectral line energy distribution changes with time. As the plot shows, the SLED initially peaks near , and moves to a cooler SLED at time passes. At the final time shown, is the dominant coolant. Note that this differs from the result shown in Figure 2 for a typical GMC because the post-shock slab we are considering has a significantly lower velocity dispersion and a significantly higher density. Both of these favor cooling through higher lines, the former because it increases the optical depth for low lines, and the latter because it helps to thermalize higher states.

4.4 Carbon-Oxygen Chemistry

Figure 9: Abundances of various carbon- and oxygen-bearing species as a function of gas column column density, computed using DESPOTIC’s chemistry module following the procedure described in Section 4.4.

The next sample application demonstrates DESPOTIC’s chemistry capability. For this test, I use the physical parameters for the MilkyWayGMC model listed in Table 1, except that I reduce the cosmic ray ionization rate to s and raise the gas temperature to 10 K. I then consider a range of column densities cm, in steps of dex, and use the Nelson & Langer (1999) chemical network implemented in DESPOTIC to calculate the equilibrium chemical state of the cloud. The column density range is chosen to model the transition from a composition dominated by H, C, and O (the so-called “dark gas" – Wolfire, Hollenbach & McKee (2010)) to one dominated by H and CO, as found in the interiors of molecular clouds. The total execution time of the calculation, which involves finding the equilibrium chemical state 201 times, was roughly a minute.

Figure 9 shows the result. As shown in the figure, at low column density the chemical state is such that the carbon is mostly C and the oxygen is mostly O. As the column density increases, the carbon shifts into C and CO as the dominant states, while the oxygen also shifts into CO and OH. The chemical transition is driven by the decreasing rate of photodissociation as the column density increases, which, following Nelson & Langer’s prescription, is handled using the tabulated shielding functions of van Dishoeck & Black (1988). This test demonstrates DESPOTIC’s ability to perform limited astrochemistry calculations.

4.5 Inverse P Cygni Profiles

As a final application, I use DESPOTIC to calculate line profiles in a collapsing protostellar core. For this example, I consider a core with a radius of pc with a velocity profile km s. The temperature profile is K, so that the temperature reaches a peak of 20 K at the center, dropping close to 8 K at large radii. The core also has a position-independent non-thermal velocity dispersion of km s. I use a uniform density cm.

For this core I use DESPOTIC to compute the profiles of the HCN(1-0) and NH lines. This combination of lines is often used to measure infall motions (e.g., Sohn et al., 2007), as their overall spatial distributions in a protostellar core are thought to be quite similar on chemical grounds, but the HCN(1-0) tends to be marginally optically thick and develop inverse P Cygni profiles, while the NH tends to be optically thin and show symmetric profiles. I adopt abundances and based on the models of Lee, Bergin & Evans (2004).

Figure 10: Brightness temperature versus velocity relative to line center for the lines HCN() and NH() produced by a collapsing protostellar core. The contribution of the CMB has been subtracted off. Details of the core parameters are given in § 4.5.

Figure 10 show the line profiles computed by DESPOTIC. As expected, the marginally optically thick HCN line produces a double-peaked asymmetric inverse P Cygni profile, indicative of infall. The NH line is optically thin and produces a symmetric profile of lower total intensity. The total time required to perform the computation is s.

5 Limitations and Caveats

While DESPOTIC provides reasonable estimates of the thermal behavior and spectra of interstellar clouds over a wide range of environments, it also has significant limitations, which I discuss here as a warning to potential users. The major limitations of the code are:

  • DESPOTIC’s treatment of dust temperatures is very crude in the regime of clouds that are optically thick to their own cooling radiation. In such clouds the dust temperature will be determined largely by the value of that the user selects. If a user requires accurate dust temperatures in such clouds, he or she is advised to use a code like dusty (Ivezic & Elitzur, 1997) to calculate the dust temperature and radiation field within the cloud, then use this to set for the purposes of a DESPOTIC calculation.

  • DESPOTIC neglects the contribution of the dust radiation field to the photon occupation number when calculating level populations, on the grounds that, because dust optical depths are small at low frequencies, such fields are often highly sub-thermal at the low frequencies where most important molecular lines lie. However, in some circumstances, e.g. protostellar disks (Krumholz, Klein & McKee, 2007), the column density is so high that dust optical depths can exceed unity even at frequencies as low as GHz. In such environments excitation and de-excitation of molecules by interaction with the infrared field is non-negligible, and DESPOTIC will not give accurate results.

  • DESPOTIC uses a one-zone model, and this is not capable of capturing effects that depend on radiative transfer. In particular, DESPOTIC cannot handle maser emission, and it cannot handle effects on the line shape that arise from spatially-variable departures from LTE.

  • The repository version of DESPOTIC includes only a single, very simple chemical network. This can be used to make reasonable predictions for carbon and oxygen chemistry in H-dominated environments, but not for other species or in other environments. It is up to the user to either input chemical abundances directly, or to implement chemical networks appropriate for the environment he or she wishes to simulate. The results DESPOTIC produces will only be as good as those abundances or networks. More subtly, DESPOTIC does include the effects of selective chemical destruction of excited states on line emission, and it does not include any heating or cooling of the gas or dust as a result of chemical reactions, such as heating of dust grains by exothermic formation of H on grain surfaces (e.g. Lesaffre et al., 2005). DESPOTIC provides a mechanism to include chemical heating and cooling, since the user can specify arbitrary additional heating and cooling terms, but it is up to the user to determine whether there are any energetically-important chemical reactions for the problem under consideration, and, if so, to implement the necessary code.

6 Summary

I introduce DESPOTIC, a Python-based, open-source software library for calculating spectra, heating and cooling rates, and time-dependent and time-independent thermal properties of optically thick interstellar clouds. DESPOTIC includes all the dominant heating and cooling processes for both gas and dust over a wide range of interstellar environments, and can be used to conduct both fast sweeps of parameter space and interactive explorations within an interactive Python environment. It is intended to allow theoretical investigators to obtain approximate values of parameters such as cloud temperatures, major heating and cooling processes, and observable line emission, without the difficulty and time investment of developing their own statistical and thermal equilibrium codes, and with significantly less investment of CPU and human time than would be required to approach such problems using a detailed PDR code. DESPOTIC is under continued development, and additional features capabilities will be released to the community as they are implemented.


I thank the creators and maintainers of the Leiden Atomic and Molecular Database, F. Schöier, F. van der Tak, E. van Dishoeck, and J. Black, for providing that valuable resource. I thank B. Draine for helpful suggestions regarding modeling of dust, and F. van der Tak for helpful suggestions on the manuscript and advice on RADEX. I acknowledge support from the Alfred P. Sloan Foundation, the NSF through CAREER grant AST-0955300, and NASA through Astrophysics Theory and Fundamental Physics Grant NNX09AK31G.


  • Anderson et al. (1999) Anderson E. et al., 1999, LAPACK Users’ Guide, 3rd edn. Society for Industrial and Applied Mathematics, Philadelphia, PA
  • Bakes & Tielens (1994) Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822
  • Bergin & Langer (1997) Bergin E. A., Langer W. D., 1997, ApJ, 486, 316
  • Black & Bodenheimer (1975) Black D. C., Bodenheimer P., 1975, ApJ, 199, 619
  • Bolatto, Wolfire & Leroy (2013) Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, in press, arXiv:1301.3498
  • Boley et al. (2007) Boley A. C., Hartquist T. W., Durisen R. H., Michael S., 2007, ApJ, 656, L89
  • Chakrabarti & McKee (2005) Chakrabarti S., McKee C. F., 2005, ApJ, 631, 792
  • Dalgarno & McCray (1972) Dalgarno A., McCray R. A., 1972, ARA&A, 10, 375
  • Dalgarno, Yan & Liu (1999) Dalgarno A., Yan M., Liu W., 1999, ApJS, 125, 237
  • Danby et al. (1988) Danby G., Flower D. R., Valiron P., Schilke P., Walmsley C. M., 1988, MNRAS, 235, 229
  • Daniel, Dubernet & Grosjean (2011) Daniel F., Dubernet M.-L., Grosjean A., 2011, A&A, 536, A76
  • de Jong (1977) de Jong T., 1977, A&A, 55, 137
  • de Jong, Boland & Dalgarno (1980) de Jong T., Boland W., Dalgarno A., 1980, A&A, 91, 68
  • de Jong, Dalgarno & Chu (1975) de Jong T., Dalgarno A., Chu S.-I., 1975, ApJ, 199, 69
  • Dislaire et al. (2012) Dislaire V., Hily-Blant P., Faure A., Maret S., Bacmann A., Pineau Des Forêts G., 2012, A&A, 537, A20
  • Draine (2003) Draine B. T., 2003, ApJ, 598, 1017
  • Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton University Press: Princeton, NJ
  • Federrath, Klessen & Schmidt (2008) Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688, L79
  • Feldmann, Gnedin & Kravtsov (2012a) Feldmann R., Gnedin N. Y., Kravtsov A. V., 2012a, ApJ, 747, 124
  • Feldmann, Gnedin & Kravtsov (2012b) Feldmann R., Gnedin N. Y., Kravtsov A. V., 2012b, ApJ, 758, 127
  • Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, Proc. Astron. Soc. Pac., 110, 761
  • Flower (1999) Flower D. R., 1999, MNRAS, 305, 651
  • Glassgold, Galli & Padovani (2012) Glassgold A. E., Galli D., Padovani M., 2012, ApJ, 756, 157
  • Glassgold & Langer (1973) Glassgold A. E., Langer W. D., 1973, ApJ, 186, 859
  • Goldsmith (2001) Goldsmith P. F., 2001, ApJ, 557, 736
  • Green (1991) Green S., 1991, ApJS, 76, 979
  • Hindmarsh (1983) Hindmarsh A. C., 1983, in IMACS Transactions on Scientific Computation, 10th IMACS world Congress on Systems Simulation and Scientific Computation, Vol. 1, Scientific Computing, Applications of Mathematics and Computing to the Physical Sciences, Stepleman, Carver, Peskin, Ames, Vichnevetsky, eds., North Holland, Amsterdam, pp. 55–64
  • Ivezic & Elitzur (1997) Ivezic Z., Elitzur M., 1997, MNRAS, 287, 799
  • Jaquet et al. (1992) Jaquet R., Staemmler V., Smith M. D., Flower D. R., 1992, Journal of Physics B Atomic Molecular Physics, 25, 285
  • Krumholz, Klein & McKee (2007) Krumholz M. R., Klein R. I., McKee C. F., 2007, ApJ, 665, 478
  • Krumholz, Leroy & McKee (2011) Krumholz M. R., Leroy A. K., McKee C. F., 2011, ApJ, 731, 25
  • Krumholz & Thompson (2007) Krumholz M. R., Thompson T. A., 2007, ApJ, 669, 289
  • Le Petit et al. (2006) Le Petit F., Nehmé C., Le Bourlot J., Roueff E., 2006, ApJS, 164, 506
  • Lee, Bergin & Evans (2004) Lee J.-E., Bergin E. A., Evans, II N. J., 2004, ApJ, 617, 360
  • Lemaster & Stone (2008) Lemaster M. N., Stone J. M., 2008, ApJ, 682, L97
  • Lesaffre et al. (2005) Lesaffre P., Belloche A., Chièze J., André P., 2005, A&A, 443, 961
  • Masunaga, Miyama & Inutsuka (1998) Masunaga H., Miyama S. M., Inutsuka S., 1998, ApJ, 495, 346
  • Meijerink & Spaans (2005) Meijerink R., Spaans M., 2005, A&A, 436, 397
  • Meijerink et al. (2011) Meijerink R., Spaans M., Loenen A. F., van der Werf P. P., 2011, A&A, 525, A119+
  • Moré, Garbow & Hillstrom (1980) Moré J. J., Garbow B. S., Hillstrom K. E., 1980, User Guide for MINPACK-1. Tech. Rep. ANL-80-74, Argonne National Laboratory
  • Muñoz & Furlanetto (2013) Muñoz J. A., Furlanetto S. R., 2013, MNRAS, submitted, arXiv:1301.0619
  • Narayanan & Davé (2012a) Narayanan D., Davé R., 2012a, MNRAS, 423, 3601
  • Narayanan & Davé (2012b) Narayanan D., Davé R., 2012b, MNRAS, submitted, arXiv:1210.6037
  • Narayanan et al. (2011) Narayanan D., Krumholz M., Ostriker E. C., Hernquist L., 2011, MNRAS, 418, 664
  • Narayanan et al. (2012) Narayanan D., Krumholz M. R., Ostriker E. C., Hernquist L., 2012, MNRAS, 421, 3127
  • Nelson & Langer (1999) Nelson R. P., Langer W. D., 1999, ApJ, 524, 923
  • Neufeld et al. (2006) Neufeld D. A. et al., 2006, ApJ, 649, 816
  • Ostriker, Stone & Gammie (2001) Ostriker E. C., Stone J. M., Gammie C. F., 2001, ApJ, 546, 980
  • Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
  • Pagani, Roueff & Lesaffre (2011) Pagani L., Roueff E., Lesaffre P., 2011, ApJ, 739, L35
  • Papadopoulos (2010) Papadopoulos P. P., 2010, ApJ, 720, 226
  • Pollack et al. (1994) Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., Fong W., 1994, ApJ, 421, 615
  • Price, Federrath & Brunt (2011) Price D. J., Federrath C., Brunt C. M., 2011, ApJ, 727, L21
  • Schöier et al. (2005) Schöier F. L., van der Tak F. F. S., van Dishoeck E. F., Black J. H., 2005, A&A, 432, 369
  • Schroder et al. (1991) Schroder K., Staemmler V., Smith M. D., Flower D. R., Jaquet R., 1991, Journal of Physics B Atomic Molecular Physics, 24, 2487
  • Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
  • Shetty et al. (2011a) Shetty R., Glover S. C., Dullemond C. P., Klessen R. S., 2011a, MNRAS, 412, 1686
  • Shetty et al. (2011b) Shetty R., Glover S. C., Dullemond C. P., Ostriker E. C., Harris A. I., Klessen R. S., 2011b, MNRAS, 415, 3253
  • Sohn et al. (2007) Sohn J., Lee C. W., Park Y.-S., Lee H. M., Myers P. C., Lee Y., 2007, ApJ, 664, 928
  • Staemmler & Flower (1991) Staemmler V., Flower D. R., 1991, Journal of Physics B Atomic Molecular Physics, 24, 2343
  • Steinacker et al. (2003) Steinacker J., Henning T., Bacmann A., Semenov D., 2003, A&A, 401, 405
  • Tielens & Hollenbach (1985) Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 722
  • Tomida et al. (2013) Tomida K., Tomisaka K., Matsumoto T., Hori Y., Okuzumi S., Machida M. N., Saigo K., 2013, ApJ, 763, 6
  • Turner et al. (1992) Turner B. E., Chan K.-W., Green S., Lubowich D. A., 1992, ApJ, 399, 114
  • van der Tak et al. (2007) van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007, A&A, 468, 627
  • van Dishoeck & Black (1988) van Dishoeck E. F., Black J. H., 1988, ApJ, 334, 771
  • Watson (1972) Watson W. D., 1972, ApJ, 176, 103
  • Wolfire, Hollenbach & McKee (2010) Wolfire M. G., Hollenbach D., McKee C. F., 2010, ApJ, 716, 1191
  • Wolfire et al. (1995) Wolfire M. G., Hollenbach D., McKee C. F., Tielens A. G. G. M., Bakes E. L. O., 1995, ApJ, 443, 152
  • Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278
  • Yang et al. (2010) Yang B., Stancil P. C., Balakrishnan N., Forrey R. C., 2010, ApJ, 718, 1062

Appendix A Specific Heats

Calculating the time evolution of the temperature requires knowing the specific heat per H nucleus at constant volume , defined by


where is the gas internal energy per unit volume, given by


where the sum runs over species , is the number density of species , and is the partition function per unit volume for that species. The latter is given by


where the terms appearing in the equation above are the partition functions for the translation, rotational, vibrational, and spin degrees of freedom of species . In principle we should also include a term describing electronic degrees of freedom, but at the relatively low temperatures for which DESPOTIC is intended, we can safely assume that these are not excited. For all the species included in DESPOTIC except molecular hydrogen (i.e. for H i, He, H, and ), the contribution of the specific heat is trivial, because all of the partition functions except translation and spin are unity, and the spin term is temperature-independent. Thus for all these species


For ortho- and para-H on the other hand, and are not unity (Black & Bodenheimer, 1975; Boley et al., 2007; Tomida et al., 2013):


where K and K. Note that the vibrational partition function is the same for ortho- and para-H, but the rotational partition functions are different. With these partition functions, the energy per unit volume including all species is


where again the sum runs over all all species.

Deriving the specific heat from this expression requires making an assumption about how the number densities of ortho- and para-H vary with temperature. At the low temperatures found in interstellar clouds, there is generally no efficient mechanism for converting between the two states, and thus the most reasonable assumption is that these number densities are temperature-independent. Observations showing that the ortho- to para- ratio in molecular clouds is far from equilibrium (e.g. Neufeld et al., 2006; Pagani, Roueff & Lesaffre, 2011; Dislaire et al., 2012) support this assumption. For temperature-independent values of and , we therefore have


Note that this expression involves the abundances ratios rather than number densities because we have normalized all quantities to the number density of H nuclei. The specific heat at constant pressure is simply .

Appendix B Heating and Cooling Processes

Here I give explicit formulae for all the heating processes included in DESPOTIC. In the following description, all heating, cooling, and energy exchange rates are given as energies per H nucleus per unit time.

b.1 Gas

b.1.1 Ionization Heating

Gas can gain energy through ionization heating; in this process primary electrons with energies produced when the gas is ionized by cosmic rays or hard x-rays thermalize, adding energy. The rate at which this process adds energy is given by


where is the energy added per primary ionization. The value of in turn depends on the bulk chemical composition of the gas, which determines how much of a primary electron’s eV of energy is lost via radiation rather than transformed into heat. This problem has been discussed by a number of authors (Dalgarno & McCray, 1972; Glassgold & Langer, 1973; Wolfire et al., 1995; Dalgarno, Yan & Liu, 1999; Wolfire, Hollenbach & McKee, 2010; Glassgold, Galli & Padovani, 2012). In predominantly atomic regions, the main pathway to thermalization is Coulomb scattering of the primary electron off other free electrons, and collisional excitation of H and He by the primary electron followed by collisional de-excitation of the excited atom. In this regime DESPOTIC uses the approximation recommended by Draine (2011),


In molecular regions the situation is far more complicated due to the additional thermalization channels provided by excitation of the rotational and vibrational levels of H (followed by collisional de-excitation), by dissociation of H, and by chemical heating, in which primary electrons produce reactive ions such as H, H, and He that subsequently undergo exothermic reactions with neutrals such as CO, HO, and O. In this case becomes a complex function of the gas density and temperature, and the abundances of various species, and ranges from eV as these quantities change (Glassgold, Galli & Padovani, 2012). Given the complexity of the problem, and the level of inaccuracy inherent in any one-zone model, DESPOTIC relies on a simple piecewise fit to the numerical results of Glassgold, Galli & Padovani (2012, their Table 6) on the density-dependence of in molecular regions:


where the values of in the above expression are in units of cm.

To handle the case where the composition includes both molecular and atomic gas, DESPOTIC assumes that the atomic and molecular regions are physically separated (which, depending on the physical situation, may or may not be a good assumption). In this case the total heating rate can be computed simply by summing the heating rates in the atomic- and molecular-dominated regions, weighted by their number fractions:


b.1.2 Photoelectric Heating

Gas can also gain energy through grain photoelectric heating, whereby a primary electron ejected from a dust grain by a far-ultraviolet (FUV) photon thermalizes with the gas. Unlike cosmic rays, the FUV photons responsible for photoelectric heating can be attenuated by dust rather easily, and the photoelectric heating rate therefore depends on four factors: the strength of the ISRF, the abundance of dust grains, the amount of dust shielding, and the energy yield per photoelectron; as with cosmic ray heating, the latter value has been estimated by numerous authors (Watson, 1972; de Jong, 1977; Tielens & Hollenbach, 1985; Bakes & Tielens, 1994; Wolfire et al., 2003). To account for dust shielding, which obviously varies from point to point within a real cloud, DESPOTIC uses the simple approximation proposed by Krumholz, Leroy & McKee (2011), whereby the eV photons responsible for photoelectron production are considered to be attenuated by half the mean extinction of the cloud. Since the dust opacity is relatively flat across this energy range ( variation in the models of Draine 2003), we can assign a single cross section , which is cm H for Milky Way dust. This value is near the middle of the range found in the models of Draine (2003). With this approximation, the photoelectric heating rate becomes


b.1.3 Gravitational Heating

A third possible source of heating is adiabatic compression. This obviously depends on the hydrodynamics of the flow, something that is not naturally included in a one-zone model like that used in DESPOTIC. However, this effect is calculable in the special case of compression due to gravitational contraction, as in protostellar cores for example. In this case the heating rate may be computed using the approximation introduced by Masunaga, Miyama & Inutsuka (1998),


where is a dimensionless constant of order unity that depends on the nature of the gravitational collapse. From their numerical calculations, Masunaga, Miyama & Inutsuka (1998) find . Since in general most interstellar clouds are not in a state of collapse, by default DESPOTIC does not include gravitational contraction heating, and sets . However, users do have the option of overriding this default.

b.1.4 Line Cooling

The primary cooling mechanism for gas is line radiation. For each emitter species , there is a rate of line cooling , so that the total line cooling rate is


I defer a calculation of to § B.3.

b.1.5 Dust-Gas Energy Exchange

Finally, gas can either heat or cool by exchanging energy with the dust via collisions. The gas-dust energy exchange rate is given by


where is the grain-gas coupling coefficient and the sign convention is that positive values correspond to heating of the gas and cooling of the dust. Note the presence of the clumping factor , since this is a collisional process. The coupling constant depends on the grain abundance, chemical composition, size distribution, and charge state. For Milky Way dust, Goldsmith (2001) recommends a value erg cm K for H-dominated regions, and Krumholz, Leroy & McKee (2011) estimate a value of erg cm K for H i-dominated ones, with the difference arising due to the change in both the number and mean mass of free particles between H i and H-dominated regions.

b.2 Dust

b.2.1 Cooling by Thermal Radiation

Dust grains can lose energy via thermal continuum radiation. To compute the cooling rate, consider a population of spherical grains with distribution of radii given by , where we normalize the distribution function such that is the total number density of dust grains. Let be the absorption efficiency for absorption of radiation of frequency , so that the cross section of the grain to radiation of frequency is . Further let be the Planck-weighted mean efficiency, where is the Planck function evaluated at temperature . Given this definition, we can write the rate of thermal radiation cooling from dust grains of temperature as


we have defined the term in square brackets to be the mean dust cross section per H nucleus . This expression assumes that the cloud is optically thin to its own cooling radiation; we treat the optically thick regime below. We approximate that will vary as a powerlaw with , and we therefore write


For Milky Way dust, typical opacities are cm H (Pollack et al., 1994; Semenov et al., 2003), and for temperatures such at cm is much larger than the typical grain size, we expect ; detailed grain models show that this expectation holds up to K (Semenov et al., 2003). DESPOTIC leaves both and as user-settable parameters. A naive expectation is that, at sub-Solar metallicities, , where is the dust abundance relative to Solar.

The above estimate is valid only as long as the cloud is optically thin to its own cooling radiation, which is true only as long as . Given the small value of for Milky Way dust, departures from the optically thin regime do not begin until extremely high column densities. However, there are circumstances, for example in the molecular clouds of starburst galaxies, where and can be high enough to render the optical depth to cooling radiation large. A truly accurate calculation of the cooling rate in this regime requires a multi-zone numerical treatment with a radiative transfer code such as dusty (Ivezic & Elitzur, 1997) or SteinRay (Steinacker et al., 2003), or a sophisticated analytic approximation (e.g. Chakrabarti & McKee, 2005). However, we can obtain a very crude treatment of the optically thick regime by noting that the maximum possible cooling rate for the cloud is simply , the blackbody rate for a sphere of radius equal to the cloud radius. Rewriting this as a rate per H nucleus, the maximum possible dust cooling rate is


DESPOTIC adopts the approximation


b.2.2 ISRF Heating

Grains can be heated by absorbing the interstellar radiation field produced by stars. To compute the rate of dust heating from the ISRF, we must perform a calculation similar to that for . In analogy to , we define , where is the energy density of the ISRF at frequency , as the ISRF-averaged absorption efficiency. In general . Thus, unlike in the case of thermal cooling where optical depth effects are important only in extreme circumstances, attenuation of the ISRF will be important even at modest column densities. As with photoelectric heating, it is clear that there is no single value that describes the rate of dust heating within an optically thick cloud; heating rates will be high at the edge and low at the center. Moreover, unlike in the case of photoelectric heating, the range of photon energies responsible for heating is quite broad, with half the heating coming from photons with wavelengths m even for the unattenuated ISRF (B. Draine, 2013, priv. comm.). As a result, the spectrum of the heating field changes as one moves into a cloud and shorter wavelength photons are selectively attenuated. Consequently, in addition to the geometric uncertainty, there is an additional one in the choice of dust cross section to assign. In order to maintain simplicity, DESPOTIC does not attempt to treat this problem in detail, but instead uses the same approximation as for photoelectric heating, i.e. that the characteristic heating rate is to be computed assuming an attenuation equal to half the mean value for the cloud, using a single grain cross section to compute the attenuation. With this approximation, the heating rate of grains due to the ISRF is


where is the dust abundance relative to the Milky Way value, is the energy density of the ISRF, is the energy density for the Milky Way’s ISRF, is the cross section we assign for ISRF attenuation, and the numerical coefficient is taken from Goldsmith (2001). The choice of is somewhat difficult for the reasons stated above, and if very high accuracy is desired it should be computed on a case-by-case basis. However, a reasonable default for Milky Way dust is cm H, which is roughly halfway between the values appropriate for the unextincted ISRF and the value expected for an ISRF extincted by an optical depth of 2 in V band (B. Draine, 2013, priv. comm.).

It is worth noting that, because the ISRF is exponentially attenuated by dust, when we are likely to find that is negligibly small even when is very large. In this circumstance, the ISRF is so thoroughly attenuated that none of it reaches the cloud interior where we are computing the temperature. However, if this happens, the hot outer parts of the cloud that are directly exposed to the ISRF will heat up and generate a background infrared field within the cloud interior. If the cloud is optically thin to IR cooling radiation the intensity of this field will be low and it can be neglected as a heat source. If the cloud is optically thick to IR, on the other hand, the background IR field will build up, and will heat the cloud interior. DESPOTIC provides a mechanism to handle this phenomenon by including an infrared radiation field (see the following section), and in circumstances where ISRF heating is negligible, heating by the infrared radiation field should take its place. As for the case of the cooling rate when the cloud is optically thick to IR, calculating the intensity of the background field in this circumstance requires a more sophisticated model than the one-zone treatment that DESPOTIC provides. However, we can solve the limiting case of an extremely optically thick cloud subject to external heating. If such a cloud absorbs all of the background ISRF incident on its surface, the total heating rate is , and the heating rate per H nucleus is


DESPOTIC adopts the approximation


Equating this with the limiting cooling rate for an extremely opaque cloud, , gives an equilibrium temperature for both the dust and the infrared radiation field


i.e. the dust and IR radiation field within the cloud reach a temperature such that the radiation energy density within the cloud is equal to the ISRF energy density outside it, as expected for a blackbody.

b.2.3 Heating by Infrared Radiation and the CMB

The final source of radiative energy for dust is the background thermal radiation field, and the CMB. Since both of these sources of radiation are thermal, they may be handled using exactly the same mechanics as thermal radiative cooling. The heating rate is therefore


b.2.4 Line Heating

In addition to emission and absorption of continuum radiation, there are two additional processes that can heat and cool dust grains. The first of these, collisional exchange with the gas, is discussed in § B.1. The other is absorption of line photons emitted by the gas. If we let


be the population-averaged grain cross section per H nucleus at frequency , then the mean optical depth of the cloud to line photons at frequency is . In principle one could use a detailed grain model to obtain at the frequencies of all the relevant lines. However, this procedure would be cumbersome, and is likely unimportant for most clouds since, not surprisingly, both cooling radiation and observable emission tend to be dominated by lines at frequencies such at clouds are optically thin. Nonetheless, to approximate the effects of clouds becoming optically thick to line radiation, DESPOTIC approximates by , where is the line frequency, and 208 GHz is multiplied by 10 K. With this approximation, and using the same expression for the line photon escape probability versus optical depth as discussed below in § B.3, we obtain the final heating rate of the dust due to absorption of line photons:


where is the frequency of the line produced by atoms / molecules of species transitioning between states and (see § B.3), is the escape probability for a photon corresponding to line computed using the dust optical depth, and the sum runs over all species and level pairs .

b.3 Level Populations and Line Radiation

b.3.1 Level Populations in Optically Thin Clouds

Calculating the line cooling rate requires determining the level populations for all emitting species. Consider an emitting species , and let be the energy of the th quantum state of that species, where the states are numbered by energy so that for all states . The degeneracy of state is , and the Einstein coefficient describing the rate of spontaneous radiative transitions from state to state is , where for . Finally, let be the rate coefficient for collisional transitions from state to state induced by collisions with some collision partner ; the upward and downward rate coefficients obey the usual relationship , where and By convention for .

For our species of interest, we wish to solve for the fraction of atoms / molecules in state , when that species is mixed with a gas of a given bulk composition, number density , and gas temperature , and the cloud is immersed in a sea of cosmic microwave background photons. If the cloud is optically thin to photons at the frequencies of the lines connecting the various states, in statistical equilibrium the various level populations are determined implicitly by the conditions that the rate of transitions into and out of each level balance:




are the photon occupation number at the frequency of the line connecting states and ,777Naively one would think that, in a cloud that builds up a significant trapped infrared radiation field, then the photon occupation number should also include a contribution from this field, of the same form as equation (44) but with replaced by . However, this is often not the case, for the following reason. Even in high column density environments where a significant dust-trapped infrared radiation field builds up, the spectrum of this radiation field is often not Planckian at low frequencies. This is because the dust opacity generally falls as at low frequencies, and so even if the dust is opaque to radiation near the peak of the spectral energy distribution, it is usually transparent at low frequencies. This results in a radiation spectrum that is Planckian at higher frequencies but very sub-Planckian at low frequencies, and thus has a much lower photon occupation number that a true blackbody like the CMB. A fully accurate calculation of level populations would account for this effect by solving for the frequency-dependent dust-mediated radiation field and using the appropriate photon occupation number to calculate the level populations. However, as noted above, it is not feasible to determine the dust radiation field accurately in a one-zone model. I therefore choose to optimize the accuracy of DESPOTIC for the case of lines at frequencies where the dust is optically thin, since these are, obviously, the lines that are most important for both cooling and observation. This choice dictates that the dust radiation field be ignored when computing the level populations, on the basis that its photon occupation number will be small. However, this choice does limit the accuracy of DESPOTIC for lines where infrared pumping is important, as discussed in more detail in § 5. and the rate of collisional transitions between states and summed over all collision partners . Here is the abundance of a given collision partner relative to , and the collision partners considered by DESPOTIC are H, He, pH, oH, , and H. As usual, collision rates are multiplied by the clumping factor . The left-hand sides of equations (43) describe the rate of transitions into state from all other states , with the first term representing the rate of collisional transitions, the second representing the rate of radiative transitions (including both spontaneous and stimulated emission), and the third term describing the rate of absorptions. The right-hand sides represent the rate of transitions from state to all other states , with the three terms again representing collisional transitions, spontaneous and stimulated emission, and absorption. These equations are supplemented by the constraint equation


and together equations (43) and (46) constitute a complete system.

For computational purposes it is convenient to rewrite this system as a matrix equation. Consider a species for which we track distinct energy levels. With some manipulation, equations (43) and (46) may be rewritten as888Note that DESPOTIC does not use the standard procedure in the stellar atmospheres community of recasting the equations in terms of departure coefficients.This choice is motivated by the fact that, for most of the calculations for which DESPOTIC is intended, most of the states of most species will be very far from LTE. This vitiates any advantage to recasting the equations in departure coefficient form.


where M is an matrix whose elements are