3D-PDR: A new three-dimensional astrochemistry code for treating Photodissociation Regions

Photodissociation regions (PDRs) define the transition zone between an ionized and a dark molecular region. They consist of neutral gas which interacts with far-ultraviolet radiation and are characterized by strong infrared line emission. Various numerical codes treating one-dimensional PDRs have been developed in the past, simulating the complexity of chemical reactions occurring and providing a better understanding of the structure of a PDR. In this paper we present the three-dimensional code, 3D-PDR, which can treat PDRs of arbitrary density distribution. The code solves the chemistry and the thermal balance self-consistently within a given three-dimensional cloud. It calculates the total heating and cooling functions at any point in a given PDR by adopting an escape probability method. It uses a HEALPix-based ray-tracing scheme to evaluate the attenuation of the far-ultraviolet radiation in the PDR and the propagation of the far-infrared/submm line emission out of the PDR. We present benchmarking results and apply 3D-PDR to i) a uniform-density spherical cloud interacting with a plane-parallel external radiation field, ii) a uniform-density spherical cloud interacting with a two-component external radiation field, and iii) a cometary globule interacting with a plane-parallel external radiation field. We find that the code is able to reproduce the benchmarking results of various other one-dimensional numerical codes treating PDRs. We also find that the accurate treatment of the radiation field in the fully three-dimensional treatment of PDRs can in some cases leads to different results when compared to a standard one-dimensional treatment.


I N T RO D U C T I O N
Photodissociation regions (PDRs; also known as photon-dominated regions) are ubiquitously present in the interstellar medium (ISM), and consist of predominantly neutral gas and dust illuminated by far-ultraviolet (FUV) radiation (6 < hν < 13.6 eV). Studies of PDRs allow us to understand the effects of FUV photons on the chemistry and structure of the neutral ISM in galaxies, as well as diagnosing the conditions within star-forming regions. PDRs are responsible for most of the infrared radiation from galaxies. The FUV photons usually arise from massive stars creating H II regions but sometimes from active galactic nuclei (AGN), which produce strong ultraviolet (UV) and X-ray emission.
Numerical models of PDRs have been around for about 30 yr or so and they have now evolved into complex computer codes accounting for a large number of physical and chemical effects (see review by E-mail: tb@star.ucl.ac.uk Röllig et al. 2007, hereafter R07). In particular, the past decade has seen a proliferation of codes capable of treating the chemistry and thermal balance within PDRs, each developed with distinct interests in mind. Some codes aim to encapsulate the detailed microphysics that describe the chemical and thermal processes at work in the gas and on the grains (e.g. Ferland et al. 1998;Abel et al. 2005Abel et al. , 2008Le Petit et al. 2006, 2009, while others focus on treating the gas-grain chemistry or other specific processes in detail whilst approximating others in order to explore large regions of parameter space (e.g. Röllig et al. 2006;Wolfire et al. 2008;Hollenbach et al. 2009). Some have been developed to model specific source structures with either one-plus-one-dimensional or fully two-dimensional geometries, including discs and outflow cavities around protostars (e.g. Kamp & van Zadelhoff 2001;Bruderer et al. 2009;Woitke, Kamp & Thi 2009); further departures from simple one-dimensional slab or spherical geometries have been pursued in models that treat PDRs as ensembles of discrete clumps, described by size and mass distribution functions, embedded within an interclump medium (e.g. Cubick et al. 2008;Kramer et al. 2008).

3D-PDR: a new 3D code for treating PDRs 2101
Efforts have also been made in modelling detailed microphysics in dynamically evolved simulations. Glover & Mac Low (2007a,b) and Dobbs et al. (2008) have included the formation of H 2 in threedimensional simulations of cloud formation. Improvements of these methods have been made by Glover et al. (2010) to additionally model the CO formation thus paving the way towards the comparison between simulations and observations.
Observations of atomic fine structure and molecular lines from PDRs have improved with the advent of infrared [e.g. Infrared Space Observatory (ISO) and Herschel] and submillimetre [e.g. James Clerk Maxwell Telescope (JCMT) and IRAM] telescopes, while models have been benefiting from increasingly accurate laboratory and theoretical data. Most models feature plane-parallel geometry, illuminated on one or both sides. This simplifies the radiative transfer problem because the illumination comes from one side only and hence only one line of sight needs to be taken into consideration (Flannery, Roberge & Rybicki 1980;R07 and references therein). In addition, some models use a spherical geometry where an isotropic FUV irradiation is taken into consideration. R07 provide a detailed account of the differences between plane-parallel and spherical models and underline the many assumptions and approximations implicit in both geometries, in particular when it comes to the treatment of the attenuation due to dust.
Neither of these two approaches can however deal with more complex geometrical issues such as, for example, clumpiness inside the clouds, multiple sources of radiation and non-spherical geometries of H II regions or galaxies. For extragalactic sources in particular, often unresolved at far-infrared (FIR) wavelengths, accounting for the total emission in fine structure and molecular lines, as well as from dust, requires the modelling of ensembles of star formation complexes where the geometrical issues raised above become important. Their modelling is best achieved by simultaneously modelling the observed spectra and structures of H II and PDR complexes. Hence three-dimensional integrated photoionization and PDR codes are essential for interpreting the wealth of data available for star-forming galaxies. However, no code currently offers a complete three-dimensional treatment of both ionized and PDR regimes (e.g. with realistic radiation fields and geometries and multiple exciting sources). Such a self-consistent code is particularly important for the modelling of external galaxies where the available angular resolution is not high enough to disentangle the different gas components. While three-dimensional gas and dust photoionization codes exist, such as MOCASSIN (Ercolano et al. 2003;Ercolano, Barlow & Storey 2005), a fully three-dimensional code for PDRs that can handle the gas-grain chemistry as well as the thermal balance is still lacking.
In this paper we present a development of the UCL_PDR code (Bell et al. 2006) to treat three-dimensional structures of arbitrary density distribution. The UCL_PDR code is an already benchmarked (R07) one-dimensional time-and depth-dependent gas-grain PDR code which includes time-varying density and radiation profiles. Our new code, 3D-PDR, 1 adopts the same features in modelling the chemistry as UCL_PDR does. It is a starting point towards the implementation of an integrated code which will treat dust, photoionized gas and PDRs together in fully three-dimensional computational domains. The integrated code aims to couple 3D-PDR and MOCASSIN, with the latter treating the propagation and attenuation of the UV radiation field as realistically as possible, including a detailed spectral energy distribution (SED) profile.
The paper is organized as follows. In Section 2 we discuss the numerical treatment, giving an overview of the code, our ray tracing scheme, treatment of the escape probability method, treatment of gas cooling and gas heating, the thermal balance and convergence criteria, as well as the approximations and assumptions we made. In Section 3 we present angular and spatial resolution tests for the requirements of our ray tracing scheme, and we benchmark our code against various one-dimensional codes discussed by R07. In Section 4 we show three examples to demonstrate the capabilities of our code in simulating one-or two-component UV fields in uniform or arbitrary density distributions. We discuss our conclusions in Section 5.

Overview
The 3D-PDR code uses the chemical model features of the fully benchmarked one-dimensional code UCL_PDR (Bell et al. 2006). It solves the chemistry and the thermal balance self-consistently within a given three-dimensional cloud of arbitrary density distribution. We note that 3D-PDR has the ability to solve any one-, two-or three-dimensional structure, however, in this paper we will present calculations from one-and three-dimensional PDRs only. The code uses a ray tracing scheme based on the HEALPIx package (see Section 2.2) to calculate the total column densities and thus to evaluate the attenuation of the FUV radiation into the region (see Section 2.3), and the propagation of the FIR/submm line emission out of the region. An iterative cycle is used to calculate the cooling rates (see Section 2.5) using a three-dimensional escape probability method (see Section 2.4), and heating rates (see Section 2.6). At each element within the cloud, it performs a depth-and time-dependent calculation of the abundances for a given chemical network (see Section 2.7) to obtain the column densities associated with each individual species. The iteration cycle terminates when the PDR has obtained thermodynamical equilibrium, in which the thermal balance criterion is satisfied (see Section 2.8), i.e. the heating and cooling rates are equal to within a user-defined tolerance parameter.
In Section 2.9 we present the approximations and assumptions we made for the three-dimensional treatment of PDRs, and in Appendix A we present a flowchart of the computational scheme used in 3D-PDR.

Ray tracing
The three-dimensional ray tracing scheme we use for calculating the column densities and the attenuation of the FUV radiation field is based on the HEALPIx algorithm (Górski et al. 2005). HEALPIx has been used in the past for similar purposes (i.e. Abel & Wandelt 2002;Alvarez, Bromm & Shapiro 2006;Abel, Wise & Bryan 2007;Krumholz, Stone & Gardiner 2007;Bisbas et al. 2009;Clark, Glover & Klessen 2012). It creates a set of N = 12 × 4 pixels uniformly distributed over a unit celestial sphere, where is the level of refinement. Each of these pixels represents the end of a vector (hereafter 'ray') emanating from the centre of a Cartesian coordinate system. A pixel defines the centre of an approximately square element of solid angle = 4π/N . Consider a cloud with an arbitrary density distribution consisting of N elem elements. Let p(x, y, z) be a random element from which Dashed lines show the extent of solid angle . Black dots are the evaluation points. Dot-dashed lines show the extent of the search cone which has as vertex the kth evaluation point and apex angle 2θ crit (0 < θ crit ≤ π/2 is user defined). The projection of an element p i on the HEALPIX ray will be taken if φ k i−1 ≡ p i k i−1 q ≤ θ crit , where q is the HEALPIX pixel, creating a new evaluation point. Every new evaluation point defines the vertex of the new search cone which however keeps the same apex angle in the sense that the cone 'moves' in parallel as we walk along the HEALPIX ray. the HEALPIX rays are emanated all over the computational domain. Thus the cloud will be divided into subvolumes, the elements of which belong to different rays of solid angle . For evaluating the integrations along each of these rays we create a discrete set of points which we dub 'evaluation points'.
The evaluation points are created by projecting the elements of the cloud which are closest to the line of sight of a specific ray (see Fig. 1). Similar schemes for creating evaluation points have been used also by other workers (i.e. Kessel-Deynet & Burkert 2000;Dale, Ercolano & Clarke 2007). The steps we follow are described below.
(i) We first sort all N elem elements with increasing distance from the element p using a HEAPSORT algorithm.
(ii) For a random element p and using the ANG2PIX_NEST HEALPIX subroutine, we find the pixel q in whose solid angle the element p belongs to.
(iii) If the angle φ p ≡ p pq is φ p ≤ θ crit , where 0 < θ crit ≤ π/2 is a user defined critical search angle, we take the projection of p on to the respective ray. This projection is the evaluation point k i and we assign it with identical properties to those of the projected element p . The search angle θ crit defines a search cone with apex angle 2θ crit and solid angle crit = 2π(1 − cos θ crit ) whose vertex is, in this case, the element p.
(iv) The evaluation point k i defines the next angle, φ k i , between another element, p , and the pixel, q, i.e. φ k i ≡ p k i q. k i also defines the new vertex of the search cone, which however keeps the same, 2θ crit , apex angle value in the sense that the cone 'moves' in parallel as we walk along the ray.
We repeat the above steps until we reach the end of the ray, for every HEALPIX ray; and for all N elem elements considering that each one of them is a HEALPIX source. We store all of these hierarchies in memory.
For a given spatial function f (r) (e.g. number density of some chemical species), we perform integrations along each ray by adopt-ing the trapezoidal rule: Here, L is the length of the ray, N eval is the total number of evaluation points along this ray and r k is the distance of the evaluation point k from the element p which defines the HEALPIX source. The distance |r k − r k−1 | is equivalent to the spatial distance between two evaluation points of the same ray which in turn defines the integration step dubbed as an 'adaptive step'. The advantage of this ray tracing scheme is that it can be applied directly to both grid-based and smoothed particle hydrodynamics (SPH) data without the necessity of implementing further modifications.

Treatment of the ultraviolet radiation field
For the purposes of this paper and for the current version of 3D-PDR, we simplify the treatment of the UV field; we neglect the contribution of the diffusive radiation by invoking the on-the-spot approximation (Osterbrock 1974). We do this in order to explore the effects introduced when one moves to a three-dimensional treatment of PDRs. As mentioned in the Introduction, our future plans include the coupling of 3D-PDR with MOCASSIN (Ercolano et al. 2003(Ercolano et al. , 2005 in a single integrated code in order to include a realistic treatment of a three-dimensional UV radiation field, and that is including the diffusive radiation component. In the current version of 3D-PDR the user is able to choose between three types of UV radiation field: plane-parallel (UNI); radial-sampling (ISO) or a field emitted spherically symmetrically by a point source (PNT); or any combination between these.
For a random element p(x, y, z), the strength of the incident radiation, χ (p), measured in units of the Draine (1978) interstellar radiation field, is calculated using the equation 2 where χ 0 (θ , φ) is the magnitude of the unattenuated field strength (Draines) at the surface of the cloud in direction θ and φ, and τ UV = 3.02 is a dimensionless factor converting the visual extinction to UV attenuation. The integration over solid angle ( d ) is approximated using a summation over all N HEALPIX rays. The term A V is the visual extinction defined as where A Vo = 6.29 × 10 −22 mag cm 2 . The integration ( L 0 n H dr) corresponds to the column of the H-nucleus number density n H along the ray of length L. This integration is approximated using a summation over the N eval evaluation points of the q HEALPIX ray and where r = |r k − r k−1 | is the adaptive step as discussed in equation (1).

Escape probability
Suppose that the random element, p (x, y, z), of the cloud is considered as a source of radiation. Assuming statistical equilibrium, the level populations and radiation field are related by where the summation is over the total number of levels included and n i (p), n j (p) are the populations of levels i, j, respectively. The lefthand side describes emission and the right-hand side absorption.
where A ij and B ij are the Einstein coefficients, C ij (p) is the collisional rate for excitation (i < j) and de-excitation (i > j) and J ij (p) is the mean integrated intensity received at p from all solid angles d . In our models, we adopt the large velocity gradient (LVG) or escape probability formalism (Sobolev 1960;Castor 1970;de Jong, Dalgarno & Chu 1975;Poelman & Spaans 2005) to describe the mean radiation field J ij (p) as The term S ij (p) is the source function due to transitions between levels i, j and is given by where ν ij is the photon frequency, and g i , g j are the statistical weights of n j , n i , respectively. The term B(ν ij ) is the total background radiation valid at FIR and submm wavelengths, including cosmic microwave background blackbody emission at T CMBR = 2.7 K and dust emission approximated as a modified blackbody at T dust .
The term β ij (p) describes the probability that a photon of frequency ν ij escapes from the element p without interacting with the rest of the cloud. In the present work we consider wavelengths in the FIR and submm range and we therefore neglect the absorption due to dust, considering only the line absorption component. We adopt the analytical expression for the escape probability β ij developed by de Jong et al. (1975): where τ L ≡ τ ij (p, q) is the line optical depth at the element p along the direction q, given by the expression where the integration is performed between the positions r 1 and r 2 which define the direction q, and u(p) is the root-mean-square of the thermal and turbulent velocities (see below). We approximate the integral over all space of equation (8) using the HEALPIX ray scheme described in Section 2.2 as Thus for the element p, the numerical expression for the escape probability is where q represents each individual HEALPIX ray. This equation provides the total escape probability corresponding to the summation of all individual escape probabilities per HEALPIX direction averaged over the total number of these directions. The numerical expression of equation (9) for the line optical depth, τ ij (p, q), is where the summation is over all N eval evaluation points along the q ray, r = |r k − r k−1 | is the adaptive step and u(p) is given by where k B is the Boltzmann constant, m H is the proton mass, T(p) is the gas temperature of element p and v TURB is the user-defined turbulent velocity. The LVG or escape probability method described here is used as an approximation to describe the mean radiation field intensity. Our future plans include the treatment of an exact line transfer from one region to another inside the PDR.

Treatment of gas cooling
The gas in molecular clouds is cooled primarily by the collisional excitation and subsequent emission of a number of key atomic and molecular species. Within PDRs this is usually dominated by emission in the [C II], [O I] and [C I] fine-structure lines and in the rotational transitions of CO. The cooling rates due to emission from these species are determined by the 3D-PDR code at every element within the cloud by calculating the emissivity (in erg s −1 cm −3 ) of each transition, having solved for the non-local thermodynamic equilibrium (LTE) excitation and radiative transfer under the LVG assumption, as described in Section 2.4. The total radiative cooling rate at each position within the cloud is then the sum of the emissivities of all radiative transitions considered. We include transitions between the ground-state fine structure levels of O ( 3 P 2 , 3 P 1 , 3 P 0 ), C ( 3 P 0 , 3 P 1 , 3 P 2 ) and C + ( 2 P 1/2 , 2 P 3/2 ) and 11 rotational levels of CO (note that up to 40 rotational levels can be treated in the code, but we limit the number here to increase computational speed for these initial models). Collisional excitation rates are taken from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005) for all available collision partners, namely H, He, H 2 , H + and e − for O and C; H, H 2 and e − for C + and H, He and H 2 for CO. The collisional rates at the required gas temperature are determined by linear interpolation between the fixed temperature values specified within these data files.
Deeper within molecular clouds, other molecular species, including the isotopologues of CO, OH and H 2 O, become important coolants, but despite being included in the UCL_PDR code we neglect their contribution within the 3D-PDR code since we are concerned with the thermal balance near the surfaces of these clouds, where gas temperatures show the greatest variations. At high visual extinctions, where these other molecular coolants become important, the gas and dust temperatures generally tend to a rather constant 8-15 K (e.g. fig. 12 in R07) in the absence of embedded heating sources and the appropriate dark cloud chemistry asserts itself, with little sensitivity to these small temperature variations.
At high densities (>10 6 cm −3 ), collisions with dust grains can also efficiently cool or heat the gas, depending on the temperature difference between the gas and the grains. This mechanism is also accounted for in the 3D-PDR code, following the treatment of Burke & Hollenbach (1983) and the accommodation fitting formula of Groenewegen (1994), assuming a standard Mathis, Rumpl & Nordsieck (1977, hereafter MRN) grain size distribution.
Cooling due to ro-vibrational emission of H 2 can also play a minor role in regulating the gas temperature close to the cloud surface, but we do not include this in the current version of the 3D-PDR code. At relatively high temperatures (>5000 K) neutral atomic gas can be efficiently cooled by excitation of the metastable 1 D levels of atomic carbon and oxygen. These processes are implemented in the code, but are disabled for the present study.

Treatment of gas heating
The mechanisms that contribute to the total heating of the gas and their treatment in the 3D-PDR code are identical to those in the UCL_PDR code, described in detail by Bell et al. (2006). Near the PDR surface, the dominant gas heating mechanism is the photoelectric ejection of electrons from small dust grains and polycyclic aromatic hydrocarbons (PAHs). We adopt the treatment of Bakes & Tielens (1994), where the total heating rate (in units of erg cm −3 s −1 ) is given by and G 0 is the local FUV flux expressed in units of the Habing (1968) field (which is related to the scaling factor for the Draine field by G 0 = 1.7χ ), n H is the total H-nucleus number density (cm −3 ), T is the gas temperature (K) and n e is the electron number density (cm −3 ). This heating rate is countered at high gas temperatures by cooling due to recombination of electrons with grains and PAHs, and we similarly adopt the analytical expression derived by Bakes & Tielens (1994) for this cooling rate: We assume the same standard MRN grain size distribution adopted by those authors in deriving these rates. We note that equation (14) describing the universal grain heating rate should be used for low intensities of radiation field such as those in the simulations presented here. In future updates of 3D-PDR the Weingartner & Draine (2001) heating rate will be used which gives better approximation at higher intensities of radiation fields. The collisional de-excitation of vibrationally excited H 2 following FUV pumping can also be an important heating mechanism in dense gas near the cloud surface. We assume a single excited pseudo-vibrational level of H 2 , denoted H * 2 , to effectively account for the full distribution of H 2 molecules in vibrationally excited levels and we describe our treatment of its formation and destruction in the next section. We adopt the associated heating rate for collisional de-excitation of H * 2 from Tielens & Hollenbach (1985): where n(H), n(H 2 ) and n(H * 2 ) are the number densities of H, H 2 and H * 2 , respectively (cm −3 ), E * is the energy of the single excited pseudo-vibrational level (2.6 eV) and γ H * 0 and γ H 2 * 0 are the collisional de-excitation rate coefficients (in units of cm 3 s −1 ) from the excited to the ground vibrational level for H and H 2 , given by We note that more complicated techniques for the H 2 treatment have been implemented and modelled showing differences in the heating rate of up to ∼1 dex (see R07).
In addition, photoionization of neutral carbon liberates about 1 eV per photoelectron (Black 1987, the rate of carbon photoionization is described in the next section). The internal energy of newly formed H 2 as it leaves the grain surface can also make a nonnegligible contribution to the gas heating near the PDR surface; following Black & Dalgarno (1976), we assume that the 4.48 eV of internal energy is distributed roughly equally between translation, vibration and rotation, so that 1.5 eV will go into kinetic energy per H 2 molecule formed. Deeper within the cloud, cosmic rays and turbulence (Black 1987) do most of the heating of the gas, with exothermic reactions also playing a minor role. Following Tielens & Hollenbach (1985), we assume a cosmic ray heating rate of where ζ is the cosmic ray ionization rate per H 2 molecule. For the rate of gas heating due to dissipation of supersonic turbulence, we assume (Black 1987) where v TURB is the turbulent velocity (km s −1 ) and l TURB is the turbulent scale length (assumed here to be 5 pc).

Model chemistry
The code determines the relative abundances of a limited number of atomic and molecular species at each cloud element in the model by solving the time-dependent chemistry of a self-contained network of formation and destruction reactions. The chemical network is a subset of the most recent UMIST data base of reaction rates (Woodall et al. 2007), consisting of 320 reactions between 33 species (including electrons), and includes photoionization and photodissociation reactions in addition to the standard gas-phase chemistry. We make use of the XDELOAD tool (kindly provided by L. Nejad) to construct the set of ordinary differential equations (ODEs) describing the formation and destruction of each species and the associated Jacobian matrix that together are used to compute the chemical abundances in the 3D-PDR code. We order the species in the ODE system according to their number of formation reactions, which has been shown to significantly speed up the computation of the abundances (Nejad 2005). In this paper we restrict our models to produce steady-state abundances by assuming a chemical evolution time of 100 Myr, sufficient for all reactions to reach equilibrium. However, the code is capable of following the full time-dependent evolution of the chemistry within a cloud, making it a powerful tool for modelling dynamically evolving structures. In our full chemical network, we include reactions involving vibrationally excited molecular hydrogen, whose internal energy can provide a means to overcome the activation barrier of certain neutral-neutral and ion-molecule reactions, thus considerably enhancing the abundances of some species, such as CH + . Significant abundances of molecular hydrogen in vibrationally excited levels can be maintained in PDRs due to the FUV pumping of H 2 to electronically excited states, followed by its decay to populate the vibrational levels of the ground electronic state. In our models, we adopt a simplified treatment of this mechanism following Tielens & Hollenbach (1985) and others in assuming that all vibrationally excited H 2 is in a single characteristic vibrational level, with effective spontaneous emission and collisional de-excitation rates, that approximates the full distribution amongst all vibrational levels. Molecular hydrogen in this pseudo-level is labelled H * 2 and the reactions that form and destroy it are taken from Tielens & Hollenbach (1985), with updated rates for some reactions taken from Agúndez et al. (2010).
We also include reactions involving neutral and singly ionized PAHs, adopting the rate coefficients proposed by Wolfire et al. (2003Wolfire et al. ( , 2008, though we omit these reactions for the reduced chemistry of the benchmark comparisons and test applications described in this paper. More detailed models to be presented in forthcoming papers will include the full reaction network, including the role of PAHs and vibrationally excited H * 2 . The formation of molecular hydrogen on grain surfaces is a key reaction in determining the properties of PDRs since, together with destruction through photodissociation (described below), it governs the transition from atomic H to H 2 within the PDR and therefore plays a critical role in the chemistry (many species form through reactions with H 2 ) and thermal balance (through heating processes such as those involving vibrationally excited H 2 , and as a collision partner for coolant species and a coolant in its own right). We have implemented the treatment of Cazaux & Tielens (2002 who have shown that the process can be described using a rate equation formalism that adequately accounts for both chemisorbed and physisorbed hydrogen atoms reacting on grain surfaces. We use their standard values for the properties of both graphitic and silicate grains to determine the formation efficiency, and the expression for the sticking coefficient of H atoms on grains from Hollenbach & McKee (1979). We describe the grains by their global population properties, assuming a standard MRN size distribution and graphite/silicate composition. We do not consider the level-specific distribution of newly formed H 2 leaving the grains, instead treating it as one species. For the benchmark tests described in Section 3.3, we have use a simplified formation rate (as described in that section) in order to better match the benchmark model specifications.
With the exception of the reaction rates for the photodissociation of H 2 and CO and the photoionization of carbon, all photoreaction rates are calculated using the standard UMIST treatment adapted to our HEALPIX-based ray tracing scheme, where the total rate at a given position in the cloud is the sum of the rates determined along each HEALPIX ray to the PDR surface. The rate (s −1 ) for a photoreaction i along a particular ray q is then given by where α i is the unattenuated photoreaction rate (s −1 ), evaluated for isotropic illumination by the standard Draine interstellar radiation field (Draine 1978), χ 0 (q) is the FUV flux incident on the PDR surface at the element intersected by ray q and specified by a scaled equivalent of the Draine field, A V (q) is the total visual extinction along that ray to the surface and γ i is a scaling factor that relates the attenuation in the visible to that in the FUV. The photodissociation rates of H 2 and CO depend sensitively on the column densities of these species along the direction of the incident FUV radiation, since their absorption of FUV photons leads to significant shielding against photodissociation. This (self-)shielding is therefore treated explicitly in the code, using the results of the detailed calculations of Lee et al. (1996) and van Dishoeck & Black (1988) for H 2 and CO, respectively. Those authors ran full radiative transfer models accounting for both self-shielding and line overlap of H, H 2 and CO lines in order to determine the degree of shielding produced under a range of cloud conditions. They found that these detailed processes could be well described by shielding functions that depend only on the total H 2 and CO column densities to the PDR surface. We have therefore adopted these treatments and use the tabulated shielding functions provided in their papers. We do not explicitly account for state-specific photodissociation of H 2 and CO, since the rates and self-shielding factors listed by Lee et al. (1996) and van Dishoeck & Black (1988) are strictly valid for the global populations of these molecules, with the exception of our inclusion of vibrationally excited H * 2 in a single excited pseudo-level, for which we adopt the effective photodissociation rate given by Röllig et al. (2006). We also neglect detailed treatment of ro-vibrational cascades following electronic excitation by the UV photons, since such a treatment would dramatically increase the computational time without significantly altering the resulting abundances. In addition, we account for the shielding of neutral carbon against photoionization using the treatment of Kamp & Bertoldi (2000). In all three cases, the column densities of H 2 , CO and C needed to calculate the shielding factors at a given point in the cloud are determined along each HEALPIX ray to the PDR surface.
The calculation of the H and H 2 abundances at the elements near the cloud surface depends critically on the shielding provided by H 2 against its own photodissociation, which, taken together with H 2 formation on grains, represents the main formation/destruction cycle for molecular hydrogen in the UV-illuminated gas. The amount of self-shielding along a given ray to the cloud surface is itself sensitive to the total column density of H 2 and therefore requires that the H 2 abundance be known at each evaluation point along the ray. A similar relation links the photodissociation of CO to its column density along each ray.
There exists, then, an interdependence between the photodissociation rates and the abundances of all elements near the cloud surface, requiring that the abundances be calculated and the resulting column densities updated a number of times before correct values can be obtained for both. We therefore perform a chemistry iteration each time that the gas temperatures are changed, in which the reaction rates and shielding factors are determined at the new temperature and for the current column densities, new abundances are calculated, the column densities are updated based on the new abundances and the process is repeated for I CHEM iterations. After numerous tests for convergence, we find that between five and ten chemical iterations are needed at the start of the code in order to correctly determine the abundances of H and H 2 , and C + , C and CO, near the surface. Following this first determination of the chemistry at the initial 'guess' temperature (see Section 2.8), we find that subsequent changes to the gas temperature require only one or two iterations of the chemistry to reach convergence. For the models presented in this paper, we have performed I CHEM = 8 chemical iterations at the start of the code and then three chemical iterations after each change to the gas temperature.
The capability exists within the 3D-PDR code to include the role of grains in the chemistry, including freeze-out of gaseous species on to grains, surface reactions and the release of grain mantle species back into the gas phase by means of evaporation, photodesorption or desorption due to cosmic ray heating. However, inclusion of the full gas-grain chemistry can dramatically increase the computational time needed to follow the evolving abundances, so reactions involving dust grains are currently omitted from the models in order to increase the speed of the code and to be able to compare with other PDR models that excluded grain chemistry for the benchmarking tests of R07.

Thermal balance and convergence criteria
3D-PDR starts the calculations by setting up an initial 'guess' temperature profile which is used in order to begin the iterative process for level populations and thermal balance convergence. The implicit assumption of an initially uniform temperature profile over the entire PDR of arbitrary density distribution may lead to unphysical values of level populations in certain parts of the PDR causing the iterative process to fail (i.e. level populations cannot converge). To avoid this we ran several one-dimensional calculations of uniform density PDRs interacting with several FUV field strengths. The density of PDRs in these calculations spans 10 2 ≤ n (cm −3 ) ≤ 10 6 and the field strength spans 1 ≤ χ (Draines) ≤ 10 6 . We find that the equation provides an acceptable estimate for an initial 'guess' temperature profile (first iteration over thermal balance) in order to begin the overall iterative process. In this equation, T guess (K) is the temperature and χ (Draines) is the attenuated FUV field strength.
Using the temperature profile of equation (24) the code begins the iteration process in order to obtain level population convergence. We assume that the level populations have converged when the change in any given population between two consecutive iterations is less than a user-defined tolerance parameter (in this paper σ err < 1 per cent). We then calculate the total cooling ( ) and heating ( ) rates. By comparing these rates we assign new gas temperatures i.e. if > γ a lower temperature than T guess is required and if < γ a higher temperature than T guess is required. The technique we use is the following. For a cloud element p, if its temperature changes monotonically from that given by equation (24), the new temperature used in the next iteration over thermal balance will differ by 30 per cent from its current value. If the change is not monotonic, then a binary chop routine is performed between the current temperature and the one obtained in the previous iteration on thermal balance.
Once the new kinetic temperatures have been calculated and updated, we start iterating again over the level populations. We repeat this process (i.e. assigning new temperatures, iterating over level populations, etc.) until we reach convergence over thermal balance. We assume we have obtained thermal balance convergence when the heating and cooling rates are equal to within some user-defined error tolerance (in this paper σ err,T ≤ 0.5 per cent) or when the change in gas temperature between iterations is negligible, i.e. smaller than a user-defined, T diff , value (in this paper T diff ≤ 0.01 K).

Approximations and assumptions
In order to make tractable the problem of simultaneously calculating the chemistry, thermal balance and radiative transfer within a three-dimensional cloud, we have necessarily made a number of simplifying assumptions and approximations in the 3D-PDR code. We list here the most important of them.
(i) Assumption: the level populations and resulting emissivities change rapidly with respect to the changes in abundance due to the chemistry.
(ii) Approximation: the radiative transfer of the FUV radiation into the model cloud is treated by considering only the attenuation by dust and neglecting the line scattering/diffusive terms by invoking the on-the-spot approximation (Osterbrock 1974). The grain temperature is not affected by the UV radiation field.

Angular resolution and search angle
In this section we explore the values of the level, , of angular resolution and the search angle, θ crit , for which 3D-PDR evaluates integrations at reasonable accuracy along the HEALPIX rays. To do this we perform a test to measure the integration accuracy of the code in the evaluation of column density. We consider a spherically symmetric cloud, the centre of which defines the centre of the coordinate system. The radius of the cloud is R = 1 pc and its spatial H-nucleus number density profile is n H (r) = n 0 e −r/R , where r is the radial distance from the centre and n 0 = 100 cm −3 . Therefore, the column density, N, in any direction as seen from the centre is N = R 0 n H (r) dr 1.95 × 10 20 cm −2 .
For the construction of the cloud we use N elem = 10 5 elements uniformly distributed. We also assume that the hierarchy of rays is emanated from the element p 0 which is positioned at the centre of the sphere.
Since the escape probability function β ij (p) of the element p (equation 11) is the average value of the escape probabilities over all HEALPIX rays, it is important to know the integration accuracy of the present method in calculating a function averaged over all rays. For this reason, we define as the average value of column densities over all HEALPIX rays, where N q is the column density of the q ray calculated as discussed in equation (3). Fig. 2 shows the error between the analytical value, N, and the calculated N N of the column density versus θ crit and for different levels of refinement, . The solid line represents the errors at = 0 (N 0 = 12 rays), the dotted line at = 1 (N 1 = 48 rays), the short dashed at = 2 (N 2 = 192 rays) and the long dashed at = 3 (N 3 = 768 rays).
Overall we see that σ err 3 per cent even for θ crit as low as π/12, so the accuracy is acceptable. We note that by increasing N elem , the error σ err is decreased. We also observe that by increasing the level, , of angular resolution refinement the error σ err is decreased, however, at the expense of computational time since more evaluation points are created. For the tests and applications described here we will use = 0 and θ crit = 0.8 π/4.

Resolution along a ray
Here we examine the resolution requirements needed along a HEALPIX ray in order to establish when our calculations are converged. We use a one-dimensional cloud of uniform H-nucleus number density n H consisting of N elem elements. Since the density  For = 0 and θ crit ≥ π/6 we achieve reasonable integration accuracy and at low computational cost. is constant, from equation (3) the length L of the cloud will be given by the equation where we set A V,max = 10. The elements are aligned with two opposite HEALPIX rays, namely 3 the rays with ID = 4 and 6 of = 0. Consequently, the search angle criterion is eliminated and the elements pre-define the evaluation points. Considering that these two rays define the x-axis of a Cartesian coordinate system, we apply a UV field of strength χ from the −x side; that is from ray ID = 6. For the rest of HEALPIX rays we assign very high optical depths, implying that the one-dimensional line represents a threedimensional semi-infinite slab. We construct the cloud by creating elements logarithmically distributed along the x-axis. We use N A V elements per A V dex and with −5 ≤ log (A V ) ≤ 1. We run three different tests with n H = 10 2 cm −3 and χ = 1 Draine (T1); n H = 10 3 cm −3 and χ = 10 3 Draines (T2) and n H = 10 4 cm −3 and χ = 10 Draines (T3). We vary the number N A V and we measure the value of visual extinction A V,trans at which the H/H 2 transition occurs. Fig. 3 plots N A V versus the relative error σ err between the A V,trans and the respective value for N A V = 200. In all tests we find that we obtain convergence for N A V as low as 20. In addition, we find that the minimum value of visual extinction needed in order to resolve the H/H 2 transition at this error must satisfy the relation log (A V,min /A V,trans ) −1.

Comparison with the other PDR codes
In order to assess the reliability of the new code, we have run a series of models designed to reproduce the benchmark tests described in the R07 comparison study. In that paper, results from a number of the most widely used PDR codes were compared for a set of models in which the capabilities of all codes were restricted to an agreed upon (and much simplified) set of parameters and treatments for processes such as the chemistry, UV attenuation and heating rates. Table 1 lists the main physical parameters that varied between the four models considered; the reader is referred to the R07 paper for full details of the other parameters used in the models. We note that, whilst we have adopted the same elemental abundances and treatment for the attenuation of the UV radiation in our tests with the 3D-PDR code, a number of differences remain between our code and those included in the R07 comparison. In particular, we are using a slightly larger and updated chemical network with rates taken from the most recent release of the UMIST data base, and updated collisional excitation rates taken from the LAMDA data base. Most importantly, we are using a modified form of the standard rate of H 2 formation on grain surfaces (in units of cm 3 s −1 ), taken from de Jong (1977), that includes an additional exponential term to reduce the formation efficiency at high temperatures: As discussed below, this additional term leads to some minor differences between our results and those of the R07 benchmark models. Whilst this treatment has been superseded in recent years by more advanced formalisms that account for both chemisorption and physisorption and for mixed grain compositions (see e.g. Cazaux & Tielens 2004), we have chosen to neglect such treatments for the simplified models we present here.
To allow comparison with the R07 results, we have restricted the cloud geometry used in the models to that of a one-dimensional slab of constant density illuminated by a UV field from only one side, as described in Section 3.2. In addition, we have restricted the escape of photons from the cloud to a single ray directed towards the illuminated surface, thereby emulating a semi-infinite slab. The gas temperature is allowed to vary and is determined by solving the thermal balance explicitly (see Section 2), but we assume a fixed dust temperature of T dust = 20 K throughout the cloud. In all models, the cloud consists of N A V = 100 elements logarithmically distributed per A V dex with −4 ≤ log (A V ) ≤ 1.3 for the n H = 10 3 cm −3 density case (so N elem = 530) and with −6 ≤ log (A V ) ≤ 1.3 for the n H = 10 5.5 cm −3 density case (so N elem = 730).
The results of the 3D-PDR test models are compared to those of the R07 paper in Figs 4-6 (for discussion on Fig. 6 see Section 4.1). We use the workshop results 4 of the following codes: UCL_PDR (Papadopoulos, Thi & Viti 2002;Bell et al. 2005Bell et al. , 2006, CLOUDY (Ferland et al. 1998;Abel et al. 2005;Shaw et al. 2005), COSTAR (Kamp & Bertoldi 2000;Kamp & van Zadelhoff 2001), HTBKW (Tielens & Hollenbach 1985;Kaufman et al. 1999;Wolfire et al. 2003), KOSMA-τ (Stoerzer, Stutzki & Sternberg 1996;Bensch et al. 2003;Röllig et al. 2006), LEIDEN (Black & van Dishoeck 1987;van Dishoeck & Black 1988;Jansen et al. 1995), MEIJERINK (Meijerink & Spaans 2005), MEUDON (Le Bourlot et al. 1993;Le Petit, Roueff & Le Bourlot 2002;Le Petit, Roueff & Herbst 2004) and STERNBERG (Sternberg & Dalgarno 1989Boger & Sternberg 2006). Fig. 4 shows the results for benchmark model V2 (n H = 10 3 cm −3 ; χ = 10 5 Draines), including the gas temperature profile, number densities of H, H 2 , C + , C and CO and emergent intensities (surface brightnesses) of the dominant cooling lines. As can be seen, the overall agreement is very good, with the results from the 3D-PDR code typically falling within the scatter of results produced by the other codes. In addition to the results from the R07 comparison study, we have also obtained results for the four benchmark models using the latest version of the UCL_PDR code by adopting identical physical parameters and the same chemical network as used in the 3D-PDR code. The results for model V2 are included in Fig. 4 (labelled as UCL_PDR11) and show excellent agreement between the one-dimensional and three-dimensional versions of the UCL_PDR code. The results are similarly identical for the other three benchmark models and are therefore not shown in the remaining figures.
The only notable differences between the 3D-PDR and the R07 results visible in Fig. 4 are the much lower H 2 abundance at the outer cloud edge, which is due to the reduced H 2 formation efficiency at high temperatures, as discussed above, and the rise in neutral carbon abundance (and corresponding drop in C + abundance) at lower A V in the 3D-PDR model, which is due to the more advanced treatment of the carbon photoionization rate that we have adopted, including shielding by lines of H 2 and C (see Section 2.8). This difference is also reflected in the [C I] local emissivity profile. Fig. 5 shows selected comparisons of the results for benchmark models V3 and V4. Although we state that the benchmarking model V3 of R07 should not be considered as a potential PDR (significant contrast between the high density and the low radiation field which leads to a temperature of ∼20 K at A V 0.1 mag), we include it in the present work as we are able to compare 3D-PDR with the other codes even under such extreme conditions. While the scatter in the range of results from all the codes is larger, the results obtained by our code nevertheless continue to show very good agreement with the rest of the codes, the only differences of note coming from the different treatments for the H 2 formation and carbon photoionization rates already discussed.
Overall, the results of these four benchmark tests demonstrate that the 3D-PDR code compares very favourably with the well-established PDR codes included in the R07 study when using similarly limited chemical networks and treatments of the various microphysical processes. We therefore consider the new code to be reliable and in the next section we go on to demonstrate some of the more advanced model applications that become possible with the fully three-dimensional geometry offered by the 3D-PDR code.
In addition, we note that 3D-PDR is approximately two times faster than UCL_PDR for running one-dimensional models. The primary source of this significant speed-up is the usage of pre-defined evaluation points which are either user specified (i.e. in the case of one-dimensional runs of 3D-PDR) or created automatically due to the projection of the elements that make up the cloud (as described in Section 2.2, i.e. in the case of two-or three-dimensional runs). Thus the evaluation points act as a fixed grid, in contrast with the adaptive grid used in UCL_PDR. However, techniques controlling an adaptive increase of resolution inside PDRs are inevitably necessary for 3D-PDR in treating complex three-dimensional structures. We plan to implement and examine these techniques in a forthcoming paper.

A P P L I C AT I O N S
In this section we present three different applications to demonstrate the capabilities of our code in simulating three-dimensional cloud structures. We explore (i) a uniform-density spherical cloud interacting with a plane-parallel external radiation field (Section 4.1), (ii) a uniform-density spherical cloud interacting with a two-component external radiation field (Section 4.2) and (iii) a cometary globule interacting with a plane-parallel external radiation field (Section 4.3). In all these applications we use = 0 levels of HEALPIX refinement and θ crit = 0.8 π/4 rad. The turbulent velocity is set to v TURB = 1 km s −1 . The dust temperature is fixed and set to T dust = 20 K. The cosmic ray ionization rate is set to ζ = 5 × 10 −17 s −1 . In addition we use I CHEM = 8 iterations over chemistry at the beginning of each simulation and I CHEM = 3 during each new iteration over thermal balance (see Section 2.7).
We note that the applications presented here are simplified examples which demonstrate however the capabilities of 3D-PDR in modelling any kind of density structure under the interaction of a UV radiation field.

Interaction of a uniform-density spherical cloud with a plane-parallel radiation field
In this application we consider a uniform-density spherically symmetric cloud with a H-nucleus number density of n H = 10 3 cm −3    Because of three-dimensional effects we find that the temperature of the limb is lower than that of the central region since in the latter case the UV field is impinging more radially. and a radius of R = 5.15 pc. The cloud has therefore radial visual extinction of A V = 10 mag at its centre, assuming that the surface of the cloud is at A V = 0 mag. A plane-parallel uniform UV radiation field of strength χ = 10 Draines is impinging from one side. The density and UV field strength used for this application correspond to the parameters used in model V1 of R07.
We construct the sphere in the following way. Using HEALPIX we create N 4 = 3072 (level = 4) pixels uniformly distributed on the surface of the sphere and we take into account only those which lie on the hemisphere on which the UV field is impinging. Each one of these pixels defines the start of a line segment which penetrates the sphere; is parallel to the direction of the UV field and consists of elements logarithmically distributed filling up the entire sphere. We use N A V = 60 elements per A V dex with A V,min = 10 −5 mag which, as discussed in Section 3.2, ensures high resolution along the direction of the UV field. Thus the total number of elements is approximately N elem 5.38 × 10 5 .
Since this application is equivalent to model V1 in R07, in Fig. 6 we compare our results with those of the benchmarked codes; we find in general a very good agreement, particularly with the UCL_PDR code; we note however that our temperature values at high A V are noisier primarily due to additional three-dimensional effects and due to our different iteration criterion which leads to less smoothing; we also find that our H/H 2 transition occurs at slightly earlier A V . Fig. 7 shows how the temperature varies when considering the limb and the equator of the sphere separately; as expected the temperature is slightly lower in the regions located around the limb of the sphere (as seen from the UV field) in comparison with the temperatures obtained in regions around the equator of the sphere. This is because the radiation field is impinging more radially in the latter than in the former leading to small differences in the attenuated field strengths.

Interaction of a uniform-density spherical cloud with a two-component radiation field
In this application we consider a uniform-density spherically symmetric cloud with a H-nucleus number density of n H = 2 × 10 3 cm −3 and a radius of R = 2.58 pc. As in Section 4.1, the cloud has a radial visual extinction of A V = 10 mag at its centre assuming that the entire surface is defined as A V = 0 mag. A two-component UV radiation field was adopted. The first component corresponds to a radial sampling field of strength χ ISO = 120 Draines, and the second component corresponds to a plane-parallel radiation field of strength χ UNI = 2 × 10 3 Draines impinging from one side.
We constructed the sphere using a combination of two different arrangements of elements. In the first arrangement we used HEALPIX to create N 4 = 3072 (level = 4) pixels uniformly distributed on the surface of the sphere. Each one of these pixels, along with the centre of the sphere, define ray segments which we constructed using N A V = 20 elements logarithmically distributed per A V dex and with A V,min = 10 −5 mag. Thus we created a sphere with approximately N elem,1 3.7 × 10 5 logarithmically radially spaced elements. In the second arrangement we create a sphere of N elem,2 = 3 × 10 4 uniformly distributed elements. Therefore the resultant combined sphere consists of approximately N elem = N elem,1 + N elem,2 4 × 10 5 elements and possesses both an approximately uniform distribution in its inner part and a logarithmic distribution in its outer parts, ensuring that the resolution requirements described in Section 3.2 are met. . At the bottom left we plot a cross-section of the cloud, showing the gas temperature. We see that the temperature at the surface of the right-hand hemisphere is ∼270 K due to the radial sampling component of the radiation field, and the temperature at the lefthand hemisphere reaches 400 K due to the additional interaction of the plane-parallel component of the radiation field. The local undulations observed here do not correspond to local differences in cooling and heating but instead are a result of numerical noise introduced by the discretization in angle of the = 0 choice of HEALPIX rays. Although a selection of > 0 values would smear out these undulations, it would increase the computational cost without offering a significant improvement in the analysis.
Overall, we see that the PDR is 'squeezed' at the left-hand side due to the plane-parallel radiation field. This is seen even for the CO (1-0) line which is embedded in the inner part. Converting to units of erg s −1 cm −2 sr −1 , we find that the strongest coolant is the [O I] 63 μm line, which on average is ∼1.7 times stronger than the second coolant, [C II] 158 μm. At the bottom right we show an RGB composite image of the CO (1-0) (red), [C I] (green) and [C II] (blue) emission maps. Here, we observe a well-defined stratification of species.
In this application we additionally explore how the calculations of the one-dimensional treatment of PDRs diverge from the corresponding calculations when a fully three-dimensional treatment is taken into consideration. To do this we perform a one-dimensional calculation of a PDR which has the same parameters of those explored above (i.e. n H = 2 × 10 3 cm −3 , N A V = 20 elements logarithmically distributed per A V dex and with A V,min = 10 −5 mag and A V,max = 10 mag and χ = 2120 Draines field strength impinging from one side). We then compare the attenuation of the UV field strengths and the number densities of C + , C and CO of this run and the corresponding values taken from two different radial directions from the spherical cloud. These directions are shown in the RGB composite image of Fig. 8. The dotted line (equator) is parallel to the direction of the plane-parallel UV field and the dashed line (diagonal) is not. In particular the equator corresponds to the HEALPIX ray ID = 6 of = 0 and the diagonal to the ray ID = 9, supposing that a 12-ray structure is emanated from the centre of the sphere. While the radial visual extinction of the sphere is quite high, we neglect the contribution of the radial sampling of the UV radiation impinging from the opposite side of the cloud. Fig. 9 shows the attenuation of the UV field (upper panel) and the number densities of C + , C and CO (lower panel) for the onedimensional calculation (solid line), along the equatorial direction (dotted line), and along the diagonal direction (dashed line). Although the unattenuated field strength at the surface of the PDR is the same in all cases, the difference due to the three-dimensional structure has an impact in the attenuation of the UV field as seen from the diagonal direction. On the other hand, due to the symmetry obtained, the equatorial direction is in a very good agreement with the one-dimensional calculation. We therefore find a difference in the distribution of the C + , C and CO abundances depending on the direction along which we perform an observation even in the present simplified case.

Interaction of a cometary globule with a plane-parallel radiation field
In this application we explore the capability of the code to simulate PDRs with arbitrary density distributions. To do this we used as initial conditions a snapshot taken from a SPH simulation presented in section 4.4 of Bisbas et al. (2009). That SPH simulation (using the SEREN code; Hubber et al. 2011) examined the interaction between an initially uniform density clump with UV radiation emitted spherically from an exciting source (by invoking the on-the-spot approximation; Osterbrock 1974) which was placed outside and far away from the clump (for full details see Bisbas et al. 2009). This interaction, referred to as 'radiation-driven implosion' (Sandford, Whitaker & Klein 1982;Bertoldi 1989;Lefloch & Lazareff 1994), drives a strong shock front into the inner part of the clump, creating a morphological structure reminiscent of the cometary globules observed in the interior of ionized regions (such as in the Helix Nebula; Matsuura et al. 2009), and which may trigger star formation (Kessel-Deynet & Burkert 2003;Gritschneder et al. 2009;Bisbas et al. 2011;Haworth & Harries 2012).
We take a snapshot at t = 0.12 Myr. A cross-section plot (at z = 0) of the density structure of the clump at t = 0.12 Myr is shown at the bottom left of Fig. 10 where the ionizing radiation is impinging from bottom to top. At this time the cloud has attained a V-shape structure which contains an approximately ellipsoidal core of density n H ∼ 10 5 -10 6 cm −3 located at its tip and directly exposed to the ionizing radiation. Although in the SPH simulation the radiation is emitted spherically symmetrically from the distant exciting source, the angular size of the clump is quite small (∼6 • ) and we therefore consider a plane-parallel radiation field here. The UV photon flux impinging is 2.18 × 10 9 cm −2 s −1 corresponding to a field strength of approximately χ = 30 Draines.
For the purpose of this application and for the sake of computational speed, we only perform calculations up to A V,max = 4 mag of visual extinction. For A V > 4 mag and for the radiation field strength and density considered for this application, the cloud cannot be treated as a PDR anymore as it has reached dark cloud conditions, i.e. the radiation field does not penetrate any longer and hence does not have any effect on the chemistry. Instead, the latter will be mainly dominated by cosmic-ray-induced reactions (independent of optical depth). For a proper treatment of dark cloud chemistry, depletion on to dust grains should be taken into consideration. The temperature has reached equilibrium values of ∼10 K;  CO lines are mainly saturated and hence not contributing much to the cooling, and no strong source of heating contributes to increasing the temperature. To avoid calculations in this dark molecular region, we use the following technique. For a random element, p, if the magnitude of A V of the HEALPIX ray with the highest attenuated flux exceeds a user-defined threshold (here A V,max = 4), then p is not considered as a PDR and calculations are omitted. Here, we also omit all regions with n H ≤ 100 cm −3 as these belong to the inner part of the H II region. The density structure shows some symmetry with no abrupt gradients, i.e. the density changes smoothly and there are no sharp density enhancements. From the emission maps we see that the species in the PDR are distributed smoothly and follow the density profile. The weakest emission is produced by [O I] 63 μm and the strongest by CO (1-0) implying that the molecular gas dominates over the atomic contribution. However, considering the transition frequencies for these maps, we find that the [O I] 63 μm and the [C II] 158 μm lines are the dominant coolants with the first being somewhat ( 10 per cent) stronger. In addition, the thickness of the PDR is very small comparing with the dimensions of the cometary globule. This is observed in the RGB composite image. Although there is a stratification of the species with [C II] emitted from the outermost part and CO (1-0) from the innermost part of the globule as expected, the transition between the species occurs in quite a thin layer. Since the colouring of these species overlap, we observe a bright white rim around the whole cometary globule.

C O N C L U S I O N S
We have presented 3D-PDR, a numerical code for simulating threedimensional PDRs of arbitrary density distribution. The code uses a ray tracing scheme based on HEALPIX in order to calculate the column densities and the escape probability in every direction for every element within the cloud. We adopted a reduced chemical network of 33 species and 320 chemical reactions. Through an iterative process the code calculates the cooling and heating rates for every cloud element adjusting the gas temperature at each iteration in an attempt to balance the heating and cooling rates. The code terminates once the difference between the heating and cooling is negligible, i.e. when the PDR has obtained thermal balance.
We tested the ray tracing scheme in calculating the column density of a particular element against its known analytical expression and we found very good agreement. We have also explored the spatial resolution requirements for simulating a PDR and we found that our code resolves one-dimensional uniform-density PDR if it is constructed using N A V = 20 elements logarithmically distributed per A V dex.
Furthermore, we repeated the benchmarking tests presented in R07 and we compared our results with the ones obtained by other one-dimensional PDR codes. Overall we find very good agreement between the one-dimensional codes and 3D-PDR. In addition, we explored the capabilities of our code in simulating three-dimensional structures exposed in one-or two-component UV radiation fields. In particular we examined the following.
(i) We examined the interaction of a uniform-density spherical cloud with a plane-parallel radiation field in which the values of density and field strength were identical to those of model V1 in R07. We found very good agreement between 3D-PDR with the onedimensional codes and in addition we observed at low A V cooler temperatures in the limb of the sphere in contrast with the higher temperatures in the equatorial regions; an effect directly related to the three-dimensional treatment in our case.
(ii) We examined the interaction of a uniform-density spherical cloud with a two-component radiation field, consisting of a radial sampling field and a plane-parallel field. We explored the differences in results obtained when a fully three-dimensional treatment of the PDR is taken into consideration in contrast with a one-dimensional simplification. We found that the results differ according to the direction at which observations are performed.
(iii) We examined the interaction of a cometary globule with a plane-parallel radiation field, where we considered as initial conditions a snapshot taken directly from an SPH simulation. We found that the PDR location follows the density profile of the globule, i.e. the abundances of species change in agreement with the density structure. We also found that the thickness of PDR is quite small in comparison with the overall size of the globule using composite RGB emission maps. This application showed also the capability of our code to model any type of density and particle distribution.
The coupling of 3D-PDR with MOCASSIN will be presented in a forthcoming paper. The integrated code should make feasible a realistic treatment of three-dimensional H II/PDR complexes including a detailed treatment of SEDs, thus offering a powerful tool in studying such structures with arbitrary density distributions and multiple exciting sources.