myXCLASS

creating single, synthetic spectra

The myXCLASS function is the core of XCLASS and models a spectrum by solving the radiative transfer equation for an isothermal object in one dimension, called detection equation [Stahler_Palla_2005]. The myXCLASS function is designed to describe line-rich sources which are often dense, where LTE is a reasonable approximation. Details of the call of the myXCLASS function are described in Sect. “myXCLASS”.

Example call of the myXCLASS function:

>>> from xclass import task_myXCLASS
>>> import os

# get path of current directory
>>> LocalPath = os.getcwd() + "/"

## define path and name of molfit file
>>> MolfitsFileName = LocalPath + "files/12C.molfit"

## define path and name of obs. xml file
>>> ObsXMLFileName = LocalPath + "files/SMF_obs.xml"

# call myXCLASS function
>>> spectra, log, te, IntOpt, JobDir = task_myXCLASS.myXCLASSCore(MolfitsFileName = MolfitsFileName, \
                                                                  ObsXMLFileName = ObsXMLFileName)

What is myXCLASS?

The myXCLASS function is able to model a spectrum with an arbitrary number of molecules where the contribution of each molecule is described by an arbitrary number of components. These could be identified with clumps, hot dense cores etc. The myXCLASS function offers the possibility to describe the position and shape of each component in a simple 1-d, see Sect. “Simple 1-d description”, or in an advanced 3-d description, see Sect. “Stacking” and Sect. “Sub-beam description”. Additionally, myXCLASS is able to take the local and non-local overlap between different lines into account, see Sect. “Local-overlap of neighboring lines”.

In general, the myXCLASS function assumes, that all components are aligned in a user-defined stack along the line-of-sight, where the background and the observer are located at the ends, respectively. The solution of the radiative transfer equation for components which are located in a layer (i=1) closest to the background is 3,

(1)\[\begin{split}T_{\rm mb}(\nu) = \sum_{m,c \in i} &\Bigg[\eta \left(\theta^{m,c}\right) \left[S^{m,c}(\nu) \left(1 - e^{-\tau_{\rm total}^{m,c}(\nu)}\right) + I_{\rm bg} (\nu) \left(e^{-\tau_{\rm total}^{m,c}(\nu)} - 1\right) \right] \Bigg] \\ &+ \left(I_{\rm bg}(\nu) - J_\mathrm{CMB} \right),\end{split}\]

where the sums go over the indices \(m\) for molecule, and \(c\) for component, respectively. In the following we will briefly describe each term in Eq. (1).

The beam filling (dilution) factor \(\eta(\theta^{m,c})\) of molecule \(m\) and component \(c\) in Eq. (1) for a source with a Gaussian brightness profile, see below, and a Gaussian beam is given by 4

(2)\[\eta(\theta^{m,c}) = \frac{(\theta^{m,c})^2}{\theta_t(\nu)^2 + (\theta^{m,c})^2},\]

where \(\theta^{m,c}\) and \(\theta_t\) represents the source and telescope beam full width half maximum (FWHM) sizes, respectively. The sources beam FWHM sizes \(\theta^{m,c}\) (in arcsec) for the different components are defined by the user in the molfit file, described in Sect. “The molfit file”. Additionally, the myXCLASS program assumes for single dish observations (indicated by Inter_Flag = False), that the telescope beam FWHM size is related to the diameter of the telescope by the diffraction limit

(3)\[\theta_t(\nu) = \left(1.22 \cdot \frac{\lambda}{D}\right) \cdot \xi = \left(1.22 \cdot \frac{c_{\rm light}}{\nu \, D}\right) \cdot \xi,\]

