---
author:
- |
    G. J. Ferland, R. L. Porter, P. A. M. van Hoof, R. J. R. Williams, N. P. Abel,\
    M. L. Lykins, Gargi Shaw, W. J. Henney, P. C. Stancil
bibliography:
- 'LocalBibliography.bib'
- './common/bibliography2.bib'
title: The 2013 Release of Cloudy
---

=1

Introduction
============

Most quantitative information we have about the cosmos comes from spectroscopy. In many astronomical environments the density is too low for equilibrium thermodynamics to apply, so the ionization, molecular state, level populations, kinetic temperature, and the resulting spectrum are the result of a host of microphysical processes. As a result the spectrum reveals much about the properties of an object, but it also means that modeling this detail is complex. Analytical results are only possible in certain limits, so numerical simulations must be used. Texts that review this field include @Spitzer1978, @DopitaSutherland03, @Tielens2005, @AGN3 (hereafter AGN3), and @Draine.B11Physics-of-the-Interstellar-and-Intergalactic-Medium.

[^1] is an open source plasma simulation code that is designed to simulate conditions in a non-equilibrium gas, and predict its spectrum. The code incorporates physical processes from first principles, as much as possible. The goal is to simulate the ionization, level populations, molecular state, and thermal state, over all extremes of density and temperature. Our approach, working from fundamental processes, means that  can be applied to such diverse regions as the corona of a star, the intergalactic medium, or the accretion disk near the supermassive black hole in a luminous quasar. As a result, the code is widely used, with nearly 200 papers citing its documentation each year. The diversity of problems it can address is a testimonial to the importance of treating the atomic physics at an elementary level.

Processor power has always limited our ability to simulate detailed microphysics. Improved computers and advances in atomic and molecular physics allow a better simulation. Improved numerical methods or coding techniques make the solutions more robust. These advancements to the fidelity of the simulation improve our insight into the inner workings of astronomical objects. Because of these changes, , like most software, goes through a development / release cycle. Our goal is to make major updates every two or three years.

This paper is a progress report on the improvements to  since the last major review, @Ferland1998, referred to as F98 in the following. Although the code’s download includes extensive documentation that is continuously updated, there has not been a recent major review. We rectify that need here.

’s development leading up to the 1998 review had emphasized the UV, optical, and IR spectra of ionized gas. The simulations have been extended to fully molecular regions with predictions of the associated IR and radio emission since then. The code can handle a broad range of physical state, from predominantly molecular to fully ionized, a broad range of densities from the low density limit to roughly $10^{15} \pcc$, and temperatures ranging from the CMB to $10^{10} \K$.

In the present paper, we summarize the major advances in the code since F98. Most of this work has been documented in past papers, which we cite. In the interests of brevity we only give references to the relevant papers with a brief summary of the advances they document. We discuss some technical details about the code, its operation, and installation. We present several calculations that show the range of applicability of . These include the properties of the ionized, atomic, and molecular gas produced by the radiation field of an Active Galactic Nucleus (AGN) and a demonstration of the range of particle and radiation density that can be considered. We conclude with an outline of future development directions.

The physics of ions, molecules, and grains
==========================================

 was originally designed to simulate the dense gas found near the black hole and accretion disk in The so-called “Broad Line Region” (BLR) mainly emits in the UV and optical, and these spectral regimes were the original focus. Some investigations are @Rees1989, @FerlandPeterson1992, and @Baldwin1995. Lower density gas, at larger distances from the center, produces the “Narrow Line Region” (NLR) spectrum, which includes a range of forbidden and permitted lines in the UV, optical, and IR. The study by @Ferland1983 is an example. Such investigations drove the development of , as summarized in the 1998 review, although Photo-dissociation Region (“PDRs”) and X-ray Dissociation Regions (“XDRs”) were also simulated [@Ferland1994; @FerlandFabian2002].

The following sections describe our improvements in the treatment of ions, molecules, and grains since F98, by citing those papers which introduced the advances. This is not meant to be a comprehensive review of the literature, rather, only a description of advances to  since the last review.

Structure of the H-like and He-like iso-electronic sequences
------------------------------------------------------------

Hydrogen and helium are the most common elements in the universe and, as a result, need to be treated with the greatest precision. Their structure is different from most heavy elements, with a first excited level at (1-1/$n^2$) = 0.75 of the ionization energy and more highly excited levels close to the ionization limit. By contrast, the heavy elements can have many low-lying levels. This means that, for most temperatures, lines from H and He – like species will have a strong recombination component, while lines of many-electron systems will be predominantly collisionally excited. The structure of the H and He – like iso sequences also means that a significant number of excited states are needed for the model atoms to correctly go to LTE in the high particle or photon limits.

Steve Cota developed the original model of , , and  emission in   as part of his PhD thesis [@1987PhDT.........7C]. He also developed an approximate treatment of three-body recombination for the heavy elements, described in the next section. This was extended by Jason Ferguson in his PhD thesis [@Ferguson.J97Spectral-Simulations-and-Abundance-Determinations; @FergusonFerland1997] to include more levels, as processor power increased. These models, which involved 15 levels with a number of higher pseudo states, were significant time sinks on the computers at that time.

Today’s unified model of the H and He iso-electronic sequences was developed as part of Ryan Porter’s thesis. Because the high charge states occur in hot gas and the line energies scale as $Z^2$, these iso-sequences sort themselves into two spectral regions. , , and   produce strong lines in the optical while , , , , , , etc produce lines in the X-ray. Despite these differences, the physics has many similarities. These sequences are treated with a common code base, which results in greater simplicity and reliability.

@Porter2005 and @Bauman2005 describe the model of  emission. @PorterFerland2007 describe ions of the He sequence, which mainly emits in the X-Ray. @LuridianaEtAl09 summarize expansions to the H-like sequence while @PorterFerlandMacAdam2007, @Porter.R09Uncertainties-in-theoretical-HeI-emissivities:-HII-regions, and @Porter.R12Improved-He-I-emissivities-in-the-case-B-approximation discuss uncertainties and more recent improvements in the atomic data for .

A schematic representation of one of the elements of the H-sequence is shown in Figure \[fig:LevelsResolvedCollapsed\]. The He-sequence has similar structure except that it is resolved into singlets and triplets with twice as many levels. Principal quantum number increases upward with the continuum at the top and $nl$ terms are indicated from left to right. Lower $n$ configurations are resolved into $nl$ terms, the “resolved” levels. Above a certain quantum number $l$-changing collisions become fast enough to guarantee that $l$ terms are populated according to their statistical weight within the $n$ configuration [@PengellySeaton1964]. Such higher levels are treated as “collapsed levels” which are $nl$ mixed.

![The iso-sequence atomic level structure consists of low-lying $nl$ resolved ($nls$ for He-like) terms with $l-$mixed ($ls$ for He-like) collapsed configurations. The number of resolved and collapsed levels can be specified when the simulation parameters are established. \[fig:LevelsResolvedCollapsed\]](LevelsResolvedCollapsed){width="\linewidth"}

This treatment includes ions of all elements up to Zinc. The models include photoionization / recombination, collisional ionization / 3-body recombination, to all levels, and collisional and radiative processes between levels, so behave correctly in the low density limit and go to LTE at high densities or exposed to a true blackbody radiation field [@Ferland1988]. Line trapping, collisions, continuum lowering, and absorption of photons by continuous opacities, are all included as general processes [@Rees1989].

The user can adjust the number of resolved and collapsed levels modeled when the simulation is specified. The spectrum is predicted with great precision, an accuracy of better than 1%, when a larger number of levels are used [@Porter.R12Improved-He-I-emissivities-in-the-case-B-approximation]. This comes at the cost of increased execution times. Smaller models are often used for simulations of clouds with significant column densities due to the compute time required. The default treatment includes the greatest number of levels for H, and increasingly smaller numbers of levels for He, common second row elements like C and O, Fe, and the remaining low abundance elements.

Structure of other ions
-----------------------

 includes all ions of the lightest thirty elements. @Lykins.M12Radiative-cooling-in-collisionally-and-photo summarize recent developments in the treatment of the ions that are not part of the H and He iso-sequences.

An equivalent two-level system is assumed in modeling the ionization balance of ions of the Li-like and multi-electron iso sequences. Photoionization and collisional ionization from the ground configuration is balanced by recombinations to all levels. This assumes that nearly all populations are in ground, an approximation which is valid for moderate densities and the low temperatures usually found in photoionization equilibrium. Photoionization cross sections are from the Verner database [@VernerFerlandKorista1996] while radiative and dielectronic recombination rates are from the Badnell web site[^2] as described in @BadnellEtAl03 and @Badnell06, supplemented by data calculated as described by @VernerFerland1996. Charge exchange ionization and recombination rates are taken from an updated version of the @Kingdon1996 database.

The treatment of inner-shell processes, including line emission following removal of an inner-shell electron, largely follows F98. As described below, the treatment of multiple electron ejection, which had followed @Weisheit1972, is now generalized to non-adjacent stages of ionization [@HenneyEtAl05]. Ionization and recombination coupling non-adjacent ion states can also be important in grain surface recombination.

This treatment is approximate at high densities and for temperatures which are a significant fraction of the ionization potential of a species. The use of summed recombination rate coefficients in effect assumes that all recombinations eventually populate the ground level. Ionization processes out of excited levels are neglected. Both are no longer true if the density is large enough for excited levels to play an important role. Expanding the treatment of these species to the approach now used for the H and He iso sequences is a high priority for future development.

Bound levels of each ion are treated with a variety of models. We are enhancing the code to use external databases as much as possible. We have the ability to read the Chianti atomic database[^3] [@Dere.K97CHIANTI---an-atomic-database-for-emission; @Landi2012] and this release of  includes Version 7.0. Chianti is used for most ions of Fe while for other species we mainly use our internally developed atomic database. @Lykins.M12Radiative-cooling-in-collisionally-and-photo shows that our original database, which is embedded in the C++ source, is in good agreement with Chianti.

Additionally we are starting to develop our own database, “Stout”. We will add new models of ions or molecules to this second database, and make it publicly available along with .

Molecular chemistry {#molchem}
-------------------

 initially included the chemistry network described by @Black1978 which was expanded to PDRs and XDRs as described by @Ferland1994. Nick Abel a massive upgrade to the heavy-element chemistry network as part of his PhD thesis, described in @Abel2004. Later refinements are discussed in @Abel2005, @Shaw2005, and @Shaw2006. Appendix A of @Abel2005 gives details of the numerical methods along with differences between UMIST and  reaction rates.  had predicted column densities for about 20 heavy element molecules, consisting of C and O atoms. It could not calculate physical conditions deep in a PDR or a molecular cloud, where most gas phase C, N, and O is in the form of molecules, due to numerical instabilities in the chemistry solver then used. The upgraded chemistry solver has no restrictions, as described in sections below.  now calculates the chemical abundance of 83 molecules using a network including $\sim 10^3$ chemical reactions involving molecules containing H, He, C, N, O, Si, S, and Cl atoms. The network adjusts automatically when elements or species are disabled.

Most reaction rates come from the UMIST 2000 database [@Le-Teuff.Y00The-UMIST-database-for-astrochemistry-1999] as updated for the Leiden workshop and described by @Roellig2007. We also predict the freeze-out of , , and  on grains, using the data given in @Hasegawa1993. Both time-steady and time-dependent chemical evolution calculations are possible.

The effects of cosmic rays and suprathermal secondary electrons can be very important in molecular regions [@1968ApJ...152..971S] and are treated as described by @FerlandMushotzky1984 and @FerlandFabianEtAl2009. Non-thermal particles can both excite and dissociate the gas. @Ferland1994 and @FerlandFabian2002 describe the treatment of chemistry in an X-ray dominated filaments in cool-core galaxy clusters. @Shaw2008 describe an application to the local ISM.

The treatment of  is fully self consistent with its surroundings [@Shaw2005]. Formation on grain surfaces is treated using the derived grain properties and the @Cazaux.S04H2-Formation-on-Grain-Surfaces catalysis rates. Destruction by photoexcitation in electronic states is treated by computing the radiation field at each point. This includes line self shielding, the attenuation by continuous opacity sources such as grains, photoelectric opacity of the heavy elements, and Rayleigh scattering, and emission by the cloud itself.  , the resonance lines, and electronic lines of  are especially important sources of higher-energy photons.

Molecular excitations and spectra
---------------------------------

We next describe our treatment of molecular emission processes. Gargi Shaw developed our model of the  molecule, the most common molecule in the Universe, as part of her PhD thesis [@Shaw2005]. All levels within the ground electronic state are included [@Dabrowski.I84The-Lyman-and-Werner-bands-of-H2], along with all electronic excited states on the Meudon web site[^4] and @Abgrall.H94The-Bprime1Sigmau-rarr-X1Sigmag-and-D1Piu. Collision rates are from @Wrathmall2007. @Lee.T05Close-coupling-calculations-of-low-energy-inelastic [@Lee.T06State-to-state-rotational-transitions-in-H2H2; @LeeH2H22008], @Shaw.G09Rotationally-Warm-Molecular-Hydrogen, and @Gay.C12Rovibrationally-Resolved-Direct-Photodissociation summarize more recent updates.

Grain formation pumping of  can be treated using several theories. @Shaw2005 provides more details. We use the @Takahashi2001 results by default.

Line emission for molecules heavier than  is predicted using the level energy, collision, and radiative data in the LAMDA database [@Schoier.F05An-atomic-and-molecular-database-for-analysis]. The molecular excitation are solved simultaneously and self-consistently with the global environment. This includes grain properties, including emission, so molecular pumping in the infrared continuum is automatically included (if the pumping lines are present in the molecular data), along with attenuation by grains and other continuous opacity sources. Line optical depths are computed for each point in the cloud, and the full radiative transfer performed as described below.

Grains
------

Grains have been included in  since the very beginning [@Martin1980], and the basic physical processes are summarized in Appendix C of @BFM. @VanHoof2004 describe the main improvements made since @BFM.

The grain model in  includes all relevant processes: absorption and scattering of light including stochastic heating effects, the photoelectric effect including Auger emissions in X-ray environments, collisional charging (electrons as well as atomic ions), thermionic emissions, collisional energy exchange between the grains and the gas, and a calculation of the grain drift velocity. We also consider molecular freeze-out on grain surfaces and grain surface reactions. These are discussed in Section \[molchem\].

The local radiation field, including the attenuated incident Spectral Eenergy Distribution (SED) and the spectrum emitted by the cloud (generally  is the most important) are all treated self-consistently, with grains providing the opacity affecting the light, and the light affecting the grain properties following absorption. The code fully treats stochastic heating of grains using a robust and efficient algorithm (which is a comprehensively upgraded version of a code originally written by K. Volk), implementing an improved version of the procedure described in @Guhathakurta1989. This includes an approximate treatment for stochastic heating by particle collisions. Combined with resolved size distributions, this will lead to a much more realistic modeling of the grain emission under all circumstances. The stochastic heating code needs a grain enthalpy function. Supported functions are taken from @Guhathakurta1989, @Dwek1997, and @Draine2001.

The grain physics includes the advances discussed by @Weingartner2001a and @Weingartner2006. The latter provides a realistic model for grains in X-ray environments, including Auger emissions by the grains. The grain model has a detailed treatment of the photoelectric effect and collisional processes, and includes thermionic emissions. The charge for each grain constituent is determined self consistently with the local radiation field and gas properties, using the hybrid grain charge model described by @VanHoof2001 and @VanHoof2004. The grain charge is included in the overall charge balance of the system, which has a significant effect on the modeling results for molecular regions [@Abel2008]. Charge exchange on grain surfaces is treated following @Draine1987 assuming that electrons will be exchanged between the grain and the colliding particle until a minimum energy state is reached. The grain drift velocity w.r.t. the gas is calculated using the theory from @Draine1979.

Extensive comparisons in collaboration with Joe Weingartner done in 2001 show that the photoelectric heating rates and collisional cooling rates predicted by  agree very well with the results from the @Weingartner2001a model for a wide range of grain sizes (between 5 Å and 0.1 $\mu$m), and using various choices for the incident radiation field. A detailed discussion of this comparison can be found in @VanHoof2001.

We include an embedded model for Mie scattering [@Mie1908], which uses refractive index data for each material type. It is based on the code described in @Hansen1974 and was later modified by P.G. Martin. New grain materials can be included by specifying the relevant refractive index data or by supplying opacity tables. The user specifies the size-distribution function and the number of grain size bins. A number of size-distribution functions are available as built-in functions, including the ones described in @Mathis1977 [hereafter MRN], @BFM, and @Abel2008.  will then “compile” the grain data to create scattering and absorption cross sections, the scattering asymmetry factor, and the inverse attenuation length for each material type and grain size as a function of frequency. A number of pre-compiled grain types are included in the code distribution.

In practice, it is difficult to introduce new grain materials because of the wavelength range considered by the code. Laboratory data can extend from the IR into the UV, but little data exists shortward of $\sim 0.1 \mu$m. Theoretical relationships between the bulk grain properties and the photoionization cross sections of the constituent atoms can be used to create complete refractive index data for grain materials when combined with the Kramers-Kronig relations. Refractive index data for astronomical silicate, amorphous carbon, and graphite, based on @Martin1991, @Rouleau1991, and @Laor1993 are included in the  download. Opacity data for PAHs from Volk (private communication, based on data from [@Bregman1989], [@Desert1990], and [@Schutte1993]), @Li2001, and @Draine2007 are built into the code.

The Mie code includes the possibility to mimic mixtures of materials in grains using effective medium theory (EMT). The following EMT recipes are supported by : @Bruggeman1935, @Stognienko1995, and @Voshchinnikov1999 (based on the theory in [@Farafonov2000]). The first two are appropriate for randomly mixed grain materials, while the latter is intended for layered grains.  includes a refractive index file for vacuum which enables modeling fluffy grains when combined with other materials using an EMT.

@Abel2008 discusses how our treatment, which is based on elementary processes as as possible, affects results for PDR simulations. Conventional PDR codes use precomputed integrals of photoionization or photodissociation cross sections over an SED representative of the Galactic starlight background. This rate is assumed to be a function of the radiation field scaled to the intensity of the Galactic background, and the visual extinction . In contrast we explicitly integrate stored cross sections over the local radiation field for those processes which have energy-specific data. These include photo rates for all atoms, ions, and grains, and some molecules, most notably . The incident radiation field is attenuated by the computed gas and dust opacity. Our treatment can handle any grain opacity distribution or abundance, or SED shape.

The grain scattering theory predicts the scattering asymmetry factor $g$, which is the average of the cosine of the scattering angle of an incident photon [@Martin1979]. When a point source like a star is viewed through an intervening cloud at a large distance from the star, photons scattered by even a small amount are lost from the beam and cannot reach the observer. So in this case scattering attenuates the radiation field by the full opacity, $\alpha_{scat}$. We refer to this, the quantity measured by observations of stellar extinction, as the point-source extinction. On the other hand, when modeling a cloud irradiated by a star, photons that are scattered in a forward direction can still propagate into the next radial zone of the model and are therefore not lost from the radiation field [@Martin1979; @BFM; @AGN3]. In  this is approximated using an effective opacity $\alpha_{scat} (1-g)$, which is said to discount forward scattering. We refer to this as the extended extinction. Both extended and point-source extinctions are reported by the code, but the user should be aware that the PDR literature always uses the point-source extinction.

Line transfer
-------------

The treatment of line transfer is largely unchanged from F98. The escape probability formalism is used, using the framework given by, among others, @Irons.F78On-the-equality-of-the-mean-escape-probability, @FerlandElitzur1984, @Rybicki.G84Escape-probability-methods-pp-21-64, @Kalkofen.W84Methods-in-radiative-transfer, @Netzer1985, @Elitzur.M85Line-fluorescence-in-astrophysics, @Elitzur1986, @Kalkofen1987, @Ferland1992, and @Elitzur1992.

All permitted, and many forbidden, lines are transferred using a common approach. With this, processes such as line trapping and thermalization, pumping by the local radiation field, line destruction by background opacities such as photoelectric or grain absorption, are included for the $\sim 10^5 - 10^6$ atomic and molecular lines considered by the code. Line processes couple into the model atoms, which can go over to the correct thermodynamic limits when exposed to a blackbody [@Ferland1988]. The simulation is done self consistently, with feedback between various constituents taken into account. For instance, thermal grain infrared emission can fluoresce atoms or molecules within the cloud, and grain opacity impedes the propagation of emission lines out of the cloud.

As part of the maintenance and improvement of the underlying atomic database we have incorporated the UTA data computed by @Kisielius2003 and updated by Ferland et al. (2013, ApJ submitted). The radiative rates are very large and the corresponding radiative damping parameters can be substantially greater than 1. To improve the accuracy of predictions for such lines with high damping parameters, the Voigt function used in line transfer calculations now uses the fast and accurate implementation of @Wells1999, with some modifications to improve performance in the limit of small damping using the theory presented by @Hjerting1938.

The code includes a number of hyperfine structure lines [@Goddard.W03Hyperfine-Structure-Emission-and-Absorption]. We have put a major emphasis into the physics of the  21 cm line [@Shaw.G05Role-of-molecular-hydrogen-in-star-forming]. The line is usually optically thick, and pumping by   is treated as in @Deguchi1985. @Pellegrini2007 show how feedback between the stellar SED and fluorescent excitation of   alters the 21 cm optical depth and spin temperature in the M17  region.

Solution of material state
--------------------------

Cloudy solves for self-consistent , and the charge state and temperature distribution of a spectrum of grains. This population balance is then used by outer solvers for the material temperature, pressure (when required) and radiation transfer, and used to predict the resulting spectrum. As described in the previous sections, the range of physics treated has expanded significantly since the last review.

At the simplest level, the method of solution for the populations remains similar to that used previously, with the populations of different systems solved for in an iterative loop. However, the increasing level of coupling between species has required significant work to be done on improving the robustness of convergence and the self-consistency of the solutions. The broad range of physics which Cloudy solves in a self-consistent fashion means that the resulting system has a very wide range of physical timescales, a significant challenge for any numerical scheme; the status of the code as an open resource for the astrophysical community means that this challenge must be met by means which require no user intervention. Details of this work will be described in detail in Williams et al. (2013, “Hierarchical physics”, in preparation), so in the present paper we will summarize the overall approach.

The molecular chemistry is now fully consistent with the ionization balance. This is done by scaling chemical reactions including atoms and positive ions to have rates dependent on a combined species which is the sum of the ions within the chemical network solver. The resultant effects of the chemistry on the ionization balance are allowed for by having the molecular network calculate the net source and sink rates for these ions from molecular chemistry, which are then included in the ionization balance solver.

The nonlinear system for the molecular chemistry is solved using an adaptive timestep implicit solver: in typical usage, this is set up to run the chemical balance to full late time equilibrium. The molecular chemistry and ionization solvers have been adapted to allow the solution of time-dependent and steady-state advective flows, using the approach described by @HenneyEtAl05. in the existing equilibrium solvers to adapt them to find the solution of a backward Euler implicit system for the time advance of the state. This allows the code to take advantage of the existing equilibrium solvers with minimal change, exploiting these to provide an implicit solution for the time advance, which is necessary given the wide range of physical timescales which operate in the systems which are modeled. by @HenneyEtAl05 to treat the structure of  region photo-evaporation fronts, and extended by @Henney2007 to the molecular knots within planetary nebulae, including the PDR ahead of the  region.

Work has also gone into improving the robustness of handling of processes which couple different parts of the system sufficiently strongly that a simple iterative scheme was slowly convergent. The particular cases where this has been found to be an issue is in the handling of the resonant O/H charge transfer process [@Stancil.P99Charge-transfer-in-collisions-of-O-with], the rate of which can by far dominate direct ionization processes for either ion, and the handling of Rydberg levels which are more strongly coupled to the continuum than the base ion. A simple solution acceleration approach has been found to be sufficient to allow rapid convergence, details of which will be given in Williams et al (in preparation).

In addition, the solvers for the electron density and temperature have been completely rewritten and are now much more robust and reliable than the previous versions. These changes are especially important for modeling extreme environments: very cold regions on the one hand (PDRs and molecular regions) and very hot on the other (extreme X-ray environments such as disks surrounding a black hole).

Momentum balance and the equation of state {#sec:eos-momentum}
------------------------------------------

allows great flexibility in specifying the spatial distribution of the gas (and grain density). In addition to constant density and various ad hoc density laws (e.g., power law), it is also possible to allow to solve for a self-consistent density distribution based on momentum balance. In the simplest case of a static configuration with no external forces this reduces to the requirement of constant total pressure, whereas with the addition of an external force, such as continuum radiation pressure [@BFM] or gravity [@AscasibarDiaz2009], it becomes hydrostatic equilibrium. also allows the further generalization to the case of dynamic equilibrium in the presence of steady-state gas flow [@HenneyEtAl05].

In addition to the thermal pressure of the gas, also considers other contributions to the total pressure. These include trapped resonance line radiation [@Elitzur1986], ordered and disordered magnetic fields [@HenneyEtAl05; @Pellegrini2007], cosmic rays [@Shaw.G09Rotationally-Warm-Molecular-Hydrogen], and the Reynolds stress due to turbulent motions.

Some computational details
==========================

An Open Source Project
----------------------

 is openly available on the web at the site [[](http://www.nublado.org)]{}. This includes the full source, the atomic and molecular data needed for  to operate, its documentation , and an extensive suite of test cases. The test suite illustrates how to use the code, its range of validity, and includes embedded monitors that confirm that the code is operating correctly. All previously released versions of  are available on the web site, with most stored in an openly accessible *Subversion* repository. The distribution is subject to a BSD style license.

Although this is the first major review since F98,  has been continuously developed, as witnessed by the papers cited above. New versions are released every two to three years, at the conclusion of a period of development which focused on particular aspects of the simulations. This is the seventh major release since F98. The [[](http://www.nublado.org)]{} site gives the full history.

At the time of F98, the code was $\sim 9\times 10^4$ lines written in Fortran 77, with some portable extensions. It transitioned to C in 1999 and C++ in 2006. Today  consists of roughly $2.1\times 10^5$ lines of C++ source code. In recent years, the size of the code has stabilized, as work to extend its scope is balanced by the move from the inclusion of physics data within the source coding towards the use of external database files.

The code has a broad user community. As described below, we ask that users reference paper if the code is used in a publication. At the time of this writing there are nearly 200 citations to F98 or the code’s documentation each year. We maintain a [discussion board](http://tech.groups.yahoo.com/group/cloudy_simulations)[^5] where users can ask questions and where we post announcements of updates to the code.

### Material available on [[](http://www.nublado.org)]{}

[[](http://www.nublado.org)]{}, the code’s web site, gives complete access to files and information about . The site, built using [trac](http://trac.edgewall.org/), gives top level links to a variety of items. These include:

-   [Step by step instructions](http://www.nublado.org/wiki/StepByStep) for downloading, building, and running the code. We provide links to a large tarball for the download, and makefiles are used to build the code.

-   [Stellar atmospheres](http://www.nublado.org/wiki/StellarAtmospheres). As discussed in Section \[sec:StellarAtmospheres\] below, it is possible to use grids of stellar atmospheres in deriving the incident radiation field. This page gives more details and links to available grids.

-   [Known problems](http://www.nublado.org/wiki/KnownProblems), and [hot fixes](http://www.nublado.org/wiki/HotFixes). No code is perfect. Users should post questions or bug reports on the [discussion board](http://tech.groups.yahoo.com/group/cloudy_simulations). We provide a list of known problems. These are deficiencies which we know about, but which have not been fixed in the current version. Hot fixes are small changes to the source code which will fix problems discovered since the last version of the code was released. They should be applied to the code source before building it.

-   The [revision history](http://www.nublado.org/wiki/RevisionHistory) gives a list of all changes to each version of the code. The current review paper gives an overview of changes but is not meant to be complete.

-   A [FAQ](http://www.nublado.org/wiki/FaqPage) page

-   A summary of [old versions](https://www.nublado.org/wiki/CloudyOld) of , including links to download them.

-   The [developer pages](https://www.nublado.org/wiki/DeveloperPages) give links to our notes about developing . You are most welcome to help!

### The [Subversion](http://viewvc.nublado.org/index.cgi/?root=cloudy) repository

The code, its documentation, test suite, and data live in a [Subversion](http://subversion.apache.org/) repository. The layout in the repository is conventional. The *trunk* is the development version and is changed on a near-daily basis.

*Branches* usually originate as copies of the trunk and can be separated into development branches (to add new functionality) and release branches. There is a C13 release branch which split off from the trunk in late 2012. This branch is updated as bugs are fixed but no new code development is done here.

*Tags* are copies of a branch or trunk version. These do not change. Released versions are tags. For instance, the first release of C13 has the tag [*C13.00*]{}.

To assure the quality of the code, we run the test suite of the trunk on a nightly basis provided there are changes. We also test the active release branches and certain key development branches on a similar basis (albeit somewhat less frequently). Additionally we test for common programming errors such as array bounds violations and the use of uninitialized variables once or twice a week. This way many errors can be caught quickly, preventing them from causing problems in a release.

As an open source project, the entire repository is open to public view and download. All versions of the code after the creation of the repository in late 2005 are available. Older versions are maintained as separate tarballs on the [old versions](https://www.nublado.org/wiki/CloudyOld) page.

Testing
-------

The  team has long participated in open comparisons of model predictions. Such comparisons are a valuable way to exchange ideas and find problems, and are the only way to validate projects as complicated as a modern spectral synthesis code [@FerlandReliability01].

Two meetings had been held by the time of F98 to compare predictions for ionized regions, and a third was held soon after. The first was organized by Daniel P[é]{}quignot in Paris in 1985 [@Pequignot.D86Comparison-of-Photoionization-and-steady-shock] but has no on-line proceedings available. Two meetings were held in Lexington, the first a satellite of the STScI meeting in honor of the 70$^{\rm th}$ birthdays of Don Osterbrock and Mike Seaton [@1995aelm.conf.....W], and a second as part of the Conference *Spectroscopic Challenges of Photoionized Plasmas* [@2001PASP..113.1024F]. These comparisons are presented in @Ferland1995 and @Pequignot2001 respectively. Agreement at the 20% - 30% level for most important quantities was achieved by codes that participated in the workshops.

A second form of testing is accomplished by running the code into well-posed physical limits. Correct behavior in limiting cases gives some assurance that intermediate regimes are valid. Examples include the Compton, LTE, molecular and low-density limits discussed in Section \[sec:GridExtreme\] below. The code distribution includes a large test suite which exercises the code over its range of validity, and includes embedded monitors that check that it obtains the expected result. The test suite is designed so that the code can be automatically validated with little effort.

The scope of the simulations has been expanded to include atomic and molecular regions in addition to ionized gas. The goal is to calculate physical conditions of adjacent , , and  regions (or regions, PDRs, and molecular clouds) self-consistently. We begin a calculation at the face of a cloud illuminated by a hot O star and end in cold regions completely shielded from UV radiation [see @Abel2005]. Such a calculation is a better representation of what actually goes on in nature, where , , and  regions are physically adjacent and the properties of each region depend on the radiative and dynamical coupling between the regions. This type of calculation is particularly advantageous in environments where the observed emission could come from more than one region.

 is well tested and in agreement with other spectral synthesis codes that specialize in PDR modeling. The 2004 Leiden PDR meeting compared the results of several PDR codes for 8 benchmark calculations. These calculations are summarized in @Roellig2007. These results show that  agrees very well with the / and // transition, the dependence of other molecules with depth, temperature structure, and FIR emission-line spectrum. The  test suite includes PDR calculations with parameters used by @Tielens1985a, @Kaufman.M99Far-Infrared-and-Submillimeter-Emission-from, and @Le-Petit.F04H3-and-other-species-in-the-diffuse-cloud. These also agree fairly well. For more information, see the test suite input scripts that come with   along with @Abel2005 and @Shaw2005.

@Abel2008 describe some differences between our predictions and those of the PDR codes discussed in @Roellig2007. These largely are the result of our use of elementary processes rather than fitting formulae in determining the physical state. We show examples of this below.

Assessing the effects of uncertainties in atomic / molecular physics rates
--------------------------------------------------------------------------

The atomic and molecular data set needed for a full simulation of the microphysics of a non-equilibrium gas is vast, including ionization, dissociation, and recombination data for all species, along with internal energies, transition probabilities, and collision rates. Many data are the results of theoretical calculations which are at the forefront of research in computational atomic/molecular physics. There will be gaps in these data, and in many cases, basic uncertainties. @Aggarwal.K13Assessment-of-atomic-data:-problems review sources of these uncertainties while @Bautista.M13Uncertainties-in-Atomic-Data-and-Their look into how they propagate through spectral simulations, both with emphasis on ions. @Wakelam.V10Sensitivity-analyses-of-dense-cloud do a similar study of molecular environments.

We have long included the ability to add a component of Monte Carlo Gaussian noise, with a specified amplitude and FWHM, in our simulations. We specify which component of the data to afflict with the noise, its amplitude, and dispersion. The data are altered when the code initializes and the disturbed data are used throughout the calculation. Each rerun of the simulation will have a different set of noise, as determined by randomly sampling the Gaussian distribution. In many cases the noise is large - uncertainties can be as large as 0.2 - 0.5 dex. The random numbers are Gaussian distributed in log space for this reason.

This capability is designed into  to make it easy to examine the effects of uncertainties. It has been applied in several studies. @Shaw2005 investigated the effects of uncertainties introduced by missing collisional rates for . @Porter.R09Uncertainties-in-theoretical-HeI-emissivities:-HII-regions and @Porter.R12Improved-He-I-emissivities-in-the-case-B-approximation documented how uncertainties in the photoionization cross sections, transition probabilities, and collisional rates affect predictions of Case B recombination coefficients.

Such studies make it possible to quantify how known uncertainties propagate into the computed physical conditions or spectrum. It is important to remember that, in many cases, the dominant uncertainties are due to physical processes which are not yet included. Early simulations of  regions and planetary nebulae failed because the importance of charge exchange and dielectronic recombination was not understood. As with any systematic error, the magnitude of the uncertainties can often only be known once they are removed.

Modeling observations {#optimize}
---------------------

Observers are often faced with the problem that they have a set of observations of a particular object and want to derive physical properties of the object from these. The observations typically comprise spectral line and continuum fluxes at various wavelengths and possibly other observables. From these they would like to derive physical properties of the irradiating source (e.g., the effective temperature) and/or the surrounding gas (e.g., the density, electron temperature, chemical abundances, etc).

One way of achieving this is to “reverse engineer” the  model by assuming values for the physical parameters, calculating the model and comparing the results to the observations. The quality of the fit is measured by a $\chi^2$ value. The problem then reduces to finding the set of input parameters that produce the best fit which is the lowest $\chi^2$ value. This is a standard mathematical problem.

 has an *optimize* command that makes carrying out this task easy. The heart of this command is the minimization algorithm for the $\chi^2$ function. There are two algorithms built in for doing this. The oldest one is the SUBPLEX algorithm [@Rowan1990]. This is a generalization of the well-known downhill simplex method AMOEBA. The second algorithm (which is the default) is PHYMIR [@VanHoof1997] which was specifically designed for use in . Both algorithms are robust against noisy functions which is a very important feature since  predictions are always noisy due to the use of adaptive stepsize algorithms and finite precision iterative schemes.

The PHYMIR algorithm has two additional advantages that make it the preferred method over the SUBPLEX algorithm. The first is that the PHYMIR optimization process can be parallelized. If $N$ input parameters are varied then up to $2N$ cores can be used simultaneously which can greatly speed up the calculation. This is discussed in more detail in Sect. \[parallel\].

The second advantage is that the PHYMIR algorithm periodically writes out state files which can be used to restart an optimization run that failed (e.g., due to a power failure) or ran into the maximum number of iterations before the minimum was reached.

Creating grids of calculations {#sec:grid}
------------------------------

The  *grid* command, initially described by @Porter2006, makes it possible to vary input parameters to create large grids of calculations. Several parameters can be varied and the result of the calculation will be predictions for each of the grid points. Figure \[fig:SaveLineRatio\] shows an example where a range of gas kinetic temperature and density were computed and the \[\] $\lambda\lambda 5007, 4363$Å lines were saved. Such diagrams can be used to deduce physical conditions in a cloud.

![ \[fig:SaveLineRatio\] The \[O III\] line $\lambda \lambda 5007 / 4363$ intensity ratio as a function of density and temperature. ](SaveLineRatio){width="\linewidth"}

Predictions are usually saved with one of the *save* commands described in 1. The predictions will normally be brought together into large files which contain the output from all grid points.

Cloudy on parallel computers {#parallel}
----------------------------

Two  commands can take advantage of multi-core computers and high performance computing (HPC) clusters. These are the *optimize* and [ *grid*]{} commands described in Sects. \[optimize\] and \[sec:grid\]. They run  as an “embarrassingly parallel” application, putting one model on each CPU core. Run this way, the code can achieve nearly 100% efficient use of parallel computers.

The parallelization is implemented using two different techniques. The oldest technique uses the *fork* system call that is available under all UNIX operating systems and Apple Darwin. The great advantages of this technique are that no external libraries are needed (i.e. it works “out of the box”) and superior fault tolerance. All the work is done by the child processes, so even in the unlikely event of a crash the parent process can continue, preventing all work from being lost. The big disadvantage of this technique is that it will only work on shared-memory machines so that it cannot take advantage of modern HPC clusters. Currently only the *optimize* command uses this technique.

In view of the fact that HPC computing is moving away from shared-memory machines towards large HPC clusters, we decided to redesign the parallel infrastructure of the code. Since version C10.00 we support parallelization under MPI version 2 or newer. Both the *optimize* and *grid* commands can use this technique. The big advantage is that we can benefit from large HPC clusters, which is important for large grid calculations which can now be run as massively parallel applications. The disadvantage is that the user may need to install an MPI environment or at the very least needs to get familiar with MPI which unfortunately is not always as intuitive as we would like.

Building complex geometries {#sec:Cloudy3D}
---------------------------

Clumping can be included and complex source geometries can be simulated. There are several general considerations.

There are powerful selection effects governing the formation of emission lines when a range of densities exist. You will tend to observe the highest-density regions because the emission per unit volume is proportional to the square density if the line is below its critical density (AGN3 Section 3.5). Only with a fine-tuned mix of densities, where the volume of material at each density exactly compensates for the change in emissivity, will an observer notice emission from a range of densities. Claims that a range of densities contribute to a single emission line should be met with some skepticism. This would require an amazing coincidence [@Ferland.G11Molecular-hydrogen-in-NLSys-and-its-implications].

But clumps do exist. If the clump size is small compared with the physical thickness of the region then they can be treated with a filling factor (see @Osterbrock1959 and AGN3 Section 5.9). In this case the gas is modeled as small clumps that are surrounded by vacuum or much lower-density gas. This is done by simply including the *filling factor* command to specify the fraction of the volume that is filled by clumps.

If the clumps are larger than the physical thickness of the region then each clump will have its own ionization structure. This is the “LOC” model of quasar emission-line clouds described by @Baldwin1995 and @FergusonKoristaBaldwin1997. The model is developed in several papers by the same team. In this case we compute grids of models and save the results. The spectra are then co-added using distribution functions to describe the range of cloud properties. The final spectrum depends on these distribution functions. @GiammancoEtAl_clumpsCloudy04 show  calculations where optically thick clumps are present in the ISM.

In practice we normally use the *grid* command (S\[sec:grid\]) for this, but there are circumstances where complex changes in parameters may be needed. The program in the programs directory in the code distribution computes a grid of models and extracts the predictions using MPI on a distributed memory machine.

Another approach is for a driving program to use  to compute differential volume elements of a large and complex structure, and then integrate to get the next emission. An example is the <span style="font-variant:small-caps;">Cloudy\_3D</span> code[^6] described in @MorissetCloudy3D06 and @MorissetStasinskaCloudy3D08. <span style="font-variant:small-caps;">Cloudy\_3D</span> was used to compute the image shown in Figure \[fig:Cloudy3D\]. The more recent <span style="font-variant:small-caps;">pycloudy</span> code by the same author[^7] is a more general tool for controlling and analyzing multiple runs via scripts. The <span style="font-variant:small-caps;">Rainy3D</span> code is another example [@MoraesDiazRAINY09].

![image](posterhazy02){width="0.7\linewidth"}

Spectral energy distributions from stellar atmosphere grids {#sec:StellarAtmospheres}
-----------------------------------------------------------

The heart of any photoionization simulation is the SED of the incident radiation field. It is this energy which is reprocessed by the cloud to produce the observed nebular emission. Several groups have created large grids of stellar SEDs using advanced stellar atmosphere codes. Other groups have used these data to create stellar population synthesis models that give the integrated spectrum of a galaxy as a function of time after a starburst.  can interpolate on SED grids having an arbitrary number of dimensions (these might include surface temperature, gravity, chemical composition, mass loss rate, age, etc) and include this in the incident radiation field.

Figure \[fig:StarsContinuum\] compares predictions for five of the 54 K SEDs that are available. These include a blackbody and atmospheres computed by @Mihalas1972, @Kurucz1979, @Kurucz1991 and @Rauch2003. All were normalized to have the same total luminosity ($10^{38}$ erg s$^{-1}$) observed from a distance of $10^{18}$ cm. Note the order of magnitude dispersion among the continua for energies around 4 Ryd. This can have a major effect on the  modeling results, showing the crucial role that the stellar SED plays.

![\[fig:StarsContinuum\]This figure shows the emergent radiation field predicted by five 54 K stars included with the code. The smoothest is the blackbody, and the @Kurucz1991 and @Rauch1997 atmospheres show the most structure.](StarsContinuum){width="\linewidth"}

Numerous stellar grids can be used with little additional work. In some cases the SED data are stored on the author’s web site while in others they are stored on the  web site. A convenient page providing links to all the necessary files can be found at [wiki.nublado.org/wiki/StellarAtmospheres](http://wiki.nublado.org/wiki/StellarAtmospheres). This web page also gives more computational details and links to the papers describing the grids.

These are the most important SED grids currently supported by :

1.  The Atlas grids. There are two versions of these, the preferred one being the new-ODF grids described in @Castelli2004. The older generation of Atlas model atmospheres described in @Kurucz1991 is also supported. These grids can be useful if you need models for extreme metalicities not covered by the new-ODF grids. Both grids contain LTE, plane-parallel, hydrostatic model atmospheres with effective temperatures ranging between 3500 and 50000 K.

2.  The Tlusty OSTAR2002 and BSTAR2006 grids described in @Lanz2003 and @Lanz2007. These grids contain non-LTE, line-blanketed, plane-parallel, hydrostatic O and B star SEDs. We also support merged OSTAR2002/BSTAR2006 grids covering a temperature range between 15000 and 55000 K.

3.  The WMbasic and CoStar O and B star grids. These are two small grids of non-LTE, line-blanketed, and wind-blanketed models. The first grid is described in @Pauldrach2001 and the second in @Schaerer1997.

4.  All PN central star SED grids computed by T. Rauch. These include the H-Ni, PG 1159, and C/O white dwarf grids, as well as the pure hydrogen, pure helium and H+He grids. The older H-Ca grids are also supported, though for most purposes they have been superseded by the H-Ni grids (unless you need models with T$_{\rm eff} >$ 190000 K). All grids contain non-LTE, line-blanketed, plane-parallel, hydrostatic model atmospheres. They are described in @Rauch1997 and @Rauch2003. The temperature range typically is between 50000 and 190000 K, though some grids have a different range.

5.  Stellar population synthesis models from Starburst99 [@Leitherer1999] and PopStar [@Molla2009]. Typically the user will do their own run and generate a  grid from the output using either  commands (Starburst99) or a script supplied on the  web site (PopStar).

Many of the grids are very large and accessing them as ASCII files would be slow. They are “compiled” to create direct access binary files as part of the installation procedure. Once complete the stellar SEDs can then be accessed with the appropriate command in the simulation control deck. Logarithmic interpolation is done to create model atmospheres with any set of specified parameters using nearby models from the original grid.

The code is very flexible and allows users to create their own SED grids, e.g., from a Starburst99 run. As a result we can also add support for new grids during a release cycle when the need arises. This is possible because no code changes are needed to do this.

Citing  and its underlying databases
------------------------------------

 is a research project that involves the creative efforts of many people. When used in publications it should be cited as follows: “*Calculations were performed with version C13.00 of , last described by Ferland et al. (2013)*”, where this paper is the reference. The specific version of the code, written as C13.00 in this example, should be given so that, in case any future questions arise, it will be possible to reproduce the calculation using the archived versions on [www.nublado.org](http://www.nublado.org).

We are now moving the atomic and molecular data to external databases. These are replacing our internal database, which had been embedded in the source. Many recombination coefficients are based on @BadnellEtAl03 and @Badnell06 and posted on the web sites <http://amdpp.phys.strath.ac.uk/tamoc/RR/> and <http://amdpp.phys.strath.ac.uk/tamoc/DR/>. Much of the molecular emission data is from LAMDA [@Schoier.F05An-atomic-and-molecular-database-for-analysis] as accessed on Dec. 18, 2010, as well as the JPL [@Pickett1998] and CDMS [@Mueller2001; @Mueller2005] databases. Much of the ionic emission data is from CHIANTI, as described by @Dere.K97CHIANTI---an-atomic-database-for-emission and @Landi2012, using version 7.0. Much of the  data is from @Wrathmall2007, @Abgrall.H94The-Bprime1Sigmau-rarr-X1Sigmag-and-D1Piu, and the Meudon web site (<http://molat.obspm.fr/index.php?page=pages/Molecules/H2/H2can94.php>).

All of these databases play a major role in most calculations. We ask that users cite both  itself, and those underlying databases, in any publications. These databases can only thrive if their role is properly acknowledged. We provide a *print citations* command that will provide the correct citations in a format that can be easily copied and pasted into papers.

Applications
============

XDRs
----

The term X-ray Dissociation Region or “XDR” was coined by @Maloney1996 to describe atomic regions near X-ray sources. (A somewhat similar calculation had been presented by @Ferland1994 in the context of optical filaments in cool-core clusters of galaxies.) In keeping with traditions established in the study of PDRs, a truncated SED, including only photons between 1 – 100 , was considered. Figure \[fig:XdrAgnSed\] shows the @Maloney1996 X-ray continuum together with the mean AGN SED derived by @Mathews1987. Both SEDs are built into .

![ \[fig:XdrAgnSed\] The thin line shows a mean AGN SED as deduced by @Mathews1987. The thick line shows the truncated SED considered in XDR calculations. Both SEDs are built into  and were normalized to have the same X-ray flux.](XdrAgnSed){width="\linewidth"}

A second Leiden meeting on radiatively excited atomic and molecular regions was held in 2012[^8]. The meeting considered five PDR and four XDR simulations. The web site[^9] gives some details along with results of the participating codes. We agreed with the PDR results, to within the considerable scatter, as had been found by @Roellig2007 and @Abel2008.

This simulation has a hydrogen density of $n_{\rm H} = 10^3 \pcc$, the hydrogen column density given above (corresponding to point-source  $\sim 10^3$ mag) and an X-ray flux of 270 . The gas ionization is proportional to the dimensionless ionization parameter $$U = \frac{\phi({\rm H})}{c \ n({\rm H})}$$ where $\phi({\rm H})$ is the flux of hydrogen ionizing photons [@AGN3]. This simulation has the highest ionization parameter of the XDR tests,

Photoionization by the incident radiation field, and by diffuse EUV line emission, emission lines produced by the XDR gas, produces a moderate level of ionization throughout the XDR2 cloud. Figure \[fig:XDR2Ionization\] shows the ionization fractions for H, He, and C as a function of the point-source . There are significant amounts of doubly ionized species. The most important of the ions shown is , which destroys CO by charge exchange.

![ \[fig:XDR2Ionization\] The ionization fractions for H, He, and C are shown as a function of depth, expressed as the point-source .](XDR2Ionization){width="\linewidth"}

The lower panel shows the absorption and scattering opacities, where the latter includes the factor $(1-g)$ which discounts forward scattering. Grains are the major contributor to the total scattering across the UV and optical, with Thomson scattering being dominant at the shortest and longest wavelengths. Rayleigh scattering is responsible for the feature at $\sim 0.1216 \micron$.

Grains are the dominant absorption opacity source across the UV, optical, and IR, with several resonant features visible. Below 0.0912  gas opacity dominates, with the strongest edge due to . Edges due to He and inner shells of the heavy elements are visible at shorter wavelengths.

The upper panel shows the local photon interaction rate, $\phi_{\nu} \alpha_{\nu}$, where $\alpha_{\nu}$ is the gas opacity \[\]. The local photon flux $\phi_{\nu} = 4 \pi \bar{J}/h\nu$ \[\] includes all components of the including the attenuated incident SED and the local diffuse line and continuous emission. The strong lines in the FUV and EUV[^10] are the result of the solution of the many-level iso-sequence atoms as described in previous sections, and have intensities fully self consistent with the opacities shown in the lower panel, the level of ionization, gas temperature, and optical depths. There are significant sources of ionizing radiation in addition to the attenuated incident XDR continuum. The most important are EUV recombination lines of and . Direct photoionization by the incident continuum produces a trace amount of while the EUV emission lines, together with the incident continuum, produce a moderate amount of  and .

![ \[fig:XDR2midplaneRadField\] Components of the radiation field at a depth corresponding to  = 5 mag in the XDR2 simulation. The lower panel shows the absorption and scattering opacities while the upper panel shows the photon interaction rate.](XDR2midplaneRadField){width="\linewidth"}

Figure \[fig:XDR2midplaneZoom\] shows a zoom into the photon interaction rate $\phi_{\nu} \alpha_{\nu}$ within the hydrogen-ionizing radiation field. The total photoionization rate of a species is the integral of $\phi_{\nu}$ over the photoionization cross section $\sigma_{\nu}$ (AGN3), while the Figure shows the total interaction rate, evaluated using the computed total opacities. The horizontal lines indicate the range of wavelengths which can ionize H and He. For reference, photoionization cross sections fall off with decreasing wavelength as a power law ranging from $\sigma_{\nu} \sim \lambda^{-3}$ for  and  to $\sigma_{\nu} \sim \lambda^{-1}$ for . The attenuated incident XDR continuum is the dominant contributor to the   photoionization rate. Recombination lines of , with a significant contribution from the XDR continuum, dominate photoionization of , and these lines, together with lines, ionize .

![ \[fig:XDR2midplaneZoom\] This shows a small portion of the radiation field around wavelengths capable of ionizing H and He. The horizontal lines indicate thresholds for photoionization of , , and .](XDR2midplaneZoom){width="\linewidth"}

The cumulative effect is that a significant amount of  exists across of the cloud (Figure \[fig:XDR2Ionization\]). Charge transfer between  and CO is the dominant CO destruction process in regions that are well-shielded from FUV radiation: $${\rm He}^+  + {\rm{CO}} \to {\rm He} + {\rm C}^+ + {\rm O}$$ [@Anicich.V77Product-distributions-for-some-thermal; @Laudenslager.J74Near-thermal-energy-charge]. The large amount of  results in efficient destruction of CO, and little CO exists as a result. The strong , , and recombination lines heat the gas through both direct photoionization, and grain electron ejection.

![ \[fig:Xdr2Con\] ](XDR2_con){width="\linewidth"}

![ \[fig:XdrCoIntensity\] ](XDR-CO-intensity){width="\linewidth"}

Table \[tab:XDR2lines\] gives the intensities of lines predicted to have emergent intensities brighter than 10% of the \[\] $\lambda 158 \micron$ line in the format used by . We intentionally present Table \[tab:XDR2lines\] in the format used by the code, as an introduction to its output. Each line is indicated by a label in the  output, given as the first column in the Table, and a wavelength, given in the second column. The line label uses the compact notation “`C  2`” for \[\] so the label is intended for identification rather than spectroscopic notation. In the  output, as in Table \[tab:XDR2lines\], “`A`” and “`m`” indicate Å and respectively. A very large number of lines are predicted by . We provide a *save line labels* command which will create a list of all emission lines with labels, wavelengths, and a comment indicating the line’s origin. There is a discussion of the various line entries in Part 2 of , the code’s documentation.

The third column gives the intensity, $4 \pi J({\rm line})$ \[\] of each line. This is the total emission radiated into $4\pi \sr$ from a unit area of cloud,

[@ cr@ cr @]{} Spectrum & & Intensity &\
`FeKa` & `1.78A` & $5.07\e{-1}$ & 5.16\
`H  1` & `1216A` & $1.73\e{-1}$ & 1.76\
`C  2` & `2326A` & $1.29\e{-1}$ & 1.32\
`O  2` & `3727A` & $1.34\e{-2}$ & 0.137\
`O II` & `3729A` & $1.63\e{-1}$ & 1.66\
`S  2` & `4074A` & $7.55\e{-2}$ & 0.768\
`N  1` & `5199A` & $1.66\e{-2}$ & 0.169\
`O  1` & `6300A` & $4.57\e{-2}$ & 0.464\
`O  1` & `6363A` & $1.48\e{-2}$ & 0.151\
`H  1` & `6563A` & $2.06\e{-2}$ & 0.209\
`H-CT` & `6563A` & $5.55\e{-3}$ & 0.056\
`N  2` & `6584A` & $2.16\e{-2}$ & 0.219\
`S II` & `6716A` & $1.77\e{-2}$ & 0.181\
`S II` & `6731A` & $1.43\e{-2}$ & 0.145\
`S  3` & `9069A` & $1.24\e{-2}$ & 0.126\
`S  3` & `9532A` & $3.22\e{-2}$ & 0.327\
`C  1` & `9850A` & $3.32\e{-2}$ & 0.338\
`S II` & `1.029m` & $1.34\e{-2}$ & 0.136\
`S II` & `1.032m` & $1.83\e{-2}$ & 0.186\
`S II` & `1.034m` & $1.30\e{-2}$ & 0.132\
`He 1` & `1.083m` & $2.34\e{-2}$ & 0.238\
`S  3` & `18.67m` & $9.64\e{-2}$ & 0.980\
`S  3` & `33.47m` & $3.37\e{-1}$ & 3.42\
`Si 2` & `34.81m` & $1.50\e{-1}$ & 1.52\
`O  1` & `63.17m` & $6.04\e{-1}$ & 6.14\
`O  1` & `145.5m` & $1.48\e{-1}$ & 1.50\
`C  2` & `157.6m` & $9.84\e{-2}$ & 1.00\
`C  1` & `369.7m` & $1.40\e{-2}$ & 0.142\
`TIR ` & `1800m` & $1.54\e{-1}$ & 1.56\
`TALL` & `10000A` & $4.83\e{+1}$ & 491.\
\

Several lines deserve special mention. “FeKa 1.78A” is the Fe K$\alpha$ X-ray line and is mainly produced by Fe locked in silicate grains. The  line is predicted to have a significant contribution by mutual neutralization excitation, indicated by the label “H-CT”. This is the process $${\rm H}^-  + {\rm{p}} \to {\rm H}^*(n=3) + {\rm H}(1s)$$ [@Peart1985; @Ferland1989]. We assume that $n=3$ is statistically populated.

XDRs and AGN
------------

Figure \[fig:XdrAgnSed\] shows that the XDR continuum is but a small part of the total SED of an AGN. The lighter line is the mean AGN SED derived by @Mathews1987 and built into . The XDR SED is meant as a way to compute the conditions in the  region that lies behind (as seen from the central object) the  region where most hydrogen-ionizing radiation is absorbed. As @Abel2005 stress, this may be a great oversimplification.

To check this, we computed two models of a “Narrow Lined Region” (NLR) cloud near an AGN. These are lower-density dusty regions that contribute to the optical spectrum and are likely to be ionized layers on the surface of larger molecular clouds (AGN3). Many different sets of chemical abundances are built into . We use a standard ISM gas-phase composition and a mixture of graphite and silicate grains combined with an MRN size distribution. The clouds have a physical thickness corresponding to a point-source  of $\sim 10^3$ mag. Galactic background cosmic rays were assumed. We now adopt the @Indriolo2007 mean  cosmic ray ionization rate of $2\times 10^{-16} \ps$ as the default Galactic background.

The AGN radiation field intensity was set with the dimensionless ionization parameter $U$, defined as the ratio of hydrogen-ionizing photon to hydrogen densities (AGN3). We adopt $\log U = -1.5$, a typical value deduced from the optical emission line spectrum [@FergusonKoristaBaldwin1997]. We normalized the XDR continuum to have the same 1 - 100  flux, 83.18 , as the AGN continuum. The hope is that an XDR computed with this SED and flux would be similar to the  region in the AGN cloud.

The equation of state relates the gas density to other physical quantities such as the kinetic temperature or radiation pressure (see Section \[sec:eos-momentum\] above). We assume constant total pressure. This is very important for the AGN continuum which produces a hot ($\sim 1\e4 \K$) layer of ionized gas on the surface of the cloud. The illuminated face of the XDR cloud is predominantly atomic and warm ($\sim 1\e3 \K$)). As a result, the gas pressures will differ by about 2 dex, the difference in temperature and particle density. We want the densities in the  region to be comparable if we are to make a meaningful comparison between the two simulations. The pressure in the AGN simulation is 5.97 . We use the same pressure in the XDR. The hydrogen density and temperature at the illuminated faces of the AGN and XDR clouds is then $10^4 \pcc$, $1.88\e4 \K$ and $2.33\e6 \pcc$, $2.29\e2 \K$ respectively. Hydrogen in the AGN is fully ionized at this point while the XDR has 29% .

Figure \[fig:XdrAgnTemperature\] shows the gas kinetic temperature and Figure \[fig:XdrAgnDensity\] shows the hydrogen density as a function of depth from the illuminated face of the layer for the two scenarios. Depth is shown in terms of the point-source  to be consistent with other literature. Figure \[fig:XdrAgnHstructure\] shows the distribution of hydrogen in its various forms.

![ \[fig:XdrAgnTemperature\] This compares the gas kinetic temperature for XDR and AGN clouds with the same X-ray flux, the SEDs shown in Figure \[fig:XdrAgnSed\], and the same total pressure.](XDR_AGN_temperature){width="\linewidth"}

![ \[fig:XdrAgnDensity\] The hydrogen density is shown for the XDR and AGN clouds.](XDR_AGN_density){width="\linewidth"}

![ \[fig:XdrAgnHstructure\] The hydrogen atomic, molecular, and ion fractions are shown for the XDR and AGN clouds.](XDR_AGN_H_structure){width="\linewidth"}

There is a very warm ionized layer in the AGN case due to the  zone where hydrogen-ionizing photons are absorbed. This layer produces most of the emission from the cloud since the AGN SED peaks in the FUV and EUV. This emission, predominantly lines in the optical and UV, may be undetectable if there is a large amount of surrounding dusty material. The UV/optical emission would then be reprocessed by other clouds into IR emission. In either case, the XDR continuum misses the majority of the continuum available for reprocessing. We shall predict the emission emergent from the cloud we model, and do not consider further reprocessing by other clouds in the system.

There are surprising large differences in the kinetic temperature in relatively shallow regions of the  region, which begins at about  $\sim 1$. The XDR produces a very flat temperature profile with $T \sim 300 - 400 \K$ being typical. The AGN produces a much warmer   layer at shallow depths with temperatures ranging from $T \sim 3000 \K$ to $T \sim 200 \K$. Soft X-rays that filter through the  layer into the  regions produce the warm gas.

Figure \[fig:XdrAgnAv15\] shows the photon interaction rate at a depth of $ = 1.5$, the region where the temperature differences between the XDR and AGN simulations are the largest. Roughly the same amount of energy is present in the region of the spectrum where the XDR incident continuum is defined, $1 - 100 \keV$, although the transmitted AGN SED extends down to lower energies. This spectral region adds additional heating to the gas. The greatest differences are at $\lambda > 0.1 \micron$, where $\phi_{\nu} \alpha_{\nu}$ is about 3 dex larger in the AGN simulation. This radiation heats the gas through grain electron photoejection, producing the much higher temperature.

![ \[fig:XdrAgnAv15\] The photon interaction rate at $ = 1.5$ is shown for the XDR and AGN simulations. This is the point where the AGN has a much warmer  region, produced by the emission from the adjacent  region.](XDR_AGN_Av15){width="\linewidth"}

Table \[tab:XDRAGNlines\] compares predicted intensities for the XDR and AGN simulations. The XDR approximation does roughly agree with the AGN case for some MIR lines. Standard XDR emission lines such as \[\], \[\], etc, are generally within factors of $0.3 - 0.5 $ dex of one another. Lines from higher ionization species, such as \[\] and \[\] are bright in the AGN but missing from the XDR due to the assumed SED. The AGN produces far more total power since an AGN SED peaks at energies that are not included in a standard XDR calculation.

[ c l l ]{} Line& AGN & XDR\
`H  1  6563A` & $2.76       $ & $4.24 \e{-3}$\
`H  1  4861A` & $8.61 \e{-1}$ & $1.01 \e{-3}$\
`H  1 1.875m` & $3.44 \e{-1}$ & $3.85 \e{-4}$\
`O  2  3727A` & $1.93 \e{-1}$ & $1.34 \e{-7}$\
`O  3  5007A` & $9.99       $ &\
`N  2  6584A` & $3.67 \e{-1}$ &\
`S  2  6720A` & $5.62 \e{-1}$ & $6.66 \e{-21}$\
`H2   2.121m` & $4.72 \e{-3}$ & $4.16 \e{-3}$\
`H2   17.03m` & $2.20 \e{-2}$ & $5.15 \e{-2}$\
`H2   12.28m` & $4.92 \e{-3}$ & $2.40 \e{-2}$\
`H2   9.662m` & $2.68 \e{-3}$ & $2.32 \e{-2}$\
`C  1 609.2m` & $1.03 \e{-4}$ & $1.80 \e{-4}$\
`C  1 369.7m` & $5.61 \e{-4}$ & $9.64 \e{-4}$\
`C  2 157.6m` & $9.48 \e{-3}$ & $3.69 \e{-3}$\
`NE 2 12.81m` & $6.15 \e{-1}$ & $2.51 \e{-1}$\
`NE 3 15.55m` & $1.35       $ & $2.99 \e{-3}$\
`SI 2 34.81m` & $1.34 \e{-1}$ & $1.82 \e{-1}$\
`S  3 18.67m` & $6.60 \e{-1}$ & $1.69 \e{-7}$\
`FE 2 1.644m` & $5.15 \e{-3}$ & $1.53 \e{-13}$\
`C  1 609.2m` & $5.48 \e{-5}$ & $8.54 \e{-5}$\
`C  1 369.7m` & $3.83 \e{-4}$ & $6.26 \e{-4}$\
`C  2 157.6m` & $1.41 \e{-2}$ & $3.37 \e{-3}$\
`O  1 63.17m` & $5.65 \e{-1}$ & $2.51 \e{-1}$\
`O  1 145.5m` & $2.71 \e{-2}$ & $1.83 \e{-2}$\
`Ne 2 12.81m` & $8.83 \e{-1}$ & $3.07 \e{-1}$\
`Ne 3 15.55m` & $2.46       $ & $3.53 \e{-3}$\
`Ne 3 36.01m` & $1.66 \e{-1}$ & $1.17 \e{-4}$\

 excitation diagrams are often used to probe physical conditions in warm molecular regions. Figure \[fig:H2Excitation\] shows this diagram for the two models. The overall distribution of As is typical for such diagrams, lower levels, which can be excited by cooler gas, indicate lower temperatures than the high levels, which are only excited in warmer regions . The populations below 2000  , a temperature something like that of the  region in Figure \[fig:XdrAgnHstructure\].

![ \[fig:H2Excitation\]  excitation diagram. The x-axis is the excitation energy of the v,J levels expressed in . The y-axis is the predicted column density of the v,J level divided by its statistical weight. The lines indicate thermal distributions at various temperatures.](H2_excitation){width="\linewidth"}

![ \[fig:XdrAgnEmergentCont\] ](XDR_AGN_EmergentCont){width="\linewidth"}

Physical conditions over an extreme range of matter and photon density {#sec:GridExtreme}
----------------------------------------------------------------------

The two previous sections highlight the types of physics that has been a particular emphasis in the code’s development since F98, dusty atomic and molecular regions.  is designed to faithfully simulate physical processes that occur in the full range of density and temperature encountered in interstellar clouds, accretion disks, or dense accretion flows. To demonstrate this range we computed (using the *grid* command described above) the properties of a unit cell of a photoionized cloud over a very wide range of density and intensity of the incident radiation field. Such tests are important because they show that predictions agree with analytical theory for asymptotic limiting cases.

The gas has solar composition (without grains) and is illuminated by a blackbody with a $T_{BB} = 10^6$ K color temperature but with a variable intensity. The hydrogen density ranges from $10^{-8}\pcc$, below that of the IGM, to $10^{18} \pcc$, a density that is typical of the atmospheres of some stars or accretion disks. The horizontal axis is the intensity of the black body given as the energy-density temperature, $T_u = (u/a)^{1/4}$ K, where $u$ is the total energy density in all wavelengths \[erg cm$^{-3}$\] and $a$ is the Stefan radiation-density constant. This range in $T_u$ includes environments extending from the Intergalactic Medium (IGM) to deep layers within a star. Most clouds encountered in astrophysics have a gas and energy density that lies somewhere in Figure \[fig:GridExtreme\].

![image](grid_extreme){width="0.8\linewidth"}

The upper left panel of Figure \[fig:GridExtreme\] shows the predicted gas kinetic temperature. This ranges from low values typical of cold molecular gas (the upper left-hand corner of the figure) to high values in the highly ionized right end. The gas temperature closely approaches $T_{BB}$ as $T_u \rightarrow T_{BB}$, as it must from thermodynamics. The right edge of the Figure corresponds to a radiation field in strict thermodynamic equilibrium (STE) since $T_u = T_{color}$. There is no lower bound to the gas kinetic temperature but we do see that generally $T_{kin} > T_u$ since the gas is not a perfect radiator. The lowest temperatures occur for the denser gas at the lowest $T_u$.

The remaining panels of Figure \[fig:GridExtreme\] show the state of hydrogen. Only fractions $n(X)/n({\rm H})> 10^{-5}$ are plotted for simplicity. Bands of constant ionization parameter $U$ are the diagonals running from lower left to upper right. The gas is highly ionized in the lower right half high-$U$ region. Moving to lower $U$, going from the lower right corner towards the upper left, the gas becomes first atomic then molecular. In low-$U$ molecular regions the chemistry occurs totally in the gas phase since grains are not present.

Figure \[fig:GridExtremeCartoon\] is an annotated version of Figure \[fig:GridExtreme\] summarizing some physical limits and showing locations of some astronomical objects. This Figure is meant to illustrate the physics occurring in various combinations of density and radiation field, and is not meant to be rigorous.

Gas in the high-density region of the Figure will be in local thermodynamic equilibrium, LTE, [@Mihalas1978; @Rutten2003] when the density is high enough for thermal collisions to control the ionization and level populations. Levels are said to be in LTE if their level populations are given by Boltzmann statistics for the local gas kinetic temperature. The radiation field may, or may not, be a black body at this temperature. Notice that only some levels of an ion may be in LTE. In general higher levels come into LTE at low densities, because of larger collision cross sections and lower transition probabilities. In C13 only ions of the H- and He-like iso-sequences have enough high levels to go to LTE. Higher densities are needed to go to LTE at larger $T_u$ for two reasons. Both the level of ionization and the illuminating radiation field increase with increasing $T_u$. Higher densities are needed if collisions are to dominate rates for level populations.

![ \[fig:GridExtremeCartoon\] This panel identifies some physical and thermodynamic limits and shows where some regions in Active Galactic Nuclei are located , for the calculations shown in Figure \[fig:GridExtreme\]. A wide range of densities, and various energy-density temperatures of the $10^6 \K$ blackbody, are shown.](grid_extreme_cartoon){width="\linewidth"}

The gas is said to be in strict thermodynamic equilibrium (STE) when the ionization, level populations, and radiation field are given by the same temperature. This occurs at the right edge of Figure \[fig:GridExtremeCartoon\], where $T_u \rightarrow T_{BB}$. Our test suite includes many cases confirming that predictions go to the LTE and STE limits where expected.

The temperature in the lower-right quadrant is determined by Compton energy exchange [@Ferland1988], which drives $T_{kin} \rightarrow T_{BB}$. Here photon - electron collisions dominate the energy exchange and the gas temperature approaches a value determined by the SED of the radiation field. Compton exchange dominates when there is little absorption opacity, which is true for the highest-$U$ regions in the lower right of the diagram.

Classical Strömgren photoionization (AGN3) operates in the mid-$T_u$, low to mid density, regions of the Figure. Here the approximation that most atoms are in the ground state and that all recombinations eventually reach ground (the equivalent two-level atom) is valid.

Grains tend to equilibrate at temperatures near $T_u$ and have sublimation temperatures around $1 - 1.5\e3 \K$. They can only exist in the left third of the diagram.

Figure \[fig:GridExtremeCartoon\] shows where some of the emission-line regions of Active Galactic Nuclei (AGN) are located. The narrow-lined region (NLR) may be molecular clouds irradiated by the radiation field of the AGN. The BLR, likely the skin of an accretion disk near the supermassive black hole, lies between the LTE and Strömgren regimes. This environment is dense enough for there to be significant populations of excited states. Photoionization and collisional ionization from these states, and the radiative transfer effects produced by their large populations, all make this a computationally challenging environment. The molecular torus, the dusty warm molecular gas that exists outside the accretion disk and creates the AGN 1 / AGN 2 distinction (AGN3), lies along the left. The intergalactic medium (IGM) has lower gas density and is illuminated by a weak radiation field.

Figures \[fig:GridExtreme\] and \[fig:GridExtremeCartoon\] show that such diverse phenomena as the IGM, AGN molecular torus, the NLR, and the BLR are simply manifestations of different regimes of atomic and molecular physics. This is the approach we take. If the microphysics is done at an elementary level the macrophysics will follow.

A look forward
==============

The development of  continues. The goal is a true simulation of the microphysics and spectrum of gas and dust over the range of conditions shown in Figure \[fig:GridExtremeCartoon\]. Our calculations have always been limited by processor power, and many important physical processes are at the forefront of research in atomic, molecular, or grain physics. New simulations, which will offer better insight into what happens in front of our telescopes, will be possible with faster computers, improved numerical methods, and better physical data.

The code infrastructure is being improved. We had developed our own database for physical processes, on an ad hoc basis, and embedded it within the C++ source code. This makes it very difficult to update the database as improvements or extensions occur. We are now well into moving physical data into external databases which are parsed when the code is initialized. This effort should be complete by the next release. This external database will be public, along with the rest of .

As described in Section \[parallel\], two commands make it possible to perform a large number of simulations in parallel using MPI. This type of “embarrassingly parallel” calculation is ideal for distributed memory systems.

Shared memory systems should be easier to program and might be used to make single models faster. As described throughout this paper, a calculation simultaneously and self-consistently solves a large number of relatively modest problems. There is no “long pole in the tent” to go after in searching for tasks to make parallel. There are many cross-dependencies between the physical parameters we calculate. Just one example: to calculate the ionization structure you need to know the radiation field, but to calculate the radiation field you need to know the ionization structure. There are many more dependencies like this one, forcing us to use iterative schemes in many places. These make parallellizing the code a lot harder. Worse, the data layout in memory is often less than optimal, resulting in poor cache utilization. Some changes have been made to improve cache locality and the potential for vectorization, but more work remains to be done. The speed of calculations on modern CPU architectures is often limited by memory bandwidth rather than compute speed. We are still considering how to take better advantage of today’s multi-core processors.

To put all this in perspective, our *pn\_paris* test case, one of the simulations from the 1985 Paris meeting, took about a minute to compute on large mainframes at the time. Today the simulation still requires about a minute, despite the astonishing increase in computer power in the past 28 years. Today’s simulation includes many more physical processes, far better emission models, and is a much more robust model of the real nebula.

The grain physics will be improved, driven by the remarkable advances from recent infrared space missions. We will include more grain surface reactions, thought to be important in forming complex molecules. Grain opacities, especially for PAHs, depend on charge and temperature and are not a constant for a particular material and size. Finally, radio emission from spinning grains can be important and is being developed.

Filaments in cool-core clusters of galaxies are thought to be excited by penetrating energetic particles from the surrounding hot intracluster medium [@FerlandFabianEtAl2009; @Fabian.A11The-energy-source-of-the-filaments-around]. Our treatment of cosmic ray or energetic particles does not now include attenuation [@FerlandMushotzky1984], which will depend on an uncertain magnetic field geometry. Theories for cosmic ray transport do exist [@Padovani.M09Cosmic-ray-ionization-of-molecular-clouds] and may be incorporated.

Tests shown in previous papers, and demonstrated in our test suite, show that species treated with our iso-sequence model go to LTE in the high radiation or particle density limits. These include all one and two-electron species. Other ions are treated assuming equivalent two-level systems, as described above, and cannot to go LTE. This is the greatest weakness in our simulations at high densities. We intend to extend the iso-sequence approach to more species, using accurate atomic databases to model lower levels. We can now use both Chianti and Stout, our new external database. Chianti does not include subordinate collisions and so cannot go to LTE. It was intended for relatively low densities. Our Stout database includes all collisions. Neither extends to high enough energy levels for collisional coupling to the continuum and LTE to occur. These models will have to be supplemented with higher Rydberg levels to allow the appropriate high-density behavior to occur.

DEDICATION {#dedication .unnumbered}
==========

*We dedicate this paper to Manuel Peimbert. Through his leadership the traditions established by Menzel, Baker, , Strömgren, Osterbrock & Seaton are being carried on, in Mexico.*

GJF acknowledges support by NSF (1108928; and 1109061), NASA (10-ATP10-0053, 10-ADAP10-0073, and NNX12AH73G), and STScI (HST-AR-12125.01, GO-12560, and HST-GO-12309). PvH acknowledges support from the Belgian Science Policy Office through the ESA Prodex program. WJH acknowledges funding from DGAPA-UNAM, Mexico, through grant PAPIIT IN102012.

[^1]: www.nublado.org

[^2]: <http://amdpp.phys.strath.ac.uk/tamoc/DATA/>

[^3]: CHIANTI is a collaborative project involving the NRL (USA), the Universities of Florence (Italy) and Cambridge (UK), and George Mason University (USA).

[^4]: <http://molat.obspm.fr/index.php?page=pages/Molecules/H2/H2can94.php>

[^5]: <http://tech.groups.yahoo.com/group/cloudy_simulations>

[^6]: <http://sites.google.com/site/cloudy3d/>.

[^7]: <https://sites.google.com/site/pycloudy/home>.

[^8]: <http://www.lorentzcenter.nl/lc/web/2012/482/info.php3?wsid=482>

[^9]: <http://home.strw.leidenuniv.nl/~loenen/LC-CO/>

[^10]:
