Userguide

The complete userguide can be downloaded here.

DustEM is a numerical tool in fortran 95 jointly developed by IAS and IRAP to compute the extinction, the emission and the polarization properties of interstellar dust grains heated by photons. The dust emission is calculated in the optically thin limit (no radiative transfer) and the default spectral range is 40 to 10\(^8\) nm. The code has been designed so that dust properties can easily be changed and mixed and to allow for the inclusion of new grain physics (Compiègne et al. 2011). The grain temperature distribution \({\rm d}P/{\rm d}T\) is derived with the formalism of Désert et al. (1986) assuming that the grain cooling is continuous. To describe the long wavelength emission of dust, the initial algorithm has been revised to cover the low temperature part of \({\rm d}P/{\rm d}T\).

After a few words on how to set up DustEM, we detail below the way to control DustEM as well as the files it needs or produces.

1. Getting started
2. Input files
  2.1 Defining a dust model
  2.2 Template dust models
  2.3 The size distribution
  2.4 Grain physics
  2.5 Grain data
  2.6 Astrophysical data
3. Output files
4. DustEM Python

 

1. Getting started

After downloading the tar ball from here:

  • Unpack the archive with "tar xvzf dustem.tar.gz". This will generate the \(\tt /dustem\) directory which should contain the following subdirectories: \(\tt /data\), \(\tt /hcap\), \(\tt /oprop\), \(\tt /out\), \(\tt /py\), and \(\tt /src\). The fortran files are in the \(\tt /src\) directory, data files in the \(\tt /oprop\), \(\tt /hcap\), and \(\tt /dat\) directories and output files are in the \(\tt /out\) directory. All data files have headers (where each line begins with #) that document their content.
  • In \(\tt /src\), edit the file "DM_constants.f90" to define the string \(\tt data\_path\) as the absolute path to the /dustem directory you just created. If you are running DustEM coupled to the Meudon PDR code, you must also define the string \(\tt dir\_PDR\) as the absolute path to \(\tt /dustem/out\).
  • Edit "Makefile" and set up the standard commands corresponding to your fortran compiler (gfortran, g95 and ifc have been tested), e.g.:
    • FC = gfortran
    • FFLAGS = -O2 -fno-second-underscore
  • Type "make" to compile, possibly preceded by "make clean", if a former attempt went wrong. A successful compilation will end with \(\tt That's\; it!\). The executable is called \(\tt dustem\).
  • From within \(\tt /src\), the run command is "./dustem". This runs the model of Compiègne et al. (2011) and should take a few seconds CPU on current computers.
  • After running DustEM, you can look at the output: a suite of useful Python routines to visualize or generate DustEM files can be found in the \(\tt /py\) directory. To use a Python routine "rname", type
  • \(> {\tt from}\;{\tt rname}\; {\tt import}\;*\)
  • \(> {\tt rname()}\).
  • An overall view of the run output can be obtained with "show_dustem.py".

    ➠ Back to top

2. Input files

The main input file for DustEM is "GRAIN.DAT" and must be located in the directory \(\tt /data\). This file contains the main parameters and keywords that define a dust model. Each dust population is defined by its grain type noted tt and its size distribution as set in "GRAIN.DAT". For a given type, DustEM requires the files "/oprop/Q_tt.DAT" and "/hcap/C_tt.DAT" containing the optical properties and heat capacities of grain of type tt, respectively. These files are described in section 2.5. DustEM can handle an arbitrary number of dust populations. In general, a type of grain has a well defined composition (bulk material), structure (porosity, mass density) and shape (spherical or spheroidal). All these properties are taken into account while generating the grain data with DustProp.

2.1 Defining a dust model

Within DustEM, a dust model is defined by the control input file "GRAIN.DAT". We show below a typical example of this file (as in Compiègne et al. 2011).

# DUSTEM: definition of grain populations
#
# run keywords
# G0 scaling factor for radiation field
# grain type, nsize, population keywords, Mdust/MH, rho, amin, amax, alpha/a0 [, at, ac, gamma (ED)] [, au, zeta, eta (CV)]
# cgs units
#
sdist
1.00
PAH0_MC10 10 logn-chrg-zm 7.80E-04 2.24E+00 3.5E-08 1.2E-07 6.4E-08 1.0E-01
PAH1_MC10 10 logn-zm 7.80E-04 2.24E+00 3.5E-08 1.2E-07 6.4E-08 1.0E-01
amCBEx 15 logn 1.65E-04 1.81E+00 6.0E-08 2.0E-06 2.0E-07 3.5E-01
amCBEx 25 plaw-ed 1.45E-03 1.81E+00 4.0E-07 2.0E-04 -2.8E+00 1.5E-05 1.5E-05 2E+00
aSil 25 plaw-ed 7.80E-03 3.50E+00 4.0E-07 2.0E-04 -3.4E+00 2.0E-05 2.0E-05 2E+00

After some documentation lines (starting with #), the run keywords are given (sdist in the above example), then a scaling factor for the intensity of the radiation field (here 1.00). The following lines define the dust populations. From left to right we have: grain type (e.g. PAH0_MC10), number of sizes (nsize), the population keywords (size distribution and physical processes), the dust-to-gas mass ratio (\(Y_i=M_i/M_H\) where \(M_H\) is the mass of protons), the grain material specific mass density (\(\rho\) in g/cm3), the minimum (\(a_{min}\)) and maximum (\(a_{max}\)) sizes in cm and other parameters for the size distribution. As seen in this example for amCBEx, several dust populations can make use of the the same grain type.

Keywords are used to control the DustEM outputs. They can be written in uppercase or lowercase characters. The run keywords control the way DustEM is executed and the output files it produces. If multiple run keywords are used, they must be separated by blanks. The available run keywords are:

  • none: if no run keyword required,
  • quiet: sets the verbose mode off at runtime,
  • res_a: writes output for each grain size,
  • sdist: for each grain size and type, writes the size distribution in the file "/out/SDIST.RES",
  • zdist: for each grain size and type, writes the charge distribution in the file "/out/ZDIST.RES",
  • temp: for each grain size and type, writes the temperature distribution in the file "/out/TEMP.RES",
  • pdr: writes a summary of DustEM output in "/out/EXTINCTION_DUSTEM.DAT". To be set when DustEM is coupled with the Meudon PDR code. This requires the "G_tt.DAT" files to be in the /data directory.

The population keywords (e.g. logn, plaw) define the size distribution and the physical processes included in DustEM for each grain population (e.g. logn, plaw, spin). Population keywords can also be multiple, they are then separated by a dash. The population keywords defining the size distribution are described in section 2.3 and those defining the grain physics in section 2.4.

➠ Back to top

2.2 Template dust models

Several reference dust models are provided. The control input files for the following reference models are included in the current DustEM distribution:

To be used, these files must be renamed to "GRAIN.DAT" within the directory \(\tt /data\).

➠ Back to top

2.3 Size distribution

When the population keyword size is not set, the size distribution is defined from the parameters in "GRAIN.DAT". The size grid has nsize points and is defined over the interval \([a_{min}, a_{max}]\) with a constant logarithmic step. For each grain of mass density \(\rho\), the radius \(a\) is defined as that of an equivalent sphere of same mass \(m = \rho\, 4\pi a^3/3\). The shape of the size distribution \({\rm d}n/{\rm d}a\) (the number of grains with radius in the range \([a, a+da]\)) is defined by the following population keywords and parameters:

  • logn: defines a log-normal size distribution with \({\rm d}n/{\rm d}\log a \sim \exp(-(\log^2(a/a_0)/2\sigma^2)\). The centroid \(a_0\) (in cm) and then width \(\sigma\) are given after \(a_{max}\) in "GRAIN.DAT".
  • plaw: defines a power-law size distribution \({\rm d}n/{\rm d}a \sim a^\alpha\). The index \(\alpha\) is given after \(a_{max}\) in "GRAIN.DAT".
  • ed: applies an exponential decay to a power law size distribution multiplying \({\rm d}n/{\rm d}a\) by the function
    $$ \begin{array}{lll} D(a) = & 1 & {\rm for \;} a\leq a_t \\ & \exp(-\left[(a-a_t)/a_c\right]^\gamma) & {\rm for \;} a>a_t. \end{array} $$ The parameters \(a_t\), \(a_c\) and \(\gamma\) are given in this order after \(\alpha\) in "GRAIN.DAT".
  • cv: applies a curvature term to a power law size distribution multiplying \({\rm d}n/{\rm d}a\) by the function
    $$ C(a) = \left[\, 1 + \left|\zeta\right| (a/a_u)^\eta \,\right]^{{\rm sgn}(\zeta)} $$ If ed is set with cv the parameters \(a_u\), \(\zeta\) and \(\eta; are given in this order after \(\gamma\) in "GRAIN.DAT". Otherwise they are given after \(\alpha\).

The keywords ed and cv thus allow the user to define a size distribution similar to that of Weingartner & Draine (2001).

Alternatively, the size distribution for a given grain population can be read from a file. When the population keyword size is set the size distribution and the size grid are read from the file "/data/SIZE_tt.DAT" which must have the following format (not including the #-lines):

# Size distribution of grain species
#
# Nbr of size bins
# [ a (cm), a^4*dn/da, rho(a) ]
#
5
3.000000E-08 2.264156E-01 2.240000E+00
3.703963E-08 2.674749E-01 2.240000E+00
4.573115E-08 3.159897E-01 2.240000E+00
5.646216E-08 3.733180E-01 2.240000E+00
6.971126E-08 4.410674E-01 2.240000E+00

The "SIZE_tt.DAT" file can be generated with the Python routine "/py/create_vdist.py".

➠ Back to top

2.4 Grain physics

Physical processes affecting the grain properties are taken into account in DustEM thanks to population keywords often associated with data files from where parameters are passed. Below is the list of population keywords for the processes included in DustEM:

  • mix: applies a size dependent weighting factor \f_{mix}\) when computing the contribution of the dust population to the extinction and emission. For each population tt featuring the mix keyword, a file "/data/MIX_tt.DAT" must be given containing a single column with weighting factors for each point of the size grid. The contribution to the emission and extinction of each size bin is then weighted by \f_{mix}\), allowing the user to generate a size dependent mixing of several dust types. For each size, the weighting factors must be consistent with \(\sum_{tt}\,f_{mix}(tt)=1\). In the example "GRAIN.DAT" shown above, neutral (PAH0_MC10) and ionized PAHs (PAH1_MC10) are mixed size by size in the model. The weighing factors (in the "MIX_tt.DAT" files) then correspond to the fraction of neutral and ionized PAHs. In a coarser way, grain types can also be mixed with the mass fraction.
  • chrg: if set, the equilibrium charge distribution of grains \(f(Z)\) per grain size and type will be computed using the Kimura (2016) model for the photoemission yield and the Weingartner & Draine (2001) formalism for the other processes. The plain Weingartner & Draine (2001) model can be used by setting the keyword \(\tt chrgwd\). The computation of the grain charge requires the following data: the gas state in the file "/data/GAS.DAT" (see section 2.6) and the quantities for grain charging in the file "/data/CHRG_tt.DAT":

    # DUSTEM : Constants for charge distribution (CM20 case)
    #
    # Wf(eV) p_uait(3) p_y(3) bulk work function, Uait(V) and Y0 for WD01
    # p_6(2) p_se(3) cex aromatic correction for a-C:H (scale and a/aR), electron sticking, charge exchange
    # p_ke(1) p_le(4) electron mean KE and le(nm)
    # p_ea(2) s_ea(2) nbg EA(nm) and cross-section, size nr for Ebg (if -1 or 1 Ebg no size dependence)
    # Ebg(eV)
    #
    4.75 3.9 0.12 2.0 9.0e-03 3.7e-02 5.0
    1.0 2.0 0.5 1.0 0.0 1.0
    0.0 143.0 -2.0 5.4e-2 0.5
    0.4 0.7 0.5 3.0 -1
    0.0
    ...

    Files with typical values for classical grain types can be found in the \(/data\) directory. In the following we detail some important control parameters. The grain charge distribution is sensitive to the electron sticking coefficient which is taken to be \(s_e=\alpha (1-{\rm e}^{-\beta a/l_e})\) where \(a\) is the grain radius and \(l_e\) is the electron mean free path Kimura (2016) averaged over the distribution of photons above photoemission threshold \(IP_{t}\). The first two values in the array \({\tt p\_se}\) correspond to the Weingartner & Draine (2001) model for \(s_e\) with \(\alpha=0.5\) and \(\beta=1\). If \(\tt p\_se(3)\) is positive, eq. (28) of Weingartner & Draine (2001) (scaled by \(\tt p\_se(3)\)) is used in \(s_e\). For a-C:H grains a correction to the photoemission thresholds (see appendix) needs to be applied and is set by the \(\tt p\_6(2)\) parameter array. The \(\tt cex\) parameter, if positive, allows to take into account charge exchange between C\(^+\) and small grains of the type considered. The rate coefficient of this process is taken from Wolfire et al. (2008) and is scaled with \(\tt cex\). The \(\tt p\_ke\) parameter (\(0 < {\tt p\_ke} \leq 1\)) sets the fraction of the maximum available energy \(h\nu-IP_{t}\) that goes into kinetic energy of the photoelectron released. If \({\tt p\_ke} \leq 0\) the Weingartner & Draine (2001) model is used. The size dependence of the band gap \({\tt Ebg}\) is controlled by the integer \({\tt nbg}\): if \({\tt nbg}\) = 0 then \({\tt Ebg} = 0\), if \({\tt nbg} = 1\) the value given is used for \({\tt Ebg}\) and if \({\tt nbg} = -1\) the band gap law for a-C:H grains will be used, namely \(E_g({\rm eV})={\rm MAX}\left[0.1,\;0.2\left({5\,{\rm nm}\over a}-1\right)\right]\) (see appendix). Other size dependence laws of the bandgap can be used with \({\tt nbg} > 0\) the number of size bins. The \({\tt Ebg}\) values are then expected as a single column with \({\tt nbg}\) lines corresponding to the number of size bins of the grain type (note that, having set the run keyword \({\tt sdist}\), you can retrieve the size bin values in the file "SDIST.RES"). Finally, output on the charge of grains is obtained by adding the run keyword \({\tt zdist}\) (or \({\tt zdistf}\) to get the charge distribution \(f(Z)\) of all grains).
  • zm: allows to mix populations of neutral and charged grains with fractions given by the charge distribution \(f(Z)\). This must be set for two consecutive grain populations in "GRAIN.DAT" (see section 2.1). For the first population (neutral grains) the keyword chrg-zm must be set and a weighting factor \(f_{mix}=f(0)\) is applied for each size to the grain emission and extinction. For the second population (charged grains), chrg-zm is set again and the weighting factor for each size is then \(f_{mix}=\sum_{|Z|>0} f(Z)=1-f(0)\). This mixing can be applied to any grain population for which \(f(Z)\) can be computed and does not require any "MIX_tt.DAT" file. It is particularly useful in the computation of the emission and extinction of neutral and charged PAHs. For the gas heating and cooling, it is assumed that the neutral and charged species have the same rates which are therefore added up for the first species with chrg-zm. In case this assumption is not verified, the rates for neutral and charged species must be obtained from distinct grain types each featuring the \(\tt chrg\) keyword (i.e. without using \(zm\)).
  • spin: turns on the emission of spinning PAHs following Draine & Lazarian (1998) for the gas-grain interactions and using the emission as derived in DustEM for the radiative IR excitation rates. The permanent electric dipole moment is assumed to be \(\mu=m\,N_{1/2}\) where \(m\) is given in the file "/data/SPIN_tt.DAT" with the mean atomic (or molecular) weight of the grain (can be size dependent) to estimate \(N\) the number of atoms (or molecules) in the grain. The gas parameters are given in "/data/GAS.DAT" (see section 2.6). This option triggers the computation of the charge distribution of each grain.
  • pol: turns on the computation of the polarized emission and extinction. This requires the files "Q1_tt.DAT" and "Q2_tt.DAT" to be in the \(\tt /oprop\) directory. These tabulated cross sections for polarization must be averaged over the grain dynamics around magnetic field lines. The fraction of aligned grains or alignment function \(f(a)\) is size dependent. We use a parametric step function $$ f(a) = {1\over 2} f_{max} \left( 1 + \tanh[ \ln(a/a_{alig})/p_{stiff}] \right) $$ where \(f_{max}\) is the maximum fraction of aligned grains, \(a_{alig}\) is the size corresponding to \(f=0.5\). and \(p_{stiff}\) controls the size range for the transition between \(f\) = 0 and 1. Parameters for the polarization model must be given in the file "/data/ALIGN.DAT", please refer to Guillet et al. (2018) for a detailed description of parameters and files.

    # DUSTEM: definition of grain alignment
    # for each grain TYPE make sure you have the following files in ../oprop
    # Q1 TYPE.DAT and Q2 TYPE.DAT
    #
    # run keywords
    # lin (for linear), univ (for universal alignment law)
    # degree of anisotropy of the radiation : keep it to 0.
    # par (parametric alignment law) followed by parameters a alig, p stiff and f_max
    #
    lin univ
    0.00E+00
    par 8.0E-02 2.7E-01 6.7E-01

    If the univ keyword is set, the same alignment function f(a) is used for all grain types and only the first parameters line (par ... ) in ALIGN.DAT will be used to define f(a). If the univ keyword is not set, then parameters lines must come in same number and order as that of the grain types for which the pol keyword was added in the file GRAIN.DAT.

  • beta: applies a correction to Qabs to introduce a temperature dependence, namely $$ Q_{abs}(a,\nu)=Q_0(a,\nu)\,\left(\nu/\nu_t\right)^{\delta (T)\,H(\nu_t/\nu)} $$ where Q0 is the grain emissivity without the temperature dependence, \(\nu_t=c/\lambda_t\) is the frequency threshold and \(\delta (T)=\beta (T)- \beta_0\) is the index correction with respect to β0 some reference value. The threshold function in frequency is \(H(x)=0.5(1+\tanh(4x/s))\) where \(x=\log(\lambda/\lambda_t\) and s controls the stiffness of the transition from 0 to 1 (the width \(\Delta x\) is proportional to s): for instance for \(s=1\), \(H(x)\) rises from 0 to 1 for \(x\) between 0.5 to 1.5 and \(\Delta x=1\). Parameters for this correction are in "/data/BETA_tt.DAT":

    # DUSTEM: parameters for BETA(T) behaviour of Qabs
    # first DBETA correction to standard index beta0
    # then lambda threshold to apply BETA(T)
    #
    # beta0 a g bmax
    # lthresh(microns) lstiff (dlambda=lthresh/lstiff)
    # nbetav (nr of beta values if 0 uses formula)
    # tbeta betav (if nbetav >= 0)
    #
    2.11E+00 1.15E+01 -6.60E-01 5.00E+00
    3.00E+01 1.00E+00
    13
    1.00E+00 1.55E+00
    ...

    where lthresh and lstiff stand for the above λt and s. If nbetav=0 then \(\beta (T)=aT^g\) (Désert et al. 2008) and we take β=MIN(β,bmax). If nbetav> 0 the β(T)-values must be given below nbetav in two columns T and β.
  • dtls: adds the low temperature opacities of amorphous grains as described in Mény et al. (2007). Uses the file "/data/DTLS_tt.DAT":

    # DUSTEM: constants for the DCD and TLS effects
    # a_dtls lc(nm) c_delta
    # vt(cm/s) P*mub$^2$ gamma_e(eV)
    # omega_m(s-1) tauO(s-1) Vo(K) Vmin(K) Vm(K)
    # ldtresh(microns)
    #
    5.00e-02 2.00e-02 1.30e+01 1.87e+03
    3.00e+05 1.40e-03 1.00e+00
    2.69e+12 1.00e-13 410.0 50.0 550.0
    1.00e+02
    where all parameters are from Mény et al. (2007), with adtls the relative weight of the DCD to TLS effects and ldtresh is the threshold above which the correction is applied (scaled to the current Q value at wavelength ldtresh). The absorption efficiency per grain type and size is thus $$ \begin{array}{lll} Q_{abs}(a,\lambda) = & Q_{DCD} + {\tt a_{dtls}}\,Q_{TLS} & {\rm for \;} \lambda\geq {\tt ldthresh} \\ & Q_0 & {\rm for \;} \lambda < {\tt ldthresh}. \end{array} $$ where Q0 is the absorption efficiency without the DCD-TLS effects read in the "Q_tt.DAT" file (see section 2.5). The DCD-TLS absorption efficiencies are defined such that for \(\lambda = {\tt ldthresh}\) one has \(Q_0 = Q_{DCD} + {\tt a_{dtls}} Q_{TLS}\).
  • ➠ Back to top

    2.5 Grain data

    For a given type of grain tt, the dust properties are provided as a function of the size a in the data files, namely: "Q_tt.DAT" files for the absorption and scattering efficiencies, Q(a,λ), and "C_tt.DAT" files for the heat capacity per unit volume, C(a,T). The Q and C values are provided over the broadest possible size range (usually 0.3 to 104 nm). A grain type can be used once these Q and C-files exist in the DustEM data file set, in the optical properties \(/oprop\) and heat capacity \(/hcap\) directories, respectively. The sizes requested in "GRAIN.DAT" or "SIZE_tt.DAT" have to fit into those of the above files (i.e. extrapolation is not allowed). When DustEM is coupled to radiative transfer codes (e.g. PDR, CRT, SKIRT), the g-factors (anisotropy factors for scattering defined as \(<\cos \theta>\), where \(\theta\) is the scattering angle) are required at all wavelengths and for all sizes: they are given in "/oprop/G_tt.DAT". In the case of spherical grains, the Q and G files are generated with the Mie theory (Bohren & Huffman 1983) while for spheroidal grains the T-matrix method is used (Mischenko & Travis 1998) with DustProp. We briefly outline below the format of these files.

    The Q-files are as follows:

    # QABS and QSCA to be used by DUSTEM
    #
    # nsize (number of sampled sizes in this file)
    # sizes (microns)
    50
    3.0000E-04 3.7104E-04 4.5891E-04 5.6759E-04 7.0200E-04 8.6824E-04 1.0738E-03 1.3281E-03 ...
    #
    #### QABS ####
    1.3461E-03 1.6642E-03 2.0571E-03 2.5419E-03 3.1398E-03 3.8762E-03 4.7819E-03 5.8952E-03 ...
    ...
    #### QSCA ####
    5.0179E-07 1.1650E-06 2.6935E-06 6.1884E-06 1.4083E-05 3.1594E-05 6.9380E-05 1.4766E-04 ...

    The first uncommented line gives the number of sizes present in the file. The line next gives the sizes (in microns). Then for each grain size (one column per size), the absorption Qabs and scattering Qsca efficiencies are given. The G-files are as per the Q-files with the QABS field replaced by g(a,λ) and no QSCA field. The Q and G-files are all generated with the common wavelength grid given in "/oprop/LAMBDA.DAT".

    The C-files look like this:

    # DUSTEM heat capacity
    #
    # nr of sizes ns
    # sizes (microns)
    # nr of T-values
    # log T(K) log C_1...log C_ns (erg/K/cm3)
    50
    3.0000E-04 3.7104E-04 4.5891E-04 5.6759E-04 7.0200E-04 8.6824E-04 1.0738E-03 1.3281E-03 ...
    30
    -1.000E+00 4.8983E-01 4.8983E-01 4.8983E-01 4.8983E-01 4.8983E-01 4.8983E-01 4.8983E-01 4.8983E-01 ...

    The first uncommented line gives the number of sizes present in the file and the second line gives the sizes in microns. The third line gives the number of points in the temperature grid. The following gives the log10 of temperatures (in Kelvin, column 1) and the log10 of heat capacities (one column for each size) in units of erg/K/cm3. The same temperature grid (30 points from 0.1 to 5000 K with a constant log step) is used for all C-files. This temperature grid is given in "/hcap/TEMP.DAT".

    The grain data can be batch generated with the DustProp package (not included in the distribution, for internal use only at the moment). We list below the grain types for which Q and C files are available in the current DustEM distribution:

    • PAH0_MC10, PAH1_MC10: PAHs as defined in Compiègne et al. (2011) for neutral and singly charged PAHs, respectively,
    • PAH0_DL01, PAH1_DL01: PAHs as defined in Draine & Li (2001) for neutral and singly charged PAHs, respectively. It takes into account the transition to graphite optical properties for a > 5 nm,
    • PAH0_DL07, PAH1_DL07: PAHs as defined in Draine & Li (2007) for neutral and singly charged PAHs, respectively. It takes into account the transition to graphite optical properties for a > 5 nm,
    • PAH1_MC08: PAHs described in Compiègne et al. (2008),
    • Gra: graphite grains. The Q-values have been generated with the Mie method for spherical particles. The refractive index and heat capacity are as described in Draine & Li (2001) and Li & Draine (2001).
    • amCBEx: amorphous carbon grains. The Q-values have been generated with the Mie method. The refractive index is from Zubko et al. (1996) and the Qabs have been extrapolated above 0.9 mm (see Compiègne et al. 2011). The heat capacity is the same as that of graphite.
    • aSilx: astronomical silicates. The Q-values have been generated with the Mie method. The refractive index and heat capacity are as in Draine & Li (2007), originally from Draine & Lee (1984). The Qabs have been extrapolated above 0.9 mm (see Compiègne et al. 2011),
    • BG_MC08: big grain as defined in Désert et al. (1990); modified by Compiègne et al. (2008).
    • CM20: carbonaceous grains as defined in Jones et al. (2013). Grains up to 20 nm consist purely of aromatic-rich H-poor amorphous carbon, a-C, whereas bigger grains have a core/mantle structure, where the core is made of aliphatic-rich H-rich carbon, a-C:H, and the mantle of a-C with a thickness of 20 nm,
    • aOlM5 + aPyM5: silicates as defined in Koehler et al. (2014). These core/mantle grains consist of amorphous silicate cores (olivine and pyroxene-normative composition, Mg-rich) and 5 nm mantles of a-C.
    • amCBE_0.3333x: spheroidal amorphous carbon grains (see above) including the data for polarized emission and extinction Guillet et al. (2018),
    • aSil_0.3333_p20B, aSil2001_0.3333_p20B and aSil2001BE6pctG_0.4x: spheroidal astronomical silicates (see above) including the data for polarized emission and extinction Guillet et al. (2018).
    Data files for the Désert et al. (1990) model are also available corresponding to the grain types PAH0_DBP90, PAH1_DBP90, VSG_DBP90 and BG_DBP90.

    Grains with size-dependent mass density \(\rho(a)\) can also be handled. The mass density values need to be calculated on the size grid of the "Q_tt.DAT" file and must be given in that file right after the size array. To use these \(\rho(a)\) values for a given grain type \(\tt tt\) you then need to turn negative the corresponding mass density in the GRAIN.DAT file. If the keyword \(\tt size\) is used for the grain type \(\tt tt\), the \(\rho(a)\) values from the Q-file will supersede those given in the "SIZE_tt.DAT".

    ➠ Back to top

    2.6 Astrophysical data

    The radiation field used by DustEM is in "/data/ISRF.DAT". This file contains 2 columns, the wavelength (microns) and the flux \(4\pi I_\nu\) in erg/s/cm2/Hz. The wavelength range must be included in that of the Q-file. It can be generated with the "create_rfield.py" routine in the \(\tt /py\) directory. The current distribution includes "ISRF_MATHIS.DAT", the standard radiation field of Mathis et al. (1983) and "ISRF_CMB.DAT", where the CMB contribution has been added to the Mathis field.

    Parameters for the gas are given in "/data/GAS.DAT" as:

    # DUSTEM : Gas quantities CNM from DL98b and Weingartner & Draine 2001
    #
    # Gas temperature (K),hydrogen density,H2 density,CR rate,G0 (if>0 supersedes GRAIN.DAT),nr of charge type (e-, H+, C+, ...)
    # ion density (cm-3), mass (amu), charge, polarizability (A^3)
    # >>>>> 1st line is electron <<<<<<<
    # >>>>> 2nd line is proton <<<<<<<
    # >>>>> 3d line is carbon <<<<<<<
    #
    1.00E+02 3.00E+01 0.00E+00 5.00E-17 1.00E+00 3
    3.00E-02 5.4858E-04 -1.00 0.00
    3.60E-02 1.00794 1.00 0.67
    9.00E-03 12.0107 1.00 1.54

    If the gas temperature value \(\tt tval\) given in "GAS.DAT" is negative, this value (turned positive) will be used to estimate the gas temperature from a thermal pressure to \(G_0\) relationship Joblin et al. (2018), namely \(T_{gas}=|{\tt tval}|\times G_0/n_H\). Similarly, if the given H\(_2\) density is negative, the density computed from the formation-dissociation equilibirum is used (with an assumed dissociation rate of \(5\,10^{-11}\) s\(^{-1}\)). If negative the \(G_0\) value in "GAS.DAT" will not be used. Otherwise, if positive, this G0 value will supersede the value given in "GRAIN.DAT".

    The abundances of the main charged species (e-, \(p\) and C\(^+\)) are important for the grain charge and the gas photoelectric heating. Values for the electron and C\(^+\) density must be given in the "GAS.DAT" file. The given proton density value will be superseded by \(n_p = n_e\) - \(n_m\) with \(n_e\) the electron density and \(n_m\) the density of cations heavier than the proton (here only C\(^+\)). Positive values for \(n_e\) and \(n_m\) as given in "/data/GAS.DAT" will be used directly.

    If \(n_e\) is negative, the value at ionization equilibrium (cosmic ray ionization balanced by proton recombination) is used complemented by the approximation of Williams et al. (1998) in dense, cold gas and $$ n_e = n_H,\text{MAX}(x,x_{dc}) $$ with \(x_{dc}=10^3\,\left(\zeta\over n(\text{H}_2)\right)\) and \(x=n_e/n_H = {1\over 2}\,\left(x_m-x_{cr}+\sqrt{(x_m+x_{cr})^2+4x_{cr}}\right)\) where \(x_{m}\) is the abundance of of ionized metals (C\(^+\), S\(^+\), Si\(^+\),...), \(x_{cr}=\zeta/(\alpha n_H)\) with \(\zeta\) the cosmic ray ionization rate per proton and \(\alpha\) the proton recombination rate (in cm\(^3\)/s) as taken from Roellig et al. (2006)

    If the given \(n(\)C\(^+)\) is negative, the electron density \(n_e\) is first computed from the coupled ionization equilibria of protons and C\(^+\) possibly involving charge exchange reactions with grains. This derivation couples two equilibrium equations: the cosmic ray/proton recombination and the photoionization/neutralization of C\(^+)\), the latter process involving reactions with H\(_2\) and charge exchange with PAHs (see Ysard et al. 2011). The charge exchange with grains of type \(\tt tt\) is controlled from the "CHRG_tt.DAT" file by setting the charge exchange parameter \(\tt cex\) to a strictly positive value (see sec. 2.4). The C\(^+\) density is then calculated assuming \(n_m = n\)(C\(^+\)). Finally the proton density is derived as described above. These results also depend on the total carbon gas abundance which is taken to be the absolute value of the (negative) number given for \(n\)(C\(^+\)) in the "GAS.DAT" file.

    Several template "GAS.DAT" files can be found in the \(\tt /data\) directory:

    ➠ Back to top

    3. Output files

    The output files produced by DustEM are stored in the \(/out\) directory. They give the dust emission and extinction. Each file begins with documentation lines (first character is #). Any DustEM run will produce the 2 following files:

    • SED.RES: contains the emission per proton for each grain population \(4\pi\nu I_\nu/N_H\) in erg/s/H as a function of wavelength (\(N_H\) is the proton column density). The last column is the total of all previous columns.
    • EXT.RES: contains the extinction opacity for a column density \(N_H = 10^{21}\) H cm\(^{-2}\) for each grain population. The last column is the sum of all previous columns, i.e. the dust optical depth:
      $$ \tau(\lambda) = n_H\,\sum_{\rm tt}\,(\sigma^{\tt tt}_{abs}(\lambda) + \sigma^{\tt tt}_{sca}(\lambda)). $$ The magnitude of extinction thus being given by \(A_\lambda = 1.086\times \tau(\lambda) \).

    If the run keywords temp, zdist or sdist are used, the files TEMP.RES, ZDIST.RES or SDIST.RES will be written and will contain, respectively, the temperature, charge or the size distribution of the dust populations. Full temperature and charge distribution per grain type and size are obtained if \(\tt tempf\) or \(\tt zdistf\) are set. All these files have a documented header.

    Using the population keywords pol or spin will generate the following files:

    • SED_POL.RES: contains the polarized emission of grains with the same format as SED.RES.
    • SPIN.RES: contains the emission of spinning grains with the same format as SED.RES.

    Using the run keyword res_a will generate the files SED_A.RES, SED_POL_A.RES and SPIN_A.RES, which contain the same information, namely, the emission per grain in size bin \([a,a+da]\) or \(4\pi\nu I_\nu/N_d(a)\) in erg/s/grain (\(N_d(a)\) is the column density of grains of size a).

    ➠ Back to top

    4. DustEM Python

    Below is a list of the Python routines in the \(\tt /py\) directory with the use for each. To use a routine \({\tt rname}\), documentation is obtained by typing
    \({\tt >>> from\; rname\; import\; *}\)
    \({\tt >>> rname()}\).

    • get_band_flux (in band_flux): returns structure (dictionary) with band or spectral flux (MJy/sr) in instrument filter,
    • create_rfield: generates the radiation field in erg/cm2/s/Hz for DustEM (ISRF.DAT),
    • create_vdist: generates a normalized dust size distribution for DustEM (SIZE_tt.DAT),
    • dm4cldy: generates a grain opacity file for Cloudy from DustEM data,
    • rebin_oprop: rebin DustEM optical properties of a given grain type from LAMBDA1.DAT to LAMBDA2.DAT,
    • show_dustem: shows (and returns) the results of DUSTEM overlaid on diffuse ISM dust observations,
    • show_tdist: shows temperature distribution of grains computed by DustEM,
    • show_zdist: shows charge distribution of grains computed by DustEM.