Userguide

The userguide can be read below or 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 108 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 dP/dT is derived with the formalism of Désert et al. (1986) assuming that the grain cooling is continuous. To correctly describe the long wavelength emission of dust, the initial algorithm has been revised to cover the low temperature part of dP/dT. The fortran files are in the /src directory, data files in the /oprop, /hcap, and /data directories and output files are in the /out directory. All data files have headers (where each line begins with #) that document their content. IDL routines useful to generate or visualize DustEM files can be found in the /pro directory.

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 Size distribution
  2.4 Grain physics
  2.5 Grain data
  2.6 Astrophysical data
3. Output files

 

1. Getting started

After downloading the tar ball from here:

  • Unpack the archive: "tar xvzf dustem.tar.gz". This will generate the /dustem directory which should contain the following subdirectories: /data, /hcap, /oprop, /out, /pro, and /src.
  • In /src, edit the file "DM_constants.f90" to define the string 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 dir_PDR as the absolute path to /dustem/out.
  • Edit "Makefile" and set up the standard commands corresponding to your fortran compiler (gfortran, g95 and ifc have been tested - ifc is about twice as fast as the others), 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 That's it!. The executable is called dustem.
  • Running: ./src/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 check your run by using "show_dustem.pro" in /pro.
  • ➠ Back to top

2. Input files

The main input file for DustEM is "GRAIN.DAT" and must be located in the directory /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 (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 (Yi = Mi/MH where MH is the mass of protons), the grain material specific mass density (ρ in g/cm3), the minimum (amin) and maximum (amax) 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 /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 [amin, amax] with a constant logarithmic step. For each grain of mass density ρ, the radius a is defined as that of an equivalent sphere of same mass m = ρ 4π a3/3. The shape of the size distribution dn/da (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 dn/dlog a ∝ exp(-(log(a/a0)/σ)2. The centroid a0 (in cm) and then width σ are given after amax in "GRAIN.DAT".
  • plaw: defines a power-law size distribution dn/da ∝ aα. The index α is given after amax in "GRAIN.DAT".
  • ed: applies an exponential decay to a power law size distribution multiplying dn/da by the function
    D(a) = 1 for a ≤ at
            exp(-[(a-at)/ac]γ) for a>a_t.
    The parameters at, ac and γ are given in this order after α in "GRAIN.DAT".
  • cv: applies a curvature term to a power law size distribution multiplying dn/da by the function
    C(a) = [ 1 + |ζ| (a/au)η]sgn (ζ)
    If ed is set with cv the parameters au, ζ and η are given in this order after γ in "GRAIN.DAT". Otherwise they are given after α.

The keywords ed and cv 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
#
# shape (t-matrix code, 0:spheres, -1:spheroids) eps (b/a, eps>1:oblate, eps<1: prolate)
# Nbr of bulk materials
# Name of bulks
Bulk densities in g/cm3
# Mass fractions for each bulk
# Nbr of size bins
# [ a (cm), dloga, a^4*dn/da, rho_eff, fv ]
# fv: volume fraction of bulks, rho_eff: volume mean density
#
0 1.00E+00
1
Gra
2.2500E+00
1.0000E+00
46
8.0876E-08 2.0939E-01 7.1858E-02 2.2500E+00 1.0000E+00
....

The "SIZE_tt.DAT" file can be generated with the IDL routine "/pro/create_vdist.pro".

➠ 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 fmix 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 fmix, allowing the user to generate a size dependent mixing of several dust types. For each size, the weighting factors must be consistent with ttfmix(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, species can also be mixed with the mass fraction.
  • beta: applies a correction to Qabs to introduce a temperature dependence, namely:
    Qabs(a, ν) = Q0(a, ν) (ν/νt)(δ(T) H(νt/ν)
    where Q0 is the grain emissivity without the temperature dependence, νt=c/λt is the frequency threshold and δ(T)=β(T)-β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(λ/λt) and s controls the stiffness of the transition from 0 to 1 (the width Δ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 Δx=1. Parameters for this correction should be 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 β(T)=aTg (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 β.
  • chrg: if set, the equilibrium charge distribution of grains per grain size and type will be computed using the Weingartner & Draine (2001) formalism with the van Hoof et al. (2004) modification for the photodetachment (see also Weingartner et al. 2006). This option 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
    #
    # Wf(eV) p_ea(2)(nm) for IP and EA
    # s_ea(2) for EA cross-section: f and dE(eV)
    # le(A) p_y(3) for PE yield
    # p_uait(3) for Uait
    # m1mass molecular mass in grain (g)
    # nr of Ebg values (if 1 Ebg constant equal to 1st value)
    # Ebg(eV)
    #
    4.4 4.0 7.0
    0.5 3.0
    10. 9.0e-03 3.7e-02 5.0
    3.9 0.12 2.0
    12.01
    1
    0.0
    ...
    The bandgap Ebg can be made size-dependent: in this case the values must be provided in one column with as many lines as the number of sizes for the given grain type. When Ebg has a single value (one line), this value is applied for all sizes. The grain charge distribution is sensitive to the electron sticking coefficient se which by default is that of the Weingartner & Draine (2001) model. The case of the Bakes & Tielens (1994) model (s_e=1) can also be used by using the composed keyword chrg-bt. The Bakes & Tielens (1994) model has also been widely used in astrophysical modeling and this prescription for se has noticeable consequences (Verstraete 2011).
    The charge distribution of small grains (a ≤ 10 nm) is computed on a charge grid around Zeq, the charge value corresponding to zero current, i.e., Jpe+Jion=Je (see Weingartner & Draine (2001) for notations). The number of charge points is MIN(Zmax-Zmin+1,nzmax) with nzmax=30 by default. In the case of larger grains the charge distribution is strongly peaked around Zeq and the charge grid has nzmin around Zeq (nzmin=20 by default).
  • zm: allows to mix populations of neutral and charged grains with fractions given by the charge distribution fZ(Z). This must be set for the two following grain populations in "GRAIN.DAT" (see section 2.1). For the first population (neutral grains) zm must be set with chrg and a weighting factor fmix=fZ(0) is applied for each size to the grain emission and extinction. For the second population (charged grains), zm is set alone and the weighting factor for each size is then fmix=∑|Z|>0 fZ(Z)=1-fZ(0). This mixing can be applied to any grain population 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.
  • 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_dtl 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:
    Qabs (λ ≥ ldtresh)= QDCD + adtls QTLS) and
    Qabs (λ < ldtresh)= Q0
    where Q0 is the absorption efficiency without the DCD-TLS effects read in the "Q_tt.DAT" file (see section 2.5).
  • 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 /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/2 fmax ( 1 + tanh(ln(a/aalig)/pstiff) )
    where fmax gives the maximum fraction of aligned grains, a/aalig is the size corresponding to f=0.5 and pstiff 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 (2017) 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.
  • 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 μ=ma3/2 where m is given in the file "/data/SPIN_tt.DAT". The gas parameters are given in "/data/GAS.DAT" (see section 2.6). This option requires the charge distribution of each grain to be computed.

➠ 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 factorq for scattering defined as cosθ where θ 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:

30
# 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 ...
-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 IDL package (not included in the /pro directory, 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_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 (Gra_Dr03) 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).
  • Gra_Dr03: graphite grains, same as Gra but with the changes from Draine 2003 at short wavelengths.
  • 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.
  • aSil: 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),
  • 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.
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.

In the case of more complex grain types with size-dependent properties (composition, porosity) the Q, C and G files are generated with the information in the "SIZE_tt.DAT" file.

➠ 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πIν 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.pro" routine in the /pro 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
#
# Gas temperature (K), hydrogen density, H2 density
# number of charge type (electrons, H+, C+, ...)
# ion density (cm-3), mass (amu), charge, polarizability (A^3)
# >>>>> 1st line is electron <<<<<<<
#
1.00E+02 3.00E+01 0.00E+00 4.50E-02
3
4.50E-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

➠ 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πνIν/NH in erg/s/H as a function of wavelength (NH 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 NH = 1021 H cm-2 for each grain population. The last column is the sum of all previous columns, i.e. the dust optical depth:
    τ(λ) = NHttabsttscatt)
    Aλ = 1.086 × τ(λ)

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. Each of these files has a documented header.

Use of 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πν Iν/Nd(a) in erg/s/grain (Nd(a) is the column density of grains of size a).

➠ Back to top