where \(D\) describes the diameter of the telescope, \(c_{\rm light}\) the speed of light, and \(\xi = 3600 \cdot 180 /\pi\) a conversion factor to get the telescope beam FWHM size in arcsec. For interferometric observations (indicated by Inter_Flag = True, the user has to define the interferometric beam FWHM size directly using parameter TelescopeSize. In contrast to single dish observations the myXCLASS program assumes a constant interferometric beam FWHM size for the whole frequency range.

In general, the brightness temperature of radiation temperature \(J(T, \nu)\) is defined as

(4)\[J(T, \nu) = \frac{h \, \nu}{k_B}\frac{1}{e^{h \, \nu / k \, T} - 1}.\]

The expression \(J_\mathrm{CMB}\) describes the radiation temperature Eq. (4) of the cosmic background \(T_{\rm cbg}\) = 2.7 K, i.e. \(J_\mathrm{CMB} \equiv J(T_{\rm cbg}, \nu)\).

In Eq. (1), the expression \(S^{m,c}(\nu)\) represents the source function and is according to Kirchhoff’s law of thermal radiation given by

(5)\[\begin{split}S^{m,c}(\nu) &= \frac{\epsilon_l^{m,c}(\nu) + \epsilon_{\rm dust}^{m,c}(\nu)}{\kappa_l^{m,c}(\nu) + \kappa_{\rm dust}^{m,c}(\nu) } \\ &= \frac{\kappa_l^{m,c}(\nu) \, J(T_\mathrm{ex}^{m,c}, \nu) + \kappa_{\rm dust}^{m,c}(\nu) \, J(T_{\rm dust}^{m,c}, \nu) }{\kappa_l^{m,c}(\nu) + \kappa_{\rm dust}^{m,c}(\nu)} \\ &= \left(1 - \delta_{\gamma, 0}\right) \times \Bigg[ \frac{\tau_l^{m,c}(\nu) \, J(T_\mathrm{ex}^{m,c}, \nu) + \tau_{\rm dust}^{m,c}(\nu) \, J(T_{\rm dust}^{m,c}, \nu) }{\tau_l^{m,c}(\nu) + \tau_{\rm dust}^{m,c}(\nu) }\Bigg] + \delta_{\gamma, 0} \, J(T_\mathrm{ex}^{m,c}, \nu),\end{split}\]

where \(\epsilon_{i}^{m,c}(\nu)\) and \(\kappa_{i}^{m,c}(\nu)\) are the emission and absorption coefficients for line and dust, respectively. Additionally, the optical depth is given by \(\tau^{m,c}_i(\nu) = \int \kappa^{m,c}_i(\nu) \, ds\), which assumes that molecules and continuum contributions are well mixed. In older versions, the background temperature could only be defined as the measured continuum offset, which corresponds to the beam-averaged continuum brightness temperature. At the same time, the dust, as agent of line attenuation, was described by column density and opacity. This is practical, because the observable \(T_{\rm bg}\) is used, but does not constitute a self-consistent and fully physical description. Therefore, we now use optionally either a physical (\(\gamma \equiv 1\), defined by the input parameter setting t_back_flag = False) or phenomenological (\(\gamma \equiv 0\), defined by the input parameter setting t_back_flag = True) description of the background indicated by the Kronecker delta \(\delta_{\gamma, 0}\), i.e. \(S^{m,c}(\nu) \equiv J(T_\mathrm{ex}^{m,c}, \nu)\) for \(\gamma \equiv 0\). Note, if \(\gamma \equiv 0\), the definition of the dust temperature \(T_d^{m,c} (\nu)\), Eq. (20), is superfluous.

The total optical depth \(\tau_{\rm total}^{m,c}(\nu)\) of each molecule \(m\) and component \(c\) is defined as the sum of the optical depths \(\tau_l^{m,c}(\nu)\) of all lines of each molecule \(m\) and component \(c\) plus the dust optical depth \(\tau_{\rm dust}^{m,c}(\nu)\) (see Sect. “Dust contribution”), i.e.

(6)\[\tau_{\rm total}^{m,c}(\nu) = \tau_l^{m,c}(\nu) + \tau_{\rm dust}^{m,c}(\nu).\]

The optical depth \(\tau_t^{m,c}(\nu)\) for a transition \(t\) of molecule \(m\) and component \(c\) is described as 5

(7)\[\tau_t^{m,c}(\nu) = \Bigg[\frac{c_{\rm light}^2}{8 \pi \nu^2} \, A_{ul}^t \, N_{\rm tot}^{m,c} \, \frac{g_u^t \, e^{-E_l^t/k_B \, T_{\rm ex}^{m,c}}}{Q \left(m, T_{\rm ex}^{m,c} \right)} \left(1 - e^{-h \, \nu^t /k_B \, T_{\rm ex}^{m,c}} \right) \cdot \phi^{m,c,t}(\nu)\Bigg],\]

where \(\phi^{m,c,t}(\nu)\) indicates the Gaussian line profile function, described in Sect. “Gaussian line profile”. The Einstein \(A_{ul}\) coefficient 6, the energy of the lower state \(E_l\), the upper state degeneracy \(g_u\), and the partition function 7 \(Q \left(m, T_{\rm ex}^{m,c} \right)\) of molecule \(m\) are taken from the embedded SQLite3 database. In addition, the values of the excitation temperatures \(T_{\rm ex}^{m,c}\) and the column densities \(N_{\rm tot}^{m,c}\) for the different components and molecules are taken from the user defined molfit file.

If local-overlap of lines is not taken into account (LocalOverlapFlag = False) the optical depth \(\tau_l^{m,c}(\nu)\) of all lines for each molecule \(m\) and component \(c\) is described as

\[\tau_l^{m,c}(\nu) = \sum_t \tau_t^{m,c}(\nu),\]

where the sum with index \(t\) runs over all spectral line transitions of molecule \(m\) within the given frequency range. The calculation procedure including local-overlap is described in Sect. “Local-overlap of neighboring lines”.

The beam-averaged continuum background temperature \(I_{\rm bg} (\nu)\) is parameterized as

(8)\[I_{\rm bg} (\nu) = T_{\rm bg} \cdot \left(\frac{\nu}{\nu_{\rm min}} \right)^{T_{\rm slope}} + J_\mathrm{CMB} + I_{\rm bg}^{\rm file} (\nu)\]

to allow the user to define the continuum contribution for each frequency range, individually. Here, \(\nu_{\rm min}\) indicates the lowest frequency of a given frequency range. \(T_{\rm bg}\) and \(T_{\rm slope}\), defined by the user, describe the background continuum temperature and the temperature slope, respectively. Additionally, XCLASS offers the possibility to describe the background continuum, i.e. the continuum seen by the layer close to the cosmic microwave background, by a so-called background file \(I_{\rm bg}^{\rm file} (\nu)\).

Finally, the last term \(J_\mathrm{CMB}\) in Eq. (1) describes the OFF position for single dish observations (defined by Inter_Flag = False) where we have an intensity caused by the cosmic background \(J_\mathrm{CMB}\). For interferometric observations, the contribution of the cosmic background is filtered out and has to be subtracted as well.

Hence, the layers, which do not belong to the largest distance parameter, have to be considered in an iterative manner where the solution of the radiative transfer equation for a certain distance \(i>1\) can be expressed as

(9)\[T_{\rm mb}^i(\nu) = \sum_{m,c \in i} \left(S^{m,c}(\nu) - T_{\rm mb}^{i-1}(\nu)\right) \, \left(1 - e^{-\tau_{\rm total}^{m,c}(\nu)}\right) + T_{\rm mb}^{i-1}(\nu) e^{-\tau_{\rm total}^{m,c}(\nu)},\]

where \(m\) and \(c\) indicates the index of the current molecule and component, respectively. Here, the expression \(T_{\rm mb}^{\rm i=1}(\nu)\), described by Eq. (1), represents the spectrum caused by layers which are closest to the background, including the beam-averaged continuum background temperature \(I_{\rm bg}^{\rm core} (\nu)\). The contribution by other components (\(i > 1\)) is considered by first calculating \(T_{\rm mb}^{\rm i=1}(\nu)\) and then use this as new continuum for lines at distance \(i\).

By fitting all species and their components at once, line blending and optical depth effects are taken into account. The modeling can be done simultaneously with isotopologues (and higher vibrational states) of a molecule assuming an isotopic ratio stored in the so-called iso ratio file, see Sect. “The iso ratio file”. Here, all parameters are expected to be the same except the column density which is scaled by one over the isotopic ratio for each isotopologue. Additionally, it is assumed that radiation emitted by all isotopologues of a molecule in a component interact with all other isotopologues, but the radiation emitted in one component does not interact with other molecules or with the same molecule in different components, i.e. their intensities are added, except local-overlap is taken into account, see Sect. “Local-overlap of neighboring lines”.

In order to correctly take instrumental resolution effects into account in comparing the modeled spectrum with observations, myXCLASS integrates the calculated spectrum over each channel. Thereby, myXCLASS assumes that the given frequencies \(\nu\) describe the center of each channel, respectively. The resulting value is than given as

(10)\[T_{\rm mb}(\nu) = \frac{1}{\Delta_c} \int_{\nu-\frac{\Delta_c}{2}}^{\nu+\frac{\Delta_c}{2}} T_{\rm mb}(\tilde{\nu}) d\tilde{\nu},\]

where \(\Delta_c\) represents the width of a channel. Due to the complexity of Eqn. (1), (9) the integration in Eq. (10) can not be done analytically. Therefore, myXCLASS performs a piecewise integration of each component and channel using the trapezoidal rule and summaries the resulting values to get the final value used in Eq. (10). In order to reduce the computation effort myXCLASS determines the minimal line width of all components of a certain distance. If a channel contains no transition, myXCLASS determines the intensities at \(\nu \pm \frac{\Delta_c}{2}\) and evaluates Eq. (10) by applying the trapezoidal rule. But, if transitions are included, myXCLASS re-samples the corresponding channel, whereby the number \(n_\nu\) of additional frequency points is given by

\[n_\nu = (m \cdot 3) \cdot \min \left[1, {\rm int} \left[\frac{\Delta_c}{\sigma_{\rm min}} \right] \right],\]

where \(m\) represents the number of transitions contained in the current channel and \(\sigma_{\rm min} = \frac{\Delta v_{\rm min}}{2 \, \sqrt{2 \, \ln \, 2} \, c_{\rm light}} \, \nu\). If local-overlap is taken into account, see Sect. “Local-overlap of neighboring lines”, \(\Delta v_{\rm min}\) indicates the minimal line width of all components of the current layer, otherwise \(\Delta v_{\rm min}\) describes the line width of the current component. Here, the function int(\(x\)) gives the integer part of \(x\). Finally, myXCLASS computes the intensities at these frequencies and evaluates Eq. (10).

The molfit file

Within the molfit file the user defines both which molecules (or radio recombination lines (RRLs)) are taken into account and how many components are used. The definition of parameters for a molecule (or RRLs) starts with a line describing the name of the molecule 8 (or RRLs), followed by the number of components \(N\). The following \(N\) lines describe the parameters for each components, separately. The number of columns and their meaning depends on the modeled contribution (molecule, or RRL). Generally, all parameters have to be separated by blanks, comments are marked with the % character.

Molecules

For molecules, the user has to define for each component the excitation temperature \(T_{\rm ex}\) in K (T_rot), the column density \(N_{\rm tot}\) in cm\(^{-2}\) (N_tot), and the velocity offset (relative to v\(_{\rm LSR}\)) in km s\(^{-1}\) (V_off). Depending of the used source description, see Sect. “Source description”, the size of the source \(\theta^{m,c}\) in arcsec (s_size) and their offset relative to the beam center is required as well. Furthermore, the Gaussian line profile requires the velocity width (FWHM) \(\Delta \nu\) in km s\(^{-1}\) (V_width). Additionally, the stacking of the components along the line-of-sight requires the cf-flag (CFFlag), see Sect. “Simple 1-d description”, or the distance parameter (distance), see Sect. “Stacking”. Finally, the column (keyword) describes so-called keywords, that are used for the sub-beam description, among other things. Note, different keywords have to be separated by “_” character.

Example of a molfit file describing two molecules CS\(_{v=0}\) and HCS\(^+_{v=0}\) with two components, respectively:

% Number of molecules = 2
% s_size:  T_rot:    N_tot:  V_width:   V_off:  CFFlag:   keyword:
CS;v=0;   3
   48.470  300.00  3.91E+17      2.86  -20.564        c
   21.804  320.00  6.96E+17      8.07   30.687        c
   81.700  208.00  1.46E+17      5.16  -10.124        c
HCS+;v=0;   2
% s_size:  T_rot:    N_tot:  V_width:   V_off:  CFFlag:   keyword:
           150.00  1.10E+18      5.00   -0.154        f
           200.00  2.20E+17      3.10   -2.154        f

Limiting the number of transitions

XCLASS offers the possibility to limit the number of transitions of molecules which enter Eq. (1). Therefore, the user has to add the following command words to the beginning of the molfit file:

"%%MinNumTransitionsSQL": defines the minimum number of transitions, i.e. XCLASS takes only those molecules into account which have at least this number of transition in the given frequency ranges.
"%%MaxNumTransitionsSQL" describes the max. number of transitions (in conjunction with "%%TransOrderSQL"). Here, "0" means all transitions are included.
"%%TransOrderSQL" (in conjunction with "%%MaxNumTransitionsSQL") defines the order of transition. (E.g. "%%MaxNumTransitionsSQL = 10" and "%%TransOrderSQL = 1" means, that only the first 10 transitions with the lowest lower energy are taken into account.)
"%%MaxElowSQL" defines the upper limit of the lower energy of a transition.
"%%MingASQL" describes the lower limit of gA (upper state degeneracy * Einstein A coefficient) of a transition.

Please note, for recombination lines only the command words "%%MinNumTransitionsSQL" and "%%MaxNumTransitionsSQL" are taken into account.

Example for a molfit file containing limiting transition parameters:

%==============================================================================
%
% define transition parameters:
%
%==============================================================================
%%MinNumTransitionsSQL = 1      % min. number of transitions
%%MaxNumTransitionsSQL = 0      % max. number of transitions,
                                % =0): all transitions are taken into account
%%TransOrderSQL = 1             % order of transitions:
                                %  (=1): by lower energy,
                                %  (=2): by gA,
                                %  (=3): gA/E_low^2, else by trans. freq.
%%MaxElowSQL = 1.0e+06          % max. lower energy, i.e. upper limit of
                                %  the lower energy of a transition
%%MingASQL = 0.0                % minimal intensity, i.e. lower limit of
                                %  gA of a transition


%==============================================================================
%
% define molecules and their components:
%
%==============================================================================
% s_size:   T_rot:    N_tot:  V_width:   V_off:  CFFlag:   keyword:
CS;v=0;   3
    48.470  300.00  3.91E+17      2.86  -20.564        c
    21.804  320.00  6.96E+17      8.07   30.687        c
    81.700  208.00  1.46E+17      5.16  -10.124        c
HCS+;v=0;   2
% s_size:   T_rot:    N_tot:  V_width:   V_off:  CFFlag:   keyword:
            150.00  1.10E+18      5.00   -0.154        f
            200.00  2.20E+17      3.10   -2.154        f

The iso ratio file

In order to reduce the number of model parameters, XCLASS offers the possibility to use a so-called iso ratio file. Here, XCLASS assumes, that both molecules (isotopologue and iso master molecule) are described by the same number of components, where the source size, the rotation temperature, the line width, and the velocity offset are identical. The corresponding column densities of the isotopologue \(\left({\rm N}_{\rm tot}^{\rm (isotopologue)}\right)\) are scaled by the so-called iso ratio, i.e.

\[{\rm N}_{\rm tot}^{\rm (isotopologue)} = {\rm N}_{\rm tot}^{\rm (iso \,\, master)} / {\rm (iso \, \, ratio)}.\]

The iso ratio file contains three columns separated by tabs or at least two blank characters, where the first two columns indicates the isotopologue and the corresponding iso master molecule, respectively. The third column defines the ratio for both molecules. Comments are marked with a % character, i.e. all characters on the right side of a % are ignored.

Example for an iso ratio file:

% isotopologue:         molecule:        ratio:
HC-33-S+;v=0;           HCS+;v=0;          75.0
HC-34-S+;v=0;           HCS+;v=0;          22.5
CS;v=4;                 CS;v=0;             2.3
CS;v=3;                 CS;v=0;             2.1
CS;v=2;                 CS;v=0;             2.1
CS;v=1;                 CS;v=0;             2.0
CS-34;v=0;              CS;v=0;            22.5
CS-34;v=1;              CS;v=0;            22.5
CS-33;v=1;              CS;v=0;            75.0
CS-33;v=0;              CS;v=0;            75.0

In the example described above we define in the first line an iso ratio of 75.0 between HC-33-S+;v=0; and HCS+;v=0;. Assuming that HCS+;v=0; is described with one component and with a column density of \(9 \cdot 10^{16}\) cm\(^{-2}\) the column density of HC-33-S+;v=0; is

\[{\rm N}_{\rm tot} = 9 \cdot 10^{16} \, {\rm cm}^{-2} / 75.0 = 1.2 \cdot 10^{15} \, {\rm cm}^{-2}.\]

Globally defined iso-ratios

In addition, the user can defined so-called globally defined iso ratios, e.g. (\(^{12}\)C / \(^{13}\)C). This defined ratio is multiplied to all ratios of isotopologues in the iso ratio file which contain e.g. \(^{13}\)C raised to the power of occurrences of the globally defined ion in each isotopologue. For example, the ratio (\(^{12}\)C / \(^{13}\)C) is set globally to 45. Additionally, we set the ratio of (CH\(_3\)CN,v=0 / \(^{13}\)CH\(_3 ^{13}\)CN,v=0) to 3.0. The final iso ratio used by XCLASS for (CH\(_3\)CN,v=0 / \(^{13}\)CH\(_3 ^{13}\)CN,v=0) is then given by \(3.0 \times 45^2 = 6075\). (Here, the exponent \(2\) is caused by the two \(^{13}\)C ions). Additionally, globally defined iso-ratios can be combined as well. In order to determine the used iso ratio for e.g. (\(^{12}\)C\(^{32}\)S)/(\(^{13}\)C \(^{33}\)S) we just have to multiply the globally defined is ratios for (\(^{12}\)C / \(^{13}\)C) and (\(^{32}\)S / \(^{33}\)S).

Globally defined ratios in the iso ratio file are indicated by the phrase “GLOBAL__” and can not be used alone without other iso ratios, e.g.

% isotopologue:         molecule:           ratio:
GLOBAL__C-13            C-12                45.0
HO-18;v=0;              OH;v=0;             500
HDO-18;v=0;             HDO;v=0;            500
C-13-H3C-13-N;v=0;      CH3CN;v=0;          1.0

Fitting iso ratios in the iso ratio file

The myXCLASSFit function offers the possibility to optimize the ratios between isotopologues as well. For that purpose, the user has to add two additional columns on the right indicating the lower and the upper limit of a certain ratio.

Example for an iso ratio file where the ratios are optimized by the myXCLASSFit function:

% isotopologue:         molecule:        ratio:     lower:     upper:
CS;v=4;                 CS;v=0;             2.3        0.1      500.0
CS;v=3;                 CS;v=0;             2.1        0.1      500.0
CS;v=2;                 CS;v=0;             2.1        0.1      500.0
CS;v=1;                 CS;v=0;             2.0        0.1      500.0
CS-33;v=1;              CS;v=0;            75.0        0.1      500.0
CS-33;v=0;              CS;v=0;            75.0        0.1      500.0
HC-33-S+;v=0;           HCS+;v=0;          75.0        0.1      500.0
HC-34-S+;v=0;           HCS+;v=0;          22.5        0.1      500.0

If the lower and upper limit are equal or if the lower limit is higher than the upper limit, the ratio is kept constant and is not optimized by the myXCLASSFit function. Note, if either the iso master molecule or a corresponding isotopologue has no transition within at least one given frequency range, the myXCLASSFit function does not optimize the corresponding iso-ratio. For example, if the isotopologue HNC-13;v=0; (used in the example described above) has no transition in at least one given frequency range, the given iso-ratio (here 60) is kept constant. Additionally, if a iso master molecule has no transition within at least one given frequency range, the iso-ratios to all of its isotopologues are kept constant.

Local-overlap of neighboring lines

XCLASS offers the possibility to take the local overlap of neighboring (molecules and recombination) lines into account. Here, we follow the derivation described by [Cesaroni_Walmsley_1991] and define an average source function \(S_l\) at frequency \(\nu\) and distance \(l\)

(11)\[S_l (\nu) = \frac{\varepsilon_l (\nu)}{\alpha_l (\nu)} = \frac{\sum_{c, t \in l} \tau_t^c (\nu) \, S_\nu (T_{{\rm rot}}^c)}{\sum_{c, t \in l} \tau_t^c (\nu)},\]

where \(\varepsilon_l\) represents the emission and \(\alpha_l\) the absorption function, respectively. Additionally, \(T_{\rm rot}^c\) indicates the excitation temperature and \(\tau_t^c\) the optical depth of transition \(t\) and component \(c\). For each frequency channel, we take all components \(c\) and transitions \(t\) into account, which belongs to the current distance and whose Doppler-shifted transitions frequencies are located within a range of 5\(\Delta v_{\rm max}\), where \(\Delta v_{\rm max}\) describes the largest line width of all components of the current layer. The iterative treatment of components at different distances, takes non-local effects into account as well. In the optically thin limit, Eq. (11) is equal to the traditional approach of convolving a line with several Gaussians. But it better describes situations in which both optically thin and optically thick lines are present: photons emitted from the optically thin transition are locally absorbed by the optically thick emission. The intensities of the lines do not simply add up like in the optically thin limit, but the intensity at the overlapping frequencies is mainly described by the optically thick emission.

Please note, that the myXCLASS function can not compute the intensity and corresponding optical depth for each component if the local overlap is taken into account. Only the intensities and optical depths at each distance can be determined.

Source description

The position and size of each component in the molfit file can be defined in two different ways: A simple 1-d description with a very simplistic geometrical structure offers a fast but not very accurate and flexible description with a minimum of parameters, wheres the advanced 3-d description is more accurate but requires more parameters and more computational effort. The different source descriptions can be applied to all kinds of components, i.e. to molecules and RRLs.

Fig. 1 Sketch of a distribution of core layers within the Gaussian beam of the telescope (black ring). Here, we assume three different core components \(1\), \(2\), and \(3\), centered at the middle of the beam with different source sizes, excitation temperatures, velocity offsets (relative to v\(_{\rm LSR}\)) etc. indicated by different colors. Additionally, we assume that all core components have the same distance to the telescope, i.e. all core layers are located within a plane perpendicular to the line of sight. Furthermore, we assume that this plane is located in front of a background layer \(4\) with homogeneous intensity \(I_{\rm bg}^{\rm core}(\nu)\) over the whole beam. Core components do not interact with each other radiatively.

Simple 1-d description

For the 1-d description the molfit file has to contain column (CFFlag) indicating the core and foreground flag and must not contain column (distance) describing the stacking parameter, see Sect. “Stacking”. The 1-d assumption imposes a very simplistic geometrical structure where we recognize two classes of components:

One, the core objects (in earlier implementations called emission component), consists of an ensemble of objects centered at the middle of the beam. These could be identified with clumps, hot dense cores etc. which overlaps but do not interact either because they do not overlap in physical or in velocity space. For computational convenience, they are assumed to be centered in the beam, as shown in Fig. 1. It is also assumed that the dust emission emanates (partly) from these components. Their intensities are added, weighted with the beam filling factor, see Eq. (2). Note, all core components are contained in the layer which is closest to the background.

The second class, foreground objects (in earlier implementations called absorption components), are assumed to be in layers in front of the core components. In the current 1-d implementation, they would have a beam-averaged intensity of the core sources as background, and would fill the whole beam. Examples for such structures would be source envelopes in front of dense cores, or foreground components along the line-of-sight.

As shown in Fig. 1, we assume that core components do not interact with each other radiatively, i.e. one core layer is not influenced by the others, if local-overlap is not taken into account. But the core layers may overlap to offer the possibility to model sources consisting of several molecules and compounds. For core components, Eq. (1) has to be slightly modified 9, so that the solution of the radiative transfer equation for core layers is 10,

(12)\[\begin{split}T_{\rm mb}^{\rm core}(\nu) = \sum_m \sum_c &\Bigg[\eta \left(\theta^{m,c}\right) \left[S^{m,c}(\nu) \left(1 - e^{-\tau_{\rm total}^{m,c}(\nu)}\right) + I_{\rm bg}^{\rm core} (\nu) \left(e^{-\tau_{\rm total}^{m,c}(\nu)} - 1\right) \right] \Bigg] \\ &+ \left(I_{\rm bg}^{\rm core} (\nu) - J_\mathrm{CMB} \right),\end{split}\]

where the sums go over the indices \(m\) for molecule, and \(c\) for (core) component, respectively. The different terms contained in Eq. (12) are described in the previous sections.

Fig. 2 Sketch of a distribution of core and foreground layers within the Gaussian shaped beam of the telescope (black ring). Here, we assume three different core components \(2a\), \(2b\), and \(2c\) located in a plane perpendicular to the line of sight which lies in front of the background layer \(1\) with intensity \(I_{\rm bg}^{\rm core}(\nu)\), see Eq. (8). The foreground layers \(3\) and \(4\) are located between the core layers and the telescope along the line of sight (black dashed line). Here, each component is described by different excitation temperatures, velocity offsets etc. indicated by different colors. The thickness of each layer is described indirectly by the total column density \(N_{\rm tot}^{m,c}\), see Sect. “Optical Depth”. For each foreground layer we assume a beam filling factor of 1, i.e. each foreground layer covers the whole beam.|

In contrast to core layers, foreground components may interact with each other, as shown in Fig. 2, where absorption takes places only, if the excitation temperature for the absorbing layer is lower than the temperature of the background.

Hence, the solution of the radiative transfer equation for foreground layers can not be given in a form similar to Eq. (12). Foreground components have to be considered in an iterative manner similar to Eq. (9). The solution of the radiative transfer equation for foreground layers can be expressed as

(13)\[\begin{split}T_{\rm mb}^{\rm fore}(\nu)_{m,c=1} &= \eta \left(\theta^{m,c=1} \right) \left(S^{m,c=1}(\nu) - T_{\rm mb}^{\rm core}(\nu)\right) \, \left(1 - e^{-\tau_{\rm total}^{m,c=1}(\nu)}\right) + T_{\rm mb}^{\rm core}(\nu) \\ T_{\rm mb}^{\rm fore}(\nu)_{m,c=i} &= \eta \left(\theta^{m,c=i} \right) \left(S^{m,c=i}(\nu) - T_{\rm mb}^{\rm fore}(\nu)_{m,c=(i-1)}\right) \, \left(1 - e^{-\tau_{\rm total}^{m,c=i}(\nu)}\right) + T_{\rm mb}^{\rm fore}(\nu)_{m,c=(i-1)},\end{split}\]

where \(m\) indicates the index of the current molecule and \(i\) represents an index running over all foreground components \(c\) of all molecules. Additionally, we assume that each foreground component covers the whole beam, i.e. \(\eta \left(\theta^{m,c=i} \right) \equiv 1\) for all foreground layer 11. Thus, Eq. (13) simplifies to

(14)\[\begin{split}T_{\rm mb}^{\rm fore}(\nu)_{m,c=1} &= \Big[ S^{m,c=1}(\nu) \left(1 - e^{-\tau_{\rm total}^{m,c=1}(\nu)}\right) + T_{\rm mb}^{\rm core}(\nu) e^{-\tau_{\rm total}^{m,c=1}(\nu)}\Big] \\ T_{\rm mb}^{\rm fore}(\nu)_{m,c=i} &= \Big[S^{m,c=i}(\nu) \left(1 - e^{-\tau_{\rm total}^{m,c=i}(\nu)}\right) + T_{\rm mb}^{\rm fore}(\nu)_{m,c=(i-1)} e^{-\tau_{\rm total}^{m,c=i}(\nu)}\Big],\end{split}\]

where \(T_{\rm mb}^{\rm core}(\nu)\) describes the core spectrum, see Eq. (12), including the beam-averaged continuum background temperature \(I_{\rm bg}^{\rm core} (\nu)\). For foreground lines the contribution by other components is considered by first calculating the contribution of core objects and then use this as new continuum for foreground lines reflecting the fact that cold foreground layers are often found in front of hotter emission sources. The myXCLASS function assumes, that the cosmic background describes together with the core components one end of a stack of layers. Additionally, the foreground components are located between this plane and the telescope, see Fig. 2. The total column density \(N_{\rm tot}^{m,c}\) depends on the abundance of a certain molecule and on the thickness of a layer containing the molecule. The order of components along the line of sight is defined by the occurrence of a certain foreground component in the molfit file.

Stacking

Fig. 3 Sketch of components located at different distances along the line of sight.

In addition to the simple 1-d description described in the previous section, XCLASS offers the possibility to define the order of each component along the line-of-sight by a so-called stacking or distance parameter which replaces column CFFlag in the molfit file, i.e. the column CFFlag is replaced by column distance right after column V_off. The stacking parameter can be interpreted as distance of a component to the observer where the component(s) with the largest distance parameter is (are) closest to the background while those components with the smallest distance parameter are closest to the observer. Additionally, the components which have the same distance parameter (i.e. within a difference of \(< 10^{-9}\)) are located in the same plane perpendicular to line-of-sight. In contrast to foreground components, described in the previous section, which can describe only one molecule at a certain distance, the stacking formalism offers the possibility to describe an arbitrary number of clouds consisting of an arbitrary number of molecules or recombination lines along the line-of-sight.

Note, the distance parameter has no unit, can not be optimized by MAGIX, and must not be used in conjunction with the old cf-flag formalism, described in the previous section.

In the following example SO\(_2\) is described by four and CH\(_3\)CN with two components where the first components of SO\(_2\) and CH\(_3\)CN are located close to the background, i.e. both components are located in the same plane perpendicular to the line of sight and are illuminated by the background continuum Eq. (8). The second and third component of SO\(_2\) are located in another plane perpendicular to the line of sight, which is closer to the observer. Both components “see” a background which is caused by the background plus the first components of SO\(_2\) and CH\(_3\)CN. The fourth component of SO\(_2\) and the second component of CH\(_3\)CN are located in front of the plane at 500, where the second component of SO\(_2\) is closer to the observer.

SO2;v=0;   4

% s_size:  T_rot:   N_tot:   V_width:   V_off: distance:  keyword:
%[arcsec]     [K]   [cm-2]     [km/s]   [km/s]        []

   500.00   50.00  1.9E+14       2.86     1.22     600.0
   500.00  190.00  3.9E+15       4.16    -1.61     500.0
   500.00  250.00  3.9E+16       8.55    -5.75     500.0
   500.00   80.00  3.9E+14       1.23    -1.82     300.0

CH3CN;v=0;   2

% s_size:  T_rot:   N_tot:   V_width:   V_off: distance:  keyword:
%[arcsec]     [K]   [cm-2]     [km/s]   [km/s]        []

   500.00   50.00  2.2E+15       2.86     6.01     600.0
   500.00   80.00  8.5E+14       1.23    -2.82     200.0

Sub-beam description

Fig. 4 Sketch of a sub-beam description for a single frequency channel. Here, the gray grid indicates the model grid, the black ellipse the beam, and the red dot the center of the beam, respectively. Additionally, the filled circles indicate different components, where the gray scale represent the location of a component along the line-of-sight, i.e. light gray describe components close the cosmic microwave background, while dark gray represents components close the observer. The projected width \(w\) and height \(h\) of the beam are described by Eq. (15).

Besides to the aforementioned component stacking XCLASS offers the possibility to take the shape and the position of each component perpendicular to the line-of-sight into account as well, see Fig. 3. The so-called sub-beam description, see Fig. 4, is used whenever an offset position of a component is defined, or if a component at a smaller distance is smaller than a component at larger distance, or if local-overlap is taken into account. Each component can have a different size and center position. In the current version only circle and square-shaped components can be described. Additionally, elliptical rotated beam sizes, see Sect. “Elliptical rotated Gaussians”, can be defined as well, but the frequency dependence of single dish beam sizes is ignored, i.e. the beam size is assumed to be constant for the whole frequency range 12. For the sub-beam description, XCLASS samples the elliptical rotated beam by an rectangular linear grid, called model grid, where the number of pixels along the x- and y-direction is defined by the user. The size of each model pixel is defined by the beam, where we assume that the grid along the x- and y-axis is three times the width \(w\) and height \(h\) of the beam, respectively, described by

(15)\[\begin{split}w &= \sqrt{\left({\tt BMIN} \cdot \cos {\tt BPA} \right)^2 + \left({\tt BMAJ} \cdot \sin {\tt BPA} \right)^2},\\ h &= \sqrt{\left({\tt BMIN} \cdot \sin {\tt BPA} \right)^2 + \left({\tt BMAJ} \cdot \cos {\tt BPA} \right)^2}.\end{split}\]

Here, \({\tt BMIN}\) and \({\tt BMAJ}\) represents the minor and major axis and \({\tt BPA}\) the position angle of the beam. Furthermore, XCLASS identifies for each component the model pixels, which are covered by the corresponding component. While the identification of all pixels covered by a component is trivial for square components, the Bresenham’s circle algorithm [Bresenham_1977] is used for circular components. After this, XCLASS computes the spectra for each model pixel based on the corresponding component configuration. Afterwards, XCLASS convolves the calculated intensity map for each frequency channel with the elliptical rotated Gaussian beam using the FFTW library 13 [Padua_2011]. Finally, XCLASS uses the spectrum at the beam center for the final output. The intensities and optical depths for each component / distance are computed in the same way.

For beam-centered components only the diameter of the component (in arcsec) is required which has to be defined in the first column s_size of the molfit file, e.g.

SO2;v=0;   4

% s_size:  T_rot:   N_tot:   V_width:   V_off: distance:  keyword:
%[arcsec]     [K]   [cm-2]     [km/s]   [km/s]        []

    50.00   50.00  1.9E+14       2.86     1.22     600.0
   110.00  190.00  3.9E+15       4.16    -1.61     500.0

Components whose center are not identical to the center of the beam require three additional parameters: The column s_off_x and s_off_y describe the x- and y-coordination of the component center position relative to the beam center (in arcsec), respectively. Additionally, column (keyword) has to contain the keyword circleoff.

In the following example, we describe SO\(_2\) with two components, where the first component has a size (diameter) of 50 arcsec, and a relative offset of 1.0 and 1.45 arcsec, respectively. The second component has a larger diameter of 200 arcsec and is centered at the middle of the beam.

SO2;v=0;   2

% s_size: s_off_x: s_off_y: T_rot: N_tot: V_width:   V_off: distance: keyword:
%[arcsec] [arcsec] [arcsec]    [K] [cm-2]   [km/s]   [km/s]        []

    50.00      1.0    1.45   50.00  1.e14     2.86     1.22     600.0 circleoff
   200.00                    90.00  3.e15     4.16    -1.61     500.0

Additionally, the user can use keywords box and boxoff to define quadratic components. For keyword box, XCLASS assumes that the corresponding component is located at the center of the beam, where column s_size indicates the length (in arcsec) of the square. Similar to keyword circleoff the keyword boxoff defines a square which center (in arcsec) is defined by columns s_off_x and s_off_y, respectively.

Radio recombination lines (RRLs)

In addition to molecules myXCLASS can deal with radio recombination lines (RRLs) as well.

RRLs in LTE

The optical depths of RRLs (in LTE) is given by

(16)\[\begin{split}\tau_{{\rm RRL}, \nu} &= \int \kappa_{n_1, n_2, \nu}^{\rm ext} \, ds \\ &= \frac{\pi \, h^3 \, e^2}{(2 \pi \, m_e \, k_B)^{3/2} \, m_e \, c} \cdot {\rm EM} \cdot \frac{n_1^2 \, f_{n_1, n_2}}{T_e^{3/2}} \, \exp \left[\frac{Z^2 \, E_{n_1}}{k_B T_e} \right] \, \left(1 - e^{-h \nu_{n_1, n_2} / k_B T_e} \right) \, \phi_\nu,\end{split}\]

where the LTE source function \(S_{{\rm RRL}, \nu}^{\rm LTE}\) can be written in terms of the brightness temperature Eq. (4) as

(17)\[S_{{\rm RRL}, \nu}^{\rm LTE} = \frac{\tau_{{\rm RRL}, \nu} \cdot J(T_e)}{\tau_{{\rm RRL}, \nu}}.\]

The terms \(n_1\), \(f_{n_1, n_2}\), \(E_{n_1}\), and \(\nu_{n_1, n_2}\) are taken from the embedded database whereas the emission measure EM (in cm\(^{-6}\) pc), electronic temperature \(T_e\) (in K), source size \(\theta^{m,c}\) (in arcsec) (or other parameters related to source description), line width \(\Delta v\) (in km/s), and velocity offset \(v\) (in km/s) used for the line profile function \(\phi_\nu\) are taken from the user-defined molfit file.

The parameters which are used for the description of radio recombination lines (RRLs) are quite similar to the parameters used for molecules. Here, the kinetic temperature is interpreted as electronic temperature (in K) and the column density as emission measure EM (in pc cm\(^{-6}\)). RRLs of hydrogen are marked as RRL-H, of Helium as RRL-He, of Carbon as RRL-C, of Nitrogen as RRL-N, and of Oxygen as RRL-O. In the following example we describe RRLs of hydrogen at two distances with a Gaussian line shape:

RRL-H    2
% s_size:     T_e:   EM_RRL:  V_width:   V_off:  distance:
500.00000    6.9e4   5.5E+16       1.3      2.0      2.e15
500.00000    2.3e4   2.3E+12       2.5      2.0      1.e15

Gaussian line profile

In order to take broadening caused by the thermal motion of the gas particles and micro-turbulence into account a normalized Gaussian line profile, i.e. \(\int_0^{\infty} \phi^{m,c,t}_{\rm Gauss}(\nu) \, d\nu\) = 1, for a spectral line \(t\) can be assumed:

\[\phi^{m,c,t}_{\rm Gauss}(\nu) = \frac{1}{\sqrt{2 \pi} \, \sigma^{m,c,t}} \cdot e^{-\frac{\left(\nu - \nu_0^{m,c,t} \right)^2} {2 (\sigma^{m,c,t})^2}},\]

where \(\nu_0^{m,c,t}\) indicates the Doppler-shifted transition frequency and \(\sigma^{m,c,t}\) the standard deviation which is related to the full width at half maximum \(f_G\) by

(18)\[f_G = \left(2 \, \sqrt{2 \, \ln 2}\right) \cdot \sigma^{m,c,t} = \frac{\Delta {\rm v}_{\rm width-Gauss}^{m,c}}{c_{\rm light}} \cdot \left(\nu^t + \delta \nu_{\rm LSR}^{m,c,t} \right).\]

Similar to the Lorentzian line profile \(\Delta {\rm v}_{\rm width-Gauss}^{m,c}\) and \(\delta {\rm v}_{\rm offset}^{m,c}\) indicates the line width and offset of molecule \(m\) and component \(c\), respectively and \({\rm v}_{\rm LSR}\) the source velocity.

The myXCLASS function assumes a Gaussian line shape for all components by default.

Continuum

Dust contribution

The dust optical depth \(\tau_d^{m,c}(\nu)\) takes the dust attenuation into account and is given by

(19)\[\begin{split}\tau_d^{m,c}(\nu) &= \tau_{d, {\rm ref}}^{m,c} \cdot \left(\frac{\nu}{\nu_{\rm ref}} \right)^{\beta^{m,c}} + \tau_d^{\rm file}(\nu) \\ &= \left(N_H^{m,c} \cdot \kappa^{m,c}_{\nu_{\rm ref}} \cdot m_{H_2} \cdot \frac{1}{\zeta_{\rm gas-dust}}\right) \cdot \left(\frac{\nu}{\nu_{\rm ref}} \right)^{\beta^{m,c}} + \tau_d^{\rm file}(\nu).\end{split}\]

Here, \(N_H^{m,c}\) describes the hydrogen column density (in cm\(^{-2}\)), \(\kappa^{m,c}_{\nu_{\rm ref}}\) the dust mass opacity for a certain type of dust (in cm\(^{2}\) g\(^{-1}\)) [Ossenkopf_Henning_1994], and \(\beta^{m,c}\) the spectral index. Additionally, \(\nu_{\rm ref}\) = 230 GHz indicates the reference frequency for \(\kappa^{m,c}_{\nu_{\rm ref}}\), \(m_{H_2}\) the mass of a hydrogen molecule, and \(1 / \zeta_{\rm gas-dust}\) the dust to gas ratio which is set here to (1/100) [Hildebrand_1983]. The equation is valid for dust and gas well mixed. Furthermore, the expression \(\tau_d^{\rm file}(\nu)\) represents the dust optical depth described in an ASCII file, which path and name is defined by the input parameter DustFileName.

In Eqn. (1) and (9), the dust continuum contribution is described through the source function \(S^{m,c}(\nu)\), Eq. (5), by an effective dust temperature \(T_d^{m,c}\) (in K) (through \(J(T_d^{m,c}, \nu)\)).

The myXCLASS function offers three different ways to define the dust parameters \(N_H^{m,c}\), \(\kappa^{m,c}_{\nu_{\rm ref}}\), \(\beta^{m,c}\), and \(T_d^{m,c}\):

On the one hand the dust parameters \(N_H^{m,c}\), \(\kappa^{m,c}_{\nu_{\rm ref}}\), and \(\beta^{m,c}\) can be defined globally for all components in the molfit file by the input parameters N_H, kappa_1300, and beta_dust of the myXCLASS function, while the dust temperature is set to the excitation temperature of the corresponding component. For the myXCLASSFit, myXCLASSMapFit, myXCLASSMapRedoFit, and LineID function the dust parameters can be defined globally in an obs. xml file as well, see Sect. “Observational xml file”.

Additionally, the dust parameters can be described in the molfit file for different distances, while the dust continuum has to be marked as cont-dust. Each component is described by the size of the dust component s_size (in arcsec), the dust temperature T_cont_dust (in K), the hydrogen column density nHcolumn_cont_dust (in cm\(^{-2}\)), the dust mass opacity kappa_cont_dust (in g cm\(^{-2}\)), the spectral index beta_cont_dust, and the distance parameter distance. In the example described below, we define four different dust components with different dust temperatures at different distances

cont-dust    4
% s_size: T_cont_dust: nHcolumn_cont_dust: kappa_cont_dust: beta_cont_dust: distance:
500.00000       50.000             3.9E+16             0.42             2.0    4.e15
500.00000       10.000             1.2E+20             0.42             2.0    3.e15
500.00000       32.000             5.3E+12             0.42             2.0    2.e15
500.00000        5.000             7.1E+18             0.42             2.0    1.e15

In order to be backward compatible with older XCLASS releases, the user can define the dust parameters for each component of a molecule as well. Here, the dust temperature \(T_d^{m,c}\) is parameterized as offset to the corresponding excitation temperature \(T_{\rm ex}^{m,c}\) as

(20)\[\begin{split}T_d^{m,c} (\nu) &= T_{\rm ex}^{m,c} + \Delta T_d^{m,c}(\nu) \\ &= T_{\rm ex}^{m,c} + T_{d, {\rm off}}^{m,c} \cdot \left(\frac{\nu}{\nu_{\rm min}} \right)^{T_{d, {\rm slope}}^{m,c}},\end{split}\]

where \(T_{d, {\rm off}}^{m,c}\) and \(T_{d, {\rm slope}}^{m,c}\) can be defined by the user for each component in the molfit. If \(T_{d, {\rm off}}^{m,c}\) and \(T_{d, {\rm slope}}^{m,c}\) are not defined for a certain component, we assume \(T_d^{m,c} (\nu) \equiv T_{\rm ex}^{m,c} (\nu)\) for all components. For a physical (\(\gamma \equiv 1\)) description of the background intensity, see Eq. (5), the user can define the dust opacity, Eq. (19), and dust temperature, Eq. (20), for each component. In order to define a dust temperature for each component which is not identical to the excitation temperature \(T_{\rm ex}\) of the corresponding component, the molfit file has to contain two additional columns between columns V_off and CFFlag (or column distance), describing the parameters \(T_{\rm d, off}^{m,c}\) (T_doff) and \(T_{\rm d, slope}^{m,c}\) (T_dSlope), respectively, and keyword t-dust-inline in column (keyword), e.g.

CS;v=0;    3
% s_size: T_rot: N_tot: V_width: V_off: T_doff: T_dSlope: CFFlag: keyword:
%[arcsec]    [K] [cm-2]   [km/s] [km/s]     [K]        []   [c/f]
    48.47  50.00   3e16     2.86  -2.56     3.0       0.0       c t-dust-inline
    40.10  56.53   2e18     4.21  -7.31     3.0       0.0       c t-dust-inline
    29.09  68.44   2e17     4.02   0.21     2.5       1.0       c t-dust-inline

Additionally, the myXCLASS program offers the possibility to define \(N_H\), \(\kappa_{\nu_{\rm ref}}\) and \(\beta\) for each component of a molecule as well. In order to define these parameters individually for each component, the molfit file has to contain three additional columns on the left side of column CFFlag (or column distance), indicating parameters \(N_H^{m,c}\) (nHcolumn), \(\kappa^{m,c}_{\nu_{\rm ref}}\) (kappa), and \(\beta^{m,c}\) (beta), and keyword nh-dust-inline in column (keyword), e.g.

CS;v=0;    3
% s_size: T_rot: N_tot: V_width: V_off: nHcolumn:  kappa: beta: CFFlag: keyword:
%[arcsec]    [K] [cm-2]   [km/s] [km/s]    [cm-2] [cm2/g]    []   [c/f]
    48.47  50.00   3e16     2.86  -2.56      3e24    0.42   2.0       c nh-dust-inline
    56.53  20.00   2e18     4.21  -7.31      2e24    0.42   2.1       c nh-dust-inline
    68.44  90.00   2e17     4.02   0.21      3e21    0.42   1.9       c nh-dust-inline

In order to define \(T_{d, {\rm off}}^{m,c}\), \(T_{d, {\rm slope}}^{m,c}\), \(N_H\), \(\kappa_{\nu_{\rm ref}}\), and \(\beta\) for each component the columns defining the dust temperature\(T_{d, {\rm off}}^{m,c}\) (T_doff), \(T_{d, {\rm slope}}^{m,c}\) (T_dslope) have to be given before the hydrogen column density, kappa and beta, i.e. in the order of s_size, T_rot, N_tot, V_width, V_off, T_doff, T_dslope, nHcolumn, kappa, beta, and CFFlag. Additionally, the column (keyword) has to contain the keyword t-dust-inline_nh-dust-inline.

Phenomenological continuum description

In addition to the the dust continuum contribution described above, XCLASS offers the possibility to describe the continuum of a given distance by phenomenological functions. Currently, the following functions are included:

Standing wave (ContPhenFuncID = 1):

\[c(\nu) = A_1 \, \sin \left( (\nu \cdot x) + \phi_1 \right) + A_2 \, \cos \left( (\nu \cdot x) + \phi_2 \right),\]

where \(A_{1, 2}\) and \(x\) describe scale parameters whereas \(\phi_{1,2}\) indicate phase shifts. Here, \(A_1\) and \(A_2\) are defined by parameters ContPhenFuncParam1 and ContPhenFuncParam2, respectively. Parameter ContPhenFuncParam3 indicates the scale parameter \(x\) and the phase shifts \(\phi_{1,2}\) are described by parameters ContPhenFuncParam4 and ContPhenFuncParam5, respectively.
Planck function (ContPhenFuncID = 2):

\[c(\nu) = \frac{h \, \nu}{k_B}\frac{1}{e^{h \, \nu / k \, T} - 1},\]

with temperature \(T\) described by input parameter ContPhenFuncParam1. The other input parameters ContPhenFuncParam2 - ContPhenFuncParam5 have to be set to zero.
beam-averaged continuum background temperature (ContPhenFuncID = 3)

(21)\[c(\nu) = T_{\rm bg} \cdot \left(\frac{\nu}{\nu_{\rm ref}} \right)^{T_{\rm slope}},\]

where \(\nu_{\rm ref}\) (described by input parameter ContPhenFuncParam1) indicates the reference frequency for \(T_{\rm bg}\) (defined by input parameter ContPhenFuncParam2). The temperature slope \(T_{\rm slope}\) is defined by input parameter ContPhenFuncParam3. If input parameter ContPhenFuncParam4 is set to zero Eq. (21) is applied to all frequency ranges. But the user can restrict the application of Eq. (21) to a given frequency range by using input parameters ContPhenFuncParam1 and ContPhenFuncParam4, where ContPhenFuncParam1 defines the first and ContPhenFuncParam4 the last frequency, respectively. The other input parameter ContPhenFuncParam5 must always be set to zero.

Note, the aforementioned continuum function \(c(\nu)\) is added to the background intensity of the corresponding distance.

In order to use a phenomenological description of the continuum in the molfit file, the phrase cont-phen has to be used. Each component is described by the source size s_size (in arcsec) (or other parameters related to source description), the function-id func-id, the input parameters p_1, p_2, p_3, p_4, and p_5, representing the input parameters ContPhenFuncParam1 - ContPhenFuncParam5 and the distance distance. In the example describe below, we use an expression for a standing wave to describe the continuum phenomenologically at a given distance:

cont-phen    1
% s_size:   func-id:   p_1:   p_2:   p_3:   p_4:   p_5:   distance:
500.00000          1    0.2    0.2    0.2    0.2    0.2       1.e20

Similar to the dust contribution, the phenomenological continuum parameters can be defined globally, by using the input parameters ContPhenFuncID, ContPhenFuncParam1, ContPhenFuncParam2, ContPhenFuncParam3, ContPhenFuncParam4, and ContPhenFuncParam5. For the myXCLASSFit, myXCLASSMapFit, myXCLASSMapRedoFit, and LineID function the phenomenological parameters can be defined globally in a obs. xml file as well, see Sect. “Observational xml file”.

Derivations

Detection Equation

In absence of scattering, the radiative transfer equation

\[\frac{\mathrm{d} I_{\nu}}{\mathrm{d} s} = -\kappa_{\nu}(s) I_{\nu} + \epsilon_{\nu}(s)\]

describes the propagation of radiation which passes through a medium. During the propagation photons are absorbed and emitted indicated by the emission and absorption coefficients \(\epsilon_\nu\) and \(\kappa_{\nu}\), respectively.

The optical depth \(\tau_{\nu}(s)\) which measures the distance in units of the mean free path, is given by

(22)\[\tau_{\nu}(s) = \int_{s'=0}^{s'=s} \kappa_{\nu} \, d s'.\]

By using the source function \(S_{\nu} = \frac{\epsilon_{\nu}}{\kappa_{\nu}}\), the radiative transfer equation can be expressed as

(23)\[\frac{\mathrm{d} I_{\nu}}{\mathrm{d} \tau} = I_{\nu}(\tau) + S_{\nu}(\tau).\]

Within the local thermodynamic equilibrium (LTE) the source function is described by Planck’s law for black body radiation \(B_{\nu}(T)\). Integrating Eq. (23) over \(\tau\) leads to

\[I_{\nu}(s) = I_{\nu}(0) \, e^{-\tau_{\nu}} + \int_{\tau'=0}^{\tau'=\tau} S_{\nu}(\tau'_{\nu}) e^{\tau'_{\nu} - \tau_{\nu}} d \tau'_{\nu}.\]

Assuming a constant source function, i.e. constant emission and absorption coefficients through the medium, the transfer equation can be written as

(24)\[I_{\nu}(s) = I_{\nu}(0) \, e^{-\tau_{\nu}} + S_{\nu} \left(1 - e^{-\tau_{\nu}} \right).\]

Additionally, we have to take into account, that the different components may not cover the whole beam, i.e. that the background behind a certain component might contribute as well. So, we have to extend Eq. (24) by introducing the beam filling factor \(\eta\), Eq. (2), which describes the fraction of the beam covert by a component

(25)\[I_{\nu}(s) = \eta \, \left[I_{\nu}(0) \, e^{-\tau_{\nu}} + S_{\nu} \left(1 - e^{-\tau_{\nu}} \right)\right] + \left(1 - \eta\right) \, I_{\nu}(0).\]

Here, the term \(\eta \, I_{\nu}(0) \, e^{-\tau_{\nu}}\) indicates the attenuated radiation from the background \(I_{\nu}(0)\). The second term \(\eta \, S_{\nu} \left(1 - e^{-\tau_{\nu}} \right)\) describes the self attenuated radiation emitted by a certain component. Finally, the last term \(\left(1 - \eta\right) \, I_{\nu}(0)\) represents the contribution of the background which is not covert by a component.

In real observations, we do not measure absolute intensities but only differences of intensities, i.e. we have to subtract the OFF-position from Eq. (25) as well, where we have an intensity caused by the cosmic background \(J_\mathrm{CMB}\). So, we achieve

(26)\[I_{\nu}(s) = \eta \, \left[I_{\nu}(0) \, e^{-\tau_{\nu}} + S_{\nu} \left(1 - e^{-\tau_{\nu}} \right)\right] + \left(1 - \eta\right) \, I_{\nu}(0) - J_\mathrm{CMB}.\]

Beam Filling Factor

Fig. 5 a)~A single two-dimensional Gaussian function, Eq. (27), with \(\sigma_x = \sigma_y = 1.9\) and \(\mu_x = \mu_y = 0\) representing a Gaussian beam of a telescope. b)~Cut through the Gaussian function described in a) at half height. c)~A single two-dimensional Gaussian function, with \(\sigma_x = \sigma_y = 1.1\) and \(\mu_x = 0.5\) and \(\mu_y = 0.9\) representing a Gaussian beam of a point source. d)~Cut through the Gaussian function described in a) at half height. e)~Convolved Gaussian of telescope and point`source as given by Eq. (29). f)~Cut through the convolved Gaussian function described in e) at half height.

The derivation of the beam filling factor (2) starts with the normalized two-dimensional Gaussian function

(27)\[g(x, y, \sigma_x, \sigma_y, \mu_x, \mu_y) = \frac{1}{\sqrt{2 \, \pi \, \left(\sigma_x^2+\sigma_y^2\right)}} \, e^{-\left(\frac{\left(x-\mu_x\right)^2}{2 \, \sigma_x^2}+\frac{\left(y - \mu_y\right)^2}{2 \, \sigma_y^2}\right)},\]

where \(\sigma_x^2\) and \(\sigma_y^2\) describe the variances and \(\mu_x\) and \(\mu_y\) the center along the \(x\) and \(y\) axis, respectively.

Observing a Gaussian shaped extended source with a telescope is described by a convolution of two two-dimensional Gaussian functions:

(28)\[\begin{split}\left(g_1 * g_2 \right) &= \int_{-\infty}^\infty \, \int_{-\infty}^\infty \, g_1(x - u, y - v, \sigma_{x,1}, \sigma_{y,1}, \mu_{x, 1}, \mu_{y, 1}) \cdot g_2(u, v, \sigma_{x,2}, \sigma_{y,2}, \mu_{x,2}, \mu_{y, 2}) \, du \, dv\\ &= \frac{1}{2 \pi \sqrt{\left(\sigma_{x,1}^2 + \sigma_{x,2}^2\right) \left(\sigma_{y,1}^2 + \sigma_{y,2}^2\right)}} \, e^{-\left(\frac{\left(\mu_{x,1} + \mu_{x,2} - x\right)^2}{2 \, \left(\sigma_{x,1}^2 + \sigma_{x,2}^2\right)} + \frac{\left(\mu_{y,1} + \mu_{y,2} - y\right)^2}{2 \, \left(\sigma_{y,1}^2 + \sigma_{ y,2}^2\right)}\right)}.\end{split}\]

Assuming that \(g_1\) describes the telescope with \(\mu_{x,1} = \mu_{y, 1} \equiv 0\) and that telescope and extended source are described by non-elliptical Gaussians, i.e. \(\sigma_{x,1} = \sigma_{y,1} \equiv \sigma_1\) and \(\sigma_{x,2} = \sigma_{y,2} \equiv \sigma_2\), Eq. (28) can be simplified to 25

(29)\[\left(g_1 * g_2 \right) = \frac{1}{2 \pi \left(\sigma_1^2 + \sigma_2^2\right)} \, e^{-\frac{\left(x - \mu_{x,2}\right)^2 + \left(y - \mu_{y,2}\right)^2}{2 \left(\sigma_1^2 + \sigma_2^2\right)}}.\]

The FWHM of the resulting Gaussian is given by

(30)\[{\rm FWHM} = 2 \sqrt{2 \, \log 2} \sqrt{\sigma_1^2 + \sigma_2^2}\]

which describe an area of

(31)\[\begin{split}A_{\rm conv}^{\rm FWHM} &= \pi \cdot 2 \, \log 2 \cdot \left(\sigma_1^2 + \sigma_2^2\right) \\ &= \frac{\pi}{4} \cdot \left(\theta_1^2 + \theta_2^2\right).\end{split}\]

In the last line we used the relation between the variances \(\sigma_{1,2}\) and the user defined FWHM of telescope (\(\theta_1 \equiv \theta_t\)) and source size (\(\theta_2 \equiv \theta^{m,c}\)) which is given by

(32)\[\theta_i = 2 \, \sqrt{2 \, \log 2} \cdot \sigma_i.\]

The beam filling factor Eq. (2) is defined as ratios of areas

(33)\[\eta_{g_1,g_2} = \frac{A_{\rm source}^{\rm FWHM}}{A_{\rm conv}^{\rm FWHM}} = \frac{\frac{\pi \, \theta_2^2}{4}}{\frac{\pi}{4} \cdot \left(\theta_1^2 + \theta_2^2\right)} = \frac{\theta_2^2}{\left(\theta_1^2 + \theta_2^2\right)},\]

which is completely independent of the position \(\mu_x\) and \(\mu_y\) of the extended source within the telescope beam.

If more extended sources are observed with the telescope, we have to convolve the already convolved Gaussian Eq. (29) with a further two-dimensional Gaussian function with variance \(\sigma_{x,3} = \sigma_{y,3} \equiv \sigma_3\) and center \(\mu_{x,3}\) and \(\mu_{y,3}\). We get

\[\left(\left(g_1 * g_2 \right) * g_3 \right) = \left(g_1 * g_2 * g_3 \right) = \frac{1}{2 \pi \left(\sigma_1^2 + \sigma_2^2 + \sigma_3^2\right)} \, e^{\left(-\frac{\left(\mu_{x,2} + \mu_{x,3} - x\right)^2 + \left(\mu_{y,2} + \mu_{y,3} - y\right)^2}{2 \, \left(\sigma_1^2 + \sigma_2^2 + \sigma_3^2\right)}\right)}.\]

The FWHM of the resulting Gaussian is given by

(34)\[{\rm FWHM} = 2 \sqrt{2 \, \log 2} \sqrt{\sigma_1^2 + \sigma_2^2 + \sigma_3^2},\]

which describe an area of

(35)\[\begin{split}A_{\rm conv}^{\rm next} &= \pi \cdot 2 \, \log 2 \cdot \left(\sigma_1^2 + \sigma_2^2 + \sigma_3^2\right) \\ &= \frac{\pi}{4} \cdot \left(\theta_1^2 + \theta_2^2 + \theta_3^2\right).\end{split}\]

In the last line we used again the relation between the variances and FWHM, Eq. (32). Now, we can again define a beam filling factor similar to Eq. (2) and get

(36)\[\eta_{g_1,g_2,g_3} = \frac{\frac{\pi \, \theta_2^2}{4}+\frac{\pi \, \theta_3^2}{4}}{\frac{\pi}{4} \cdot \left(\theta_1^2 + \theta_2^2 + \theta_3^2\right)} = \frac{\theta_2^2+\theta_3^2}{\left(\theta_1^2 + \theta_2^2 + \theta_3^2\right)}.\]

When we observe more than three extended sources with a telescope we have to iteratively convolve the resulting Gaussian with a further two-dimensional Gaussian function and we can generalize Eq. (36) to

(37)\[\eta_{\rm general} = \frac{\sum_{i > 1} \theta_i^2}{\left(\sum_i \theta_i^2\right)} = \frac{\sum_{i > 1} \theta_i^2}{\theta_t^2 + \sum_{i > 1} \theta_i^2},\]

where \(\theta_1 \equiv \theta_t\) describes the FWHM of the Gaussian shaped beam of the telescope, defined by the diffraction limit, Eq. (3).

Optical Depth

Fig. 6 Transitions between the lower \(l\) and the upper \(u\) level with the corresponding Einstein coefficients.

In order to derive Eq. (7) we consider a system which involves radiative transitions between a lower \(l\) and an upper \(u\) level only. As shown in Fig. 6 , the lower level has an energy \(E_l\) and the upper level an energy \(E_u > E_l\). With

(38)\[h \, \nu_{u,l} = E_u - E_l\]

describing the energy difference between these two levels we can express the emissivity due to spontaneous radiative decay as

(39)\[\epsilon_{l,u,\nu} = \frac{h \, \nu}{4 \, \pi} \, n_u \, A_{u,l} \, \phi_{l,u} (\nu),\]

where \(A_{u,l}\) describes the Einstein A-coefficient, or radiative decay rate for the transition from the lower \(l\) to the upper \(u\) level. The expression \(1 / A_{u,l}\) gives the averaged time, that a quantum mechanical system can stay in level \(u\) before radiatively decaying to level \(l\), where we assume no collisional (de-)excitation. The expression \(\phi_{l,u} (\nu)\) describes the line profile of the transition of photons of frequency \(\nu\) and is normalized to 1, i.e. \(\int_0^\infty \phi(\nu) \, d \nu = 1\).

Similar to Eq. (39) we can write the extinction coefficient, which describes the radiative excitation from the lower to the upper level

(40)\[\kappa_{l,u,\nu}^{\rm ext} = \frac{h \, \nu}{4 \, \pi} \, n_l \, B_{l,u} \, \phi_{l,u} (\nu),\]

where \(B_{l,u}\) describes the Einstein B-coefficient for extinction.

In addition to spontaneous emission and extinction we have to take the stimulated emission into account, which can be described by adding a negative opacity contribution to Eq. (40):

(41)\[\kappa_{l,u,\nu} = \frac{h \, \nu}{4 \, \pi} \, \left(n_l \, B_{l,u} - n_u \, B_{u,l}\right) \, \phi_{l,u} (\nu),\]

where \(B_{u,l}\) represents the Einstein B-coefficient for stimulated emission. So, for \(n_l B_{l,u} < n_u \, B_{u,l}\) we get laser (maser) emission.

The different Einstein coefficients are related to each other by the Einstein relations:

(42)\[\begin{split}A_{u,l} &= \frac{2 \, h \, \nu^3}{c_{\rm light}^2} \, B_{u,l},\\ g_l \, B_{l,u} &= g_u \, B_{u,l}.\end{split}\]

Following Eq. (22), the differential optical depth \(\tau_{\nu}\) is defined as

(43)\[\begin{split}d \tau_{\nu} = \kappa_{\nu} \, ds &= \left(\frac{h \, \nu}{4 \, \pi} \, \left(n_l \, B_{l,u} - n_u \, B_{u,l}\right) \, \phi_{l,u} (\nu) \right) \, ds \\ &= \left(\frac{c_{\rm light}^2}{8 \, \pi \, \nu^2} \, A_{u,l} \, \left(n_l \, \frac{g_u}{g_l} - n_u \right) \, \phi_{l,u} (\nu) \right) \, ds,\end{split}\]

where we used the Einstein relations Eqn. (42) in the last line. By assuming LTE conditions and therefore Boltzmann population distribution

\[\frac{n_u}{n_l} = \frac{g_u}{g_l} \, e^{(-E_u - E_l)/k_B T_{\rm ex}} = \frac{g_u}{g_l} \, e^{-h \, \nu_{u,l}/k_B T_{\rm ex}}\]

we can rewrite Eq. (43) by using Eq. (38)

(44)\[d \tau_{\nu} = \left(\frac{c_{\rm light}^2}{8 \, \pi \, \nu^2} \, A_{u,l} \, n_u \, \left(e^{h \, \nu_{u,l} / k_B \, T_{\rm ex}} - 1 \right) \, \phi_{l,u} (\nu) \right) \, ds.\]

Finally, we have to integrate along the line of sight and obtain the optical depth \(\tau_{\nu}\)

(45)\[\tau_{\nu} = \frac{c_{\rm light}^2}{8 \, \pi \, \nu^2} \, A_{u,l} \, N_u \, \left(e^{h \, \nu_{u,l} / k_B \, T} - 1 \right) \, \phi_{l,u} (\nu),\]

where \(N_u = \int n_u \, ds\) describes the column density of a certain molecule in the upper state.

In order to express Eq. (45) in terms of the total column density \(N_{\rm tot} = \sum_{i=0}^\infty n_i\) we start again with the Boltzmann population distribution

(46)\[N_{\rm tot} = \sum_{i=0}^\infty n_i = n_j \, \sum_{i=0}^\infty \frac{n_i}{n_j} = n_j \, \sum_{i=0}^\infty \frac{g_i}{g_j} \, e^{(-E_i+E_j)/k_B T_{\rm ex}} = \frac{n_j}{g_j} \, e^{E_j/k_B T_{\rm ex}} \cdot Q \left(T_{\rm ex} \right),\]

where we used the partition function \(Q(T_{\rm ex})\), which is defined as sum over all states

(47)\[Q(T_{\rm ex}) = \sum_{i=0}^\infty g_i \, e^{-E_i/k_B \, T_{\rm ex}}.\]

For our purpose we rewrite Eq. (46) in terms of \(N_u\) and \(N_l\):

(48)\[N_{\rm tot} = \frac{N_u}{g_u} \, e^{E_u/k_B T_{\rm ex}} \cdot Q \left(T_{\rm ex} \right).\]

Inserting this expression into Eq. (45) gives

(49)\[\begin{split}\tau_{\nu} &= \frac{c_{\rm light}^2}{8 \, \pi \, \nu^2} \, A_{u,l} \, \left[\frac{N_{\rm tot} \, g_u}{Q \left(T_{\rm ex} \right)} \, e^{-E_u/k_B T_{\rm ex}} \right] \, \left(e^{h \, \nu_{u,l} / k_B \, T} - 1 \right) \, \phi_{l,u} (\nu) \\ &= \frac{c_{\rm light}^2}{8 \, \pi \, \nu^2} \, A_{u,l} \, N_{\rm tot} \, \frac{g_u \, e^{-E_l/k_B T_{\rm ex}}}{Q \left(T_{\rm ex} \right)} \, \left(1 - e^{-h \, \nu_{u,l} / k_B \, T} \right) \, \phi_{l,u} (\nu),\end{split}\]

where we used Eq. (38) to achieve Eq. (7) for a single line of a certain component and molecule.

Elliptical rotated Gaussians

An arbitrary rotated elliptical 2D-Gaussian can be written as

(50)\[f(x,y) = \frac{1}{V} \cdot A \cdot \exp \left[ -\left( a(x - x_0)^2 + 2 b (x - x_0) (y - y_0) + c (y - y_0)^2 \right) \right],\]

where the coefficient \(A\) describes the amplitude, \(x_0\), \(y_0\) the center and \(\sigma_x\), \(\sigma_y\) are the \(x\) and \(y\) spreads of the blob. Additionally, \(V\) indicates the volume of the Gaussian and is given by

(51)\[V = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y) \, dx \, dy = 2 \pi A \sigma_{x} \sigma_{y}.\]

The coefficients \(a\), \(b\), and \(c\) are

\[\begin{split}a &= {\frac {\cos^{2} \theta} {2 \sigma_{x}^{2}}} + {\frac{\sin^{2} \theta} {2 \sigma_{y}^{2}}}, \\ b & = -{\frac{\sin 2\theta } {4 \sigma_{x}^{2}}} + {\frac{\sin 2 \theta } {4 \sigma_{y}^{2}}}, \\ c &= {\frac{\sin^{2} \theta }{2 \sigma_{x}^{2}}} + {\frac{\cos^{2} \theta }{2\sigma_{y}^{2}}} .\end{split}\]

where the blob is rotated by a clockwise angle \(\theta\).

Footnotes

3: A derivation of the expression can be found in the Sect. “Detection Equation”.
4: Derivations of the beam filling factor Eq. (2), are described in the Sect. “Beam Filling Factor”.
5: A derivation of the expression is given in the Sect. “Optical Depth.
6: The indices \(u\) and \(l\) represent upper and lower state of transition \(t\), respectively.
7: Because the database usually does not describe the partition functions at the given excitation temperature \(T_{\rm ex}^{m,c}\), the value of \(Q \left(m,T_{\rm ex}^{m,c} \right)\) is computed from a linear interpolation in log space. (With the new catalog, extrapolation should not be necessary for most conditions encountered in molecular cores.)
8: For molecules and RRLs the name declared in the molfit file must be identical to the name of the molecule included in the database.
9: Here, \(i = 1\) indicates the core component, i.e. components which are closest to the background.
10: A derivation of the expression can be found in the Sect. “Detection Equation”.
11: Note, the definition of the source size \(\theta^{m,c}\) for an foreground component will be ignored, because we assume that all foreground layers cover the whole beam.
12: For single dish observations the telescope beam FWHM size is related to the diameter of the telescope by the diffraction limit. For the sub-beam description, Eq. (3) is used to compute the beam size for single dish observations for the lowest frequency of each frequency range.This value is used for the whole range, respectively.
13: https://www.fftw.org/
25: Here “*” indicates the convolution of two functions \(g_1\) and \(g_2\).

References

Bresenham_1977: Bresenham, J. 1977, Commun. ACM, 20, 100–106
Cesaroni_Walmsley_1991: Cesaroni, R. & Walmsley, C. 1991, Astronomy and Astrophysics, 241, 537
Hildebrand_1983: Hildebrand, R. 1983, Molster, FJ, Dorschner, J., Henning, T., & Mutschke, H.
Ossenkopf_Henning_1994: Ossenkopf, V. & Henning, T. 1994, Astronomy and Astrophysics, 291, 943
Padua_2011: Padua, D. 2011, Encyclopedia of parallel computing (Springer Science & Business Media)
Stahler_Palla_2005: Stahler, S. & Palla, F. 2005, Stahler, SW & Palla, F, 3, 9