© A. Roueff et al. 2020

1 Introduction

Spectroscopic measurements are commonly used to probe astrophysical objects. In the interstellar medium, the moderate temperatures and densities of diffuse and molecular clouds (T_kin ~ 10−100 K, and n ~ 10²−10⁵ cm⁻³, Draine 2011) are well-suited for the emission in the low-energy rotational lines of molecules such as carbon monoxide, which are accessible at millimeter wavelengths. The advent of sensitive broadband heterodyne receivers provides homogeneous data sets of various CO isotopologues and other species with high signal-to-noise ratios (S/N) over large fields of view. The Outstanding Radio-Imaging of OrioN B (ORION-B IRAM; co-PIs: J. Pety and M. Gerin) 30-m large program is aimed at imaging five square degrees towards the southern part of the Orion B molecular cloud over most of the 3 mm atmospheric window. Carbon monoxide is of particular interest because it is one of the most abundant molecules after molecular hydrogen. Using the unsupervised meanshift clustering method on the intensities of the CO isotopologues, Bron et al. (2018) have shown that it is possible to cluster the emission line data across the analyzed field of view into a few classes of increasing (column) densities. In two empirical studies, Gratier et al. (2017, 2020) show that both qualitatively and quantitatively the ¹²CO(1− 0), ¹³CO(1− 0), and C¹⁸O (1−0) lines are indeed tracing the molecular gas well. Their quantitative comparisons show that the H₂ column density deduced from the dust emission can be accurately estimated from the ¹²CO(1− 0), ¹³CO(1− 0), and C¹⁸O (1−0) lines in the column density range from 10²¹ to ≳ 10²² cm⁻².

Assuming identical excitation temperatures, the opacity of the ground state transitions is expected to be smaller for ¹³CO than for ¹²CO, and even smaller for C¹⁸O, because of the difference in elemental abundances, ¹²C∕¹³C ~ 60 and ¹⁶O ∕¹⁸O ~ 500 (Langer & Penzias 1990; Wilson & Rood 1994). These three lines can thus be used to probe progressively higher gas column densities, provided the relative elemental abundances are constant and the CO isotopologue abundances track the elemental abundances. However, chemical models and observations show that selective photodissociation and carbon isotopic fractionation can significantly modify the relative abundances of carbon monoxide isotopologues, as compared to elemental abundances (Visser et al. 2009; Liszt 2017; Roueff et al. 2015). Fractionation via the exchange reaction between ¹³C⁺ and ¹²CO leads to an enhancement of the ¹³CO abundance in the diffuse or translucent regions where CO and C⁺ coexist and the kinetic temperature remains moderate (≲50 K, Liszt & Pety 2012). This mechanism widens the ¹³CO emitting region and brings it closer to that of ¹²CO, which favors the simultaneous detection of both isotopologues across wide fields of view. However, the ratio of isotopologue abundances can be significantly different from the ratio of elemental abundances, which complicates the determination of the ¹³C elemental abundance from CO observations only. Up to now, the most reliable determinations of the ¹²C∕¹³C elemental abundance ratio have been obtained using C⁺ or C observations in regions without significant fractionation (e.g., Keene et al. 1998; Ossenkopf et al. 2013) or they have involved C¹⁸O and the doubly isotopic species ¹³C¹⁸O (Langer & Penzias 1990).

No such fractionation reaction exists for oxygen. However, the more abundant CO isotopologues shield themselves from the destructive effect of UV photons more efficiently than less abundant isotopologues because the photodissociation of carbon monoxide is governed by line absorption. This effect, called selective photodissociation, plays an important role here. It has been studied in detail through laboratory experiments (e.g., Stark et al. 2014) and in models of photo-dissociation regions (e.g., Visser et al. 2009). In observations, it is clearly seen as an offset between the threshold for the apparition of ¹²CO (near A_V = 0.5 mag) and C¹⁸O (1.5 mag) in the Taurus molecular cloud and this offset is not due to a difference in the detection sensitivity (Frerking et al. 1989; Cernicharo & Guelin 1987). Typically, the ¹³CO abundance is enhanced through fractionation in the same regions where the C¹⁸O abundance decreases due to selective photodissociation. This leads to a broad range of the ¹³CO/C¹⁸O abundance ratio for a given set of elemental abundances.

Determining the ratio of elemental abundances of the C and O isotopes is interesting because it provides information on the stellar populations which have produced these elements. Some external galaxies exhibit CO isotopologue ratios that significantly differ from the expected value based on the mean elemental abundances in the solar neighborhood. Such differences can trace differences in elemental abundances, hence in stellar populations and IMF shape (Sliwa et al. 2017; Martín et al. 2019). However, a proper account of isotopic chemistry described above must be performed in order to use the information on the relative abundances of the CO isotopologues.

Finally, Orkisz et al. (2019) have shown in their analysis of the filamentary structure of the Orion B molecular cloud that the gas velocity dispersion determined from C¹⁸O reaches a minimum value in the filament ridges and that it is always lower than the velocity dispersion determined by ¹³CO. This suggests that this variation of velocity dispersion between CO isotopologue traces the dissipation of turbulence when entering the dense filaments inside molecular clouds.

Constraining all these astrophysical effects relies on a precise derivation of physical conditions and chemical composition from spectroscopic observations. This, in turn, relies on the resolution of the radiative transfer equation because the line intensities and profiles bear information on the line emission mechanisms. The large data volumes provided by observational programs like ORION-B require new statistical analysis methods using the information in an optimal way and a derivation of the physical parameters and their associated errors with a rigorous methodology. For instance, the emission of the lowest rotational transitions of the three major isotopologues of carbon monoxide, ¹²CO, ¹³CO, and C¹⁸O, is commonly used to determine the molecular gas column density and evaluate the mass of molecular gas. Because these lines can now be observed simultaneously, leading to an homogeneous flux calibration and therefore precise relative calibration, it is essential to have a good estimation of the precision on the mass estimate.

In estimation theory, the Cramer Rao bound (CRB) provides a precision of reference that does not depend on a specific estimator of the searched quantity, but only on the physical model and the statistical properties of the noise (see, e.g., Bonaca & Hogg 2018; Espinosa et al. 2018). The CRB further allows for the quantification of the loss of precision due to degeneracies between the estimated parameters (for instance, column density and excitation temperature). Hence, a large value of this bound indicates insufficient data or knowledge with respect to a given physical model. We apply this technique here in the simplest possible model framework, that is, the emission of lines in the local thermodynamic equilibrium (LTE), which can be fully expressed using analytical equations.

The level populations of interstellar molecules result from the balance of collisional (and possibly radiative) excitation and radiative and collisional de-excitation. Therefore the level populations often deviate from LTE conditions because the collisions are not efficient enough to populate all energy levels according to a Boltzmann distribution. With its low dipole moment (0.1 Debye) and high abundance relative to H₂, the low energy rotational lines of carbon monoxide are bright and easily thermalized in collisions with H₂, H, and He. This means that the LTE model is still a good approximation for this molecule; in other words, the rotational level populations can be described by a Boltzman distribution at a single excitation temperature (Liszt 2006; Leung & Liszt 1976; Goldsmith & Langer 1999; Goldreich & Kwan 1974). Deviations from the LTE model have been studied theoretically. For instance, using non-local, non-LTE radiative transfer models of a uniform (constant density and temperature) spherical cloud, Bernes (1979) has shown that the excitation temperatures of the ¹²CO(1− 0) and ¹²CO(2− 1) lines exhibit moderate spatial variations from edge to center. It has been concluded that the LTE model is mostly valid for the ground state transition and deviations from this approximation increase with the quantum number of the upper level (van der Tak et al. 2007).

With a wide range of physical conditions, from bright far-UV illuminated regions to cold and shielded regions through diffuse and translucent gas irradiated by a moderate radiation field, the Orion B molecular cloud is an ideal place for probing the extent to which fractionation and selective photodissociation can modify the elemental abundance ratio. It is also a good region for probing the differences in excitation between isotopologues as the simple hypothesis of equal excitation temperatures for ¹²CO, ¹³CO, and C¹⁸O may not be valid, as discussed in Bron et al. (2018).

The article is organized as follows. Section 2 presents the data used in this paper. Section 3 summarizes the mathematical formulation of the LTE radiative transfer. Section 4 presents the computation and analysis of the precision achievable for this theoretical framework. Section 5 illustrates the proposed methodology on actual data sets. Section 6 focuses on the astrophysical interpretations of these results. Appendix A details the calculation of some gradients necessary to compute the Fisher matrix. Appendix B describes our implementation of the maximum likelihood estimator and Appendix C presents a discussion of the performance of this estimator.

2 Description of the data

We attempted to estimate the velocity field, the column density, and the excitation temperature of the CO isotopologues from the analysis of the ¹³CO(1−0), ¹³CO(2− 1), C¹⁸O (1−0), C¹⁸O (2−1) and ¹²CO(1− 0) lines towards parts of the Orion B molecular cloud. We compared our results with the dust-traced H₂ column density and dust temperature. This section describes the associated data sets.

2.1 IRAM-30 m observations

2.1.1 3 mm CO lines from the ORION-B large program

The 3 mm data were obtained with the IRAM-30 m as part of the ORION-B large program. Pety et al. (2017) present in detail the acquisition and reduction of the dataset used in this study. In short, the applied data were acquired at the IRAM-30m telescope using the EMIR receiver and Fourier transform spectrometer from August 2013 to November 2014. The frequency range from 84 to 116 GHz was completely sampled at 200 kHz spectral resolution. The J = 1−0 lines of the CO isotopologues analyzed here are observed in a single receiver tuning. These lines are thus well inter-calibrated. The absolute flux calibration at 3 mm for the IRAM-30 m telescope is estimated to be better than 5%.

2.1.2 1 mm CO lines

The ¹³CO(2−1) and C¹⁸O (2−1) lines were also observed at the IRAM-30 m in 2006 (PI: N. Peretto) using the ABCD generation of receivers and the VESPA auto-correlator. The two lines were observed simultaneously ensuring an excellent inter-calibration.

Data reduction was carried out using the GILDAS¹/CLASS software. The contribution of the atmosphere was first removed (ON–OFF procedure) and the data were calibrated to the scale using the standard chopper-wheel method (Penzias & Burrus 1973). The data were then converted to main-beam temperatures using the standard forward (0.94) and main-beam (0.62) efficiencies for the ABCD receiver around 220 GHz². The resulting absolute flux calibration is estimated to be better than 10%. We subtracted a first-order baseline from every spectrum, excluding the velocity range from 5 to 15 km s⁻¹ in the local standard of rest (LSR) frame. Finally, the spectra were gridded into a data cube through a convolution with a Gaussian kernel of full width at half maximum (FWHM) at ~ 1∕3 of the IRAM-30 m telescope beamwidth at the line rest frequency.

2.2 Herschel observations

In order to get independent constraints on the physical conditions in the Orion B cloud, we used the dust continuum observations from the Herschel Gould Belt Survey (André et al. 2010; Schneider et al. 2013) and from the Planck satellite (Planck Collaboration I 2011). The fit of the spectral energy distribution by Lombardi et al. (2014) gives us access to the spatial distributions of the dust opacity at 850 μm and of the dust temperature. As in Pety et al. (2017), we converted τ_{850 μm} to visual extinctions using A_V = 2.7 × 10⁴ τ₈₅₀ mag, and the visual extinction into H₂ column density using N(H₂)∕A_V = 0.9 × 10²¹H cm⁻² mag⁻¹.

2.3 Field of view

We aimed to jointly analyze the J = 1−0 and J = 2−1 lines of the CO isotopologues. We thus restricted the field of view to the region that was observed at 3 and 1 mm. This covers 19′ × 26′ towards the Orion B molecular cloud part that contains the Horsehead nebula, and the HII regions NGC 2023 and IC 434. The cubes used here are rotated counterclockwise by 14° around the projection center (05^h40^m54.270^s, −02°28′00.00″) in the RA/Dec J2000 reference frame (see Fig. 1). The coordinates are given in offsets (δx, δy) in arcseconds from this projection center. The IRAM-30 m angular resolution ranges from 11.5″ at 220 GHz to 23.5″ at 110 GHz. The position-position-velocity cubes of each line were smoothed to a common angular resolution of 23.5″ to avoid resolution effects during the comparison. At a distance of 400 pc (Menten et al. 2007), the sampled linear scales range from ~0.045 pc to ~ 3 pc.

The spectral and spatial axes were resampled in order to share the same spatial grid and velocity axis for all lines. The spectroscopic observations thus provide position-position-velocity cubes³ of 129 × 170 × 80 pixels, each pixel covering 9″ × 9″ × 0.5 km s⁻¹ (Nyquist sampling at 3 mm). Figure 1 shows the maps of the intensity integrated between 0 and 20 km s⁻¹ for the five lines of interest.

Fig. 1 Spatial distributions of the integrated intensity (in K km s−1) of the considered lines. The maps have been rotated counterclockwise by 14 degrees from the RA/Dec J2000 reference frame. The spatial offsets are given in arcsecond from the projection center located at 05h 40m 54.270s, −02°28′00.00′′. Red crosses stand for two particular lines of sight, which are analyzed in Fig. 10.

2.4 Noise

In this paper, the standard deviation of the noise σ_b is estimated only on negative values of each spectrum using: (1)

where T is the intensity in Kelvin, and K_neg the number of channels that have a negative value of the intensity. This allows us to compute it without a priori information on the velocity range where the line appears, but assuming, however, that the baselining removes any intensity offset. Table 1 lists the median noise estimated after spectral resampling and angular smoothing.

2.5 Line profiles

A fraction of the studied field of view shows spectra that can only be modeled with more than one velocity component along the line of sight. While we tend to adapt our formalism to handle such cases, it is not obvious that a robust statistical test should be devised to deduce the optimal number of components that must be used. This is particularly true at transitions between regions where the number of required velocity components changes in order to get a good fit. To address this issue, we used the ROHSA (Marchal et al. 2019) algorithm that makes a Gaussian decomposition based on a multi-resolution process from coarse to fine grid. Here we only used the spectra denoised by ROHSA to provide a spatially coherent estimation of the number of components and some initial estimation of their associated central velocities for each pixel.

3 Radiative transfer in local thermodynamic equilibrium

The molecular line emission and absorption in the case of local thermodynamic equilibrium (LTE) are well-known (see, e.g., Mangum & Shirley 2015). In this section, we mainly summarize the associated notations and equations so that we can more easily explain the precision analysis framework on this case in the next section. For the sake of simplicity, we focus on a single chemical species and a single velocity component along one line of sight. The observed spectrum as a function of frequency ν is defined as: (2)

where b is a (thermal) Gaussian noise and s is the spectrum associated to the species of interest. The specific intensity, s, and the associated measurement noise, b, are expressed in Kelvins, following the standard convention in radioastronomy. The data reduction (atmospheric ON-OFF calibration and spectrum baselining to subtract the slowly-varying continuum residual from the receiver and the atmosphere) delivers a noise b that is centered (i.e., with zero-mean) and whose variance can be considered constant over each line profile.

We assume that two lines (l ∈{1, 2}) from the same species are being observed. The photons of each line are emitted at the rest frequency of the line, ν_l, and redshifted in frequency because of the Doppler shift due to the motion of the gas along the line of sight in the observation frame, typically the local standard of rest (LSR). The photon is thus received at the redshifted frequency , where V is the velocity of the emitting cell of gas in the LSR frame and c is the speed of light. This equation is the radio low-velocity approximation of the Doppler effect. The Doppler effect due to the motion of the observer relative to the LSR is automatically taken into account in the data acquisition process. Therefore, each line of the dataset is analyzed in the LSR frame. In this frame each line is centered around a typical velocity, noted Δ_V. This velocity is related to the redshifted centroid frequency of the line, , through a particular case of the previous equation, (3)

We assume that the only background source of emission is the cosmic microwave background (CMB). In this case, the intensity s at observed frequency ν around can be written as (4)

where T_CMB is the known CMB temperature (T_CMB = 2.73 K, Mather et al. 1994), T_ex is the unknown excitation temperaturealong the line of sight and J is a measureof intensity at a given temperature, (5)

where B(T, ν) is the spectral distribution of the radiation of a black body at temperature T. The term in Eq. (4) represents the emission or absorption by the emitting or absorbing medium along the line of sight, considered as a uniform slab. The function Ψ is the profile that corresponds to the integrated opacity through the whole slab. For each line l, it can be written as: (6)

In this equation, σ_V is the velocity dispersion of the source along the line of sight. It varies as a function of the local physical conditions (higher temperatures and higher turbulence lead to larger values). The function ϕ is a Gaussian profile: (7)

where σ_ν is the frequency dispersion in the source rest frame. It is related to σ_V by σ_ν = ν_l σ_V∕c, because of the Doppler effect. Finally, the amplitude α_l associated to the Gaussian profile ϕ and line l is: (8)

where A_l is the Einstein spontaneous emission rate for line l, g_up is the degeneracy of the upper level of the line, E_up is its energy (in units of Kelvin), and N is the column density of the species along the line of sight. The partition function Q(T_ex) is tabulated in molecular databases (e.g., CDMS, Müller et al. 2001, or JPL, Pickett et al. 1998), and its temperature dependence can be interpolated for each species. The partition function is computed as the sum of the populations of all energy levels E_k. If the energy levels are expressed in Kelvin, Q(T_ex) can be written as: (9)

The parameter α_l is related to the line opacity τ_l (10)

which is dimensionless. The excitation temperature is defined from the ratio of the population in the upper (n_up) and lower (n_low) levels of the studied line, (11)

When the molecules are in thermal equilibrium with their environment, the temperature T_ex is equal to the gas kinetic temperature. The kinetic temperature is not known and must be estimated.

In the previous equations, the physical characteristics of the gas (N, T_ex, σ_V, and Δ_V) depend on the specific line of sight on the sky. Moreover, while observers try to get a uniform noise when observing the source, this is never perfect and it is important to assume that the noise standard deviation σ_b also depends on the specific line of sight on the sky. These considerations imply that α_l, τ_l, , Ψ_l, s, and x also depend on the sky position.

4 Cramer-Rao bound analysis

In this paper, we aim to estimate the physical parameters of the LTE model presented in Sect. 3 based on the ORION-B data (see Sect. 2). Even when the LTE model is perfectly verified, the presence of an additive Gaussian noise induces some uncertainty on the estimation. For each physical parameter θ estimated as , the estimation error can be quantified with the mean square error (MSE) , where ⟨.⟩ represents the statistical mean over the different realizations of the noise b. MSE can be estimated with Monte Carlo simulations, but the result then depends on the choice of the implemented estimator. For example, when the MSE is large, it is not clear whether it is due to the choice of the estimator or to a lack of information within the data. In estimation theory, the Fisher matrix allows to quantify the amount of information in the considered problem. It provides a reference precision, named Cramer-Rao bound (CRB), which does not depend on the choice of a specific estimator of the searched quantity, but only on the physical model and the statistical properties of the noise (see, e.g., Bonaca & Hogg 2018; Espinosa et al. 2018).

Mathematically speaking, the CRB noted is simply a lower bound on the MSE of unbiased estimators (see Eq. (15)). Indeed, the MSE is equal to the estimation variance for an unbiased estimator because the estimation MSE is in general equal to the sum of its variance and its bias squared, that is, (12)

Therefore, a high CRB value implies that any unbiased estimator will necessarily have a high dispersion around the true value θ. A high CRBcan be understood as a lack of information provided in the underlying model with respect to the considered level of noise. Whenit occurs, one solution can be to introduce additional a priori knowledge or to make another measurement with a better S/N. In contrast, a low CRB value does not necessarily imply that there exists an unbiased estimator with a low dispersion around the true value θ. CRB is only a lower bound and it can be overly optimistic. It is thus necessary to build an estimator, which can be tested by comparing its variance with the CRB. If the estimator is unbiased and its variance is equal to the CRB, then one knows that there does not exist any better unbiased estimator. In this section, we analyze the CRB (i.e., a bound on the variance of all unbiased estimators) and in Sect. 5.1 we check with Monte Carlo simulations that an efficient estimator (i.e., one whose variance reaches the CRB) exists.

Here we compute the Fisher matrix and the associated CRB precisions for the LTE radiative transfer. We then study the variations of these reference precisions for the different unknown physical parameters (Δ_V, σ_V, T_ex, and N) as a function of the excitation temperature and column density. We finally use the CRB precision to answer two questions. The first question relates to what is the maximum noise tolerable to get a given relative precision on these parameters. Here we compare the cases where only one (¹³CO(1−0)) or two (¹³CO(1−0) and ¹³CO(2− 1)) lines are available. Second, we consider which ¹²CO line should be observed to improve the precision achieved when only ¹²CO(1− 0) observations are available.

Fig. 2 Variations of the square root of the Cramer-Rao bound (CRB) of the centroid velocity (ΔV, left panel) and velocity dispersion (σV, middle panel) in km s−1 as afunction of the column density and the excitation temperature. Right panel: variations of a function of the product of the number of channels (K) and the signal-to-noise ratio . In all three cases, the data are simulated assuming that 13CO(1−0) and 13CO(2− 1) are measured and the unit of the image contours are km s−1. In this simulation, σb = 100 mK, ΔV = 1.1 km s−1 and σV = 0.61 km s−1.

4.1 Computing the CRB from the Fisher matrix for a single line and a single velocity component

For a given line l, a sampled version of Eq. (2) can be written over K discrete frequency channels as: (13)

When b is a centered white Gaussian noise of standard deviation, σ_b,l, and the physical model, s, is expressed as a function of a set of unknown parameters, (θ_i), the Fisher matrix, I_F, which represents the amount of information provided by the line, l, can simply be computed as (Stoica & Moses 2005): (14)

where stands for the term (i, j) of the matrix A.

In our case, the physical model for s is the LTE radiative transfer introduced in Sect. 3 and the vector of unknown parameters is , which we also write to simplify the expression of the Fisher matrix in Eq. (14). In this vector of parameters, we chose to analyze the precision of the logarithm⁴ of the column density, N, instead of directly analyzing the precision of N because the column density can vary over orders of magnitudes in giant molecular clouds.

It can be shown (Garthwaite et al. 1995) that the variance of any unbiased estimator (here , , or ) is bounded by: (15)

Each diagonal term of the inverse of the Fisher matrix is called the Cramer Rao bound of the corresponding parameter. We note them as , , and . These CRBs do not depend on the choice of theestimation algorithm and they are usually asymptotically reached by the maximum likelihood estimator (Garthwaite et al. 1995). Thus, the CRB can be considered to be providing reference precision of the estimation problem. The calculation of the gradients is detailed in Appendix A.

4.2 Generalization to two lines and two velocity components

We use the CRB analysis on the case where we observe two different lines (l ∈{1, 2}) of the same species. We assume that these lines are well separated in frequency so that their frequency supports are disjointed: (16)

The Fisher matrix of the set is simply the sum of the Fisher matrices of each transition because we assume that the unknown parameters θ are identical for the two lines.

We also use the CRB framework in the case where each observed line is emitted from two independent velocity components, that is, from two gas components characterized by different values of the unknown parameters (θ_m with m ∈{1, 2}). Equation (16) that encodesthe spectrum for line l can then be written as: (17) (18)

This means that the composite line profile is considered to be the simple sum of two velocity components that do not radiatively interact. This assumption is only correct if the two velocity components are adequately separated in velocity. This case with two components is the most complex model we study in this paper. In this case, the number of unknown parameters is eight (instead of four) and, thus, the size of the Fisher Matrix is 8 × 8 (instead of 4 × 4).

4.3 CRB variations as a function of T_ex and N

As the inversion of the Fisher matrix is done numerically, we do not have a simple explicit expression of the CRBs. In this section, we empirically analyze, thus, their evolution as a function of the physical properties of the analyzed medium in a particular case taken from the ORION-B project. To generate Figs. 2–4, we assume that the two measured lines are ¹³CO(1− 0) and ¹³CO(2− 1). The two corresponding opacities are noted τ₁ and τ₂. The number of samples is K = 80 for each line at a spectral resolution of 0.5 km s⁻¹. Only one velocity component is assumed in the remainder of this section. The amount of noise is fixed and identical for both lines at σ_b,1 = σ_b,2 = 100 mK.

The values of the Cramer-Rao bounds of the unknown parameters (i.e., T_ex, log N, Δ_V and σ_V) are then computed as a function of the values of T_ex and N. The Fisher matrices are computed following Eq. (14), and then numerically inverted to obtain . The excitation temperature T_ex is sampled logarithmically between 3 and 99 K, the column density N is sampled logarithmically between 10¹³ cm⁻² and 10¹⁹ cm⁻², and the other two parameters are kept constant at Δ_V = 1.1 km s⁻¹ and σ_V = 0.61 km s⁻¹ (arbitrarily chosen). This leads us to Figs. 2–4, where T_ex varies horizontally and N varies vertically. In these Figures, the variations of are shown instead of the variations of because the square root of the CRB is homogeneous to the estimation standard deviation. It can thus be interpreted as errorbars on the estimated parameter, θ_i. Varying T_ex and N changes not only the S/N, but also the amount of information measured by the Fisher matrix because of the nonlinearity in the radiative transfer equation.

Fig. 3 Top: variations of the square root of the CRB of Tex in Kelvin as a function of the column density and the excitation temperature. Bottom: functions of the line opacities. Left: only 13CO (1−0) is analyzed. Middle: only 13CO(2−1) is analyzed. Right: both 13CO(1−0) and 13CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to Tex and N values whichlead to singular Fisher matrices. For this analysis, σb = 100 mK, ΔV = 1.1 km s−1 and σV = 0.61 km s−1. The constant σV,0 = 1 km s−1 is introduced to have expressions in (d–f) that depend on σV, but remain homogeneous to a temperature.

4.4 Precision of the estimation of the centroid velocity Δ_V and the associated velocity dispersion σ_V

Figures 2a, b shows variations of and . For N ≥ 10¹⁶ cm⁻² and T_ex ≥ 12 K, the square root of both CRBs are smaller than 0.01 km s⁻¹. This means that any efficient unbiased estimator is expected to have a small dispersion around the actual values. This can be written and .

Figure 2c shows a function of , where is the signal-to-noise ratio defined by: (19)

This expression is used in signal processing to quantify the S/N on the “energy” of the signal. In our case, we empirically find, as a rule of thumb, (20)

While the dependency on σ_V is not presented in Fig. 2, we checked that Eq. (20) remains valid when σ_V = 0.3, 1.31, and 2 km s⁻¹. The estimation precision on Δ_V and σ_V depends on the S/N andon σ_V. This is expected because of the similarity with the problem of delay estimation in radar for which le Chevalier (1989) obtained an analytic formulation similar to Eq. (20).

4.5 Precision of the estimation of the excitation temperature T_ex

In this section, we begin to quantitatively evaluate the gain in precision when two lines are observed instead of a single one. Figure 3 compares the variations of when only ¹³CO(1−0) or ¹³CO(2− 1) is available to constrain the excitation temperature, and when both lines are available. To interpret this figure, we first remark that for low column densities, the uncertainty quickly increases leading to large values of the CRB, especially for N < 10¹⁶ cm⁻². While the variation of the CRB as a function of the excitation temperature for a given column density is monotonic in the considered range for the ¹³CO(1− 0) line, it shows a different behaviour for the ¹³CO(2−1) line, with a minimum near 6 K, and an increase of the CRB for lower values of the excitation temperature. This different behaviour is related to higher energy of the upper state of the 2− 1 transition. The emerging ¹³CO(2−1) signal, which is proportional to the population of the upper level of the transition, approaches zero and becomes close to the noise level.

For the considered example, the analysis of a single line (see Figs. 3a,b) gives a reference precision on T_ex of for typically N > [10¹⁶−10^17.5] cm⁻². The dependence on T_ex is such that the same CRB is also reached at higher column densities for higher values of T_ex. This behavior of the CRB can be qualitatively understood as resulting from the increase of the line opacity. When the opacity becomes larger than about 3, the peak temperature only depends on T_ex as the factor in Eq. (4) approaches unity. The CRB almost linearly depends on log T_ex above 6 K. The analysis of two lines (see Fig. 3c) greatly improves the situation. One reaches the same precision on T_ex for column densities that are between one and two orders of magnitude lower, that is, for typically N > [10¹⁴−10^16.5] cm⁻². Here, once again the precision almost linearly depends on logT_ex above 6 K.

The second row of Fig. 3 shows functions of the opacities. When trying for several values of σ_V, we empirically obtain: (21)

when a single line is available, either ¹³CO(1−0) or ¹³CO(2− 1). In this equation, σ_V,0 = 1 km s⁻¹ is a constant fixed so that the expression depends on σ_V, but remains homogeneous to a temperature. When these two lines are available, we obtain: (22)

Hence, according to Eqs. (21) and (22), when opacities are close to one, the gain in precision (in standard deviation) is around 20 when we observe two lines of the same species instead of a single one. These relations suggest that the parameters that control the difficulty of the estimation problem are the amount of noise σ_b, the velocity dispersions σ_V, and the opacities.

Fig. 4 Top row: variations of the square root of the CRB of logN as a function of the column density and the excitation temperature. Second and third row: correlation coefficients between efficient estimators of (Tex, logN), and (logN, σV) in the second and third rows. respectively (defined in Eqs. (23) and (24)). Bottom row: variations of functions of the opacities. Left: only 13CO(1− 0) is analyzed. Middle: only 13CO(2−1) is analyzed. Right: both 13CO(1−0) and 13CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to Tex and N values whichlead to singular Fisher matrices. For this analysis, σb = 100 mK, ΔV = 1.1 km s−1 and σV = 0.61 km s−1.

4.6 Precision of the estimation of the column density N

The top row of Fig. 4 shows that the precision on the estimation of N has a complex behavior when only one line is available. To interpret this, we note that even efficient estimators of θ have correlated components described by the correlation coefficients of the Fisher matrix. The correlation coefficient between T_ex estimations and logN estimations is given by: (23)

where is the non-diagonal element of the inverse Fisher matrix, see Eq. (14). We also introduce the correlation coefficient between log N estimations and σ_V estimations (24)

where ; see Eq. (14).

Correlation coefficients are built such that their value ranges from − 1 to 1. As long as |γ| < 1, CRBs remain finite and thus estimating parameters usually remains possible. There is a complete ambiguity between estimations of the pairs (T_ex and logN) or (log N and σ_V), only when values of |γ| = 1. In this case, the variance of these estimations becomes infinite. Figure 5 shows a simulated example where log (N) and σ_V can be accurately estimated even though they are highly (but not completely) anti-correlated. Starting from a modeled spectrum with log (N cm⁻²) = 17.5, σ_V = 0.61 km s⁻¹, and T_ex = 18 K, we built one thousand realizations of the observed spectrum with a Monte Carlo simulation, and we fitted the LTE model using the estimator proposed in Sect. 5.1. This Monte Carlo simulation allows us to numerically estimate the standard deviation on the estimated parameters and the correlation coefficient between log (N) and σ_V. This coefficient is−0.94, implying that the parameters are highly anti-correlated. However, the standard deviation on the log (N cm⁻²) and σ_V estimations are 0.017 and 0.005 km s⁻¹, respectively.This corresponds to typical relative errors of 4.0 and 0.8%, respectively. Hence some high (anti-)correlation does not necessarily imply that the model parameters can not be estimated, in contrast with a widespread intuition. While we illustrated this property with a given estimator, this statement is true for the CRB analysis. This emphasizes another of its interests. It provides standard deviations and coefficient of correlations without requiring to implement any Monte Carlo simulation.

The second and third row of Fig. 4 show important ambiguities (i.e., correlation coefficients close to one or minus one) between estimations of logN and T_ex and even more ambiguities between estimations of logN and σ_V (in particular for large values of N). When a single line is observed (see Figs. 4 d, e, g, h), the results for |γ(T_ex, logN)| and |γ(logN, σ_V)| are mostly larger than 0.9 for small N (in Figs. 4d, e, yellow pixelscorrespond to γ > 0.99 and in Figs. 4g, h to γ > 0.9). For high N, the ambiguity with T_ex decreases, but not the one with σ_V (in Figs. 4d, e, light blue values are − 0.5 < γ < 0 while in Figs. 4g, h dark blue values correspond to γ < −0.9 and γ < −0.99). The horizontal asymptote on the left of the CRB maps corresponds to a very sharp change of sign of correlation coefficients γ. Figures 4 f, i, show that, although some ambiguities remain between log N and σ_V for high N and small T_ex, having two lines mitigates these ambiguities in most cases.

As a rule of thumb, with a single line (see Figs. 4 a,b), for N > [10¹⁶−10¹⁷] cm⁻² (depending on T_ex), and T_ex ≥ [6−12] K (depending on N). With two lines (Fig. 4c), the situation greatly improves: for N > [10¹⁵−10¹⁷] cm⁻² (depending on T_ex) and T_ex ≥ 6 K. Figures 4a–c also show local minima of the when T_ex increases and N ≥ 10¹⁷ cm⁻², and when N increases and T_ex ≥ 6 K. To interpret these, the last row of Fig. 4 shows functions of the opacities. Comparing these with the variations of shows that the smallest values of (i.e., the best achievable precision) are mainly located at the area where τ₁ and τ₂ are close to 1. With two lines (see Fig. 4l), the best precision is when τ₁ < 1 and τ₂ > 1.

Fig. 5 Illustration of the correlation between N and σV estimations when a single line (13CO(1−0)) is available. The blue points in the scatter plots show the estimations of N and σV obtained with a Monte Carlo simulation of individual spectra that share the same physical parameters and different realizations of a white Gaussian noise with standard deviation σb = 100 mK. The parameters are Tex = 18 K, N = 1017.5 cm−2, ΔV = 1.1 km s−1, and σV = 0.61 km s−1.

Fig. 6 Noise standard deviation σb,ρ in mK whichensures that relative precisions are better than ρ% (for details see Eq. (25)). Top: analysis of a single line. Bottom: analysis of two lines 13CO(1− 0) and 13CO(2− 1). The contours for 10 and 100 mK are highlighted because these σb values bracket the values reached during typical observations at the IRAM-30 m. For this analysis, ΔV is fixed at 1.1km s−1, but the computations are done for four different values of σV (0.3, 0.6, 1.3,and 2.0 km s−1) and then projected on the (N, Tex) plane (see text for details).

4.7 Maximum noise tolerable to get a given relative precision on the different parameters

In the previous section, the standard deviation of the noise σ_b was fixed to 100 mK. Conversely, we now derive the amount of noise that guarantees a given relative CRB precision for T_ex, log N, σ_V and Δ_V. We compute σ_b,ρ the maximal value of σ_b that satisfies the following inequalities: (25)

where ρ is a fixed threshold. In other words, instead of analyzing the precision for a given amount of noise, it is also possible to analyze the tolerable level of noise to ensure an intended precision (herein described by ρ). Such an analysis is useful for designing an observation program and optimize the telescope time needed to reach the scientific goal.

Up to this point in the paper, we have checked the variations of the quantities as a function of T_ex and N with fixed values of Δ_V and σ_V. Figure 6 shows the variations of σ_b,ρ as a function of T_ex and N. As the computation of σ_b,ρ includes the computation of maximum values, it is possible to make the computations in three dimensions (with varying values of T_ex, N, and σ_V), and to project these on the (T_ex, N) plane. That is what is shown in Fig. 6 for different values of the relative precision ρ (5, 10, and 20%).

The IRAM- time estimator for the EMIR receivers⁵ indicates that we can achieve a sensitivity of 100 mK in 30–120 s at 110 and 220 GHz for a spectral resolution of 0.5 km s⁻¹. Similarly, we can achieve a sensitivity of 10 mK in one to three hours depending on the frequency and the observing mode (frequency or position switching). The colored contours thus correspond to the “fast/slow” acquisition mode at the IRAM-30 m. The comparison between the top and bottom lines allows us to see the gain in precision when analyzing the two lowest J lines of ¹³CO instead of a single one. In particular, the surface of reachable combinations of column density and excitation temperature more than doubles when analyzing two lines.

Instead of analyzing the level of noise, we can also analyze the peak-signal-to-noise ratio defined for one transition, l, by: (26)

and for two transitions by . Figure 7 shows the variations of minimum peak-signal-to-noise ratio for similar conditions as in Fig. 6. Analyzing only the ¹³CO(1−0) line requires at least a S/N of 100 to get a relative precision of 20%. Adding the ¹³CO(2− 1) line in the analysis reduces the minimum S/N by a factor of up to ten to reach the same relative precision. The required S/N increases at high columndensities because the lines become optically thick, and at a combination of low column density and high excitation.

Fig. 7 Same as Fig. 6, except that the contours show the variations of the peak-signal-to-noise ratio, , required to reach a given relative accuracy (ρ%).

4.8 Complementing ¹²CO(1–0) observations

All the previous analyses were done for the ¹³CO isotopologue because it enabled us to study the case of low J lines that experience the transition from optically thin to thick regime over the range of column densities and excitation temperatures that are found in molecular clouds. However, the targeted transition when observing the molecular gas of a new astronomical source is usually ¹²CO(1− 0) because it is the strongest line in the easily observable 3mm atmospheric window (Wilson et al. 1970; Pety et al. 2017).

Here we asked ourselves two questions. The first question deals with what J line is best for carrying out observations in order to reach a relative precision of 20% on all the estimated parameters. Figure 8 shows the variations of the maximum noise σ_b,ρ when a single line is observed among the first six rotational transitions of ¹²CO. A global pattern is seen, especially for the higher energy transitions ¹²CO(4−3), ¹²CO(5− 4), ¹²CO(6− 5), For instance, if N lies in the interval [10¹⁷, 10¹⁹] cm⁻² and T_ex > 24 K, the ¹²CO (6−5) line seems the best choice (from the point of view of the CRB) because it allows us to reach 20% accuracy over this broad range of parameters for a noise level of 100 mK. However, this line is not easy to access from ground-based telescopes because of the limited atmospheric transmission at the line frequency of 690 GHz.

The second question considers the best J line for carrying out observations to complement the J = 1−0 line to reach the same relative precision of 20%. Figure 9 seems to indicate that observing ¹²CO(6−5) would be the most useful as it would allow to tolerate a noise level of σ_b = 300 mK and maintain a good level of precision for N in the interval [10^16.0, 10^18.5] cm⁻². We stress that this result only applies to the case where all transitions have the same excitation temperature. In practice, deviations from a Boltzmann population may be present leading to different excitation temperatures for the ¹²CO transitions (van der Tak et al. 2007) because of the higher critical densities of the higher-J transitions. Nevertheless, the usefulness of mildly excited lines remains valid. Non-LTE approaches will be developed in the future that will provide a quantitative assessment of the diagnostic power of these lines.

5 Application to the ORION-B data

The CRB is only a lower bound on the variance of any unbiased estimator. Once the orders of magnitude of the CRBs have been analyzed, the next step is to find a good estimator of physical parameters. In this section, we first propose such an estimator and analyze its performance for a realistic amount of noise (herein chosen to σ_b = 100 mK), before applying it to the ORION-B data.

5.1 Proposed estimator

The maximum likelihood estimator (MLE) is a good candidate because under mild conditions, it reaches the CRB asymptotically, that is, when σ_b ↦ 0 (Garthwaite et al. 1995). Appendix B details the computation of this estimator, its initialization, and the iterative algorithm used to yield the estimation that is noted . We also briefly discuss its computational efficiency.

Appendix C analyzes the performance of the proposed estimator on simulated data. This appendix shows that this estimator performs optimally for pairs of (T_ex, N) values such that a relative precision of reference is reached for all estimated parameters (that is, the conditions of Eqs. (25) are satisfied with ρ = 20%). Interestingly, this appendix also shows that the obtained estimation can be injected in the CRB computation to detect whether or not this estimation is accurate.

Fig. 8 Noise standard deviation σb,ρ in mK which ensures relative precisions better than 20% when a single line of 12CO is analyzed. Other details are identical to Fig. 6.

Fig. 9 Noise standard deviation σb,ρ in mK which ensures relative precisions better than 20% when a couple of 12CO lines are observed: J = 1−0 and a higher J line. Other details are identical to Fig. 6.

5.2 Estimation of the number of velocity components and initialization of the parameters

The number of velocity components is a priori unknown. While we could have tried to use our maximum likelihood estimator to fit the data with either one or two components, we would then have had to devise a statistical test to determine which assumption to choose. Instead, it is simpler to check for the presence of several local maxima in the spectrum denoised with ROHSA, the technique mentioned in Sect. 2. In particular, this allows us to have a spatially coherent detection of the number of components. The ROHSA algorithm is applied separately on the ¹³CO(1− 0), C¹⁸O(1−0), and ¹²CO(1− 0) lines. If one of the denoised spectra for ¹³CO(1−0), C¹⁸O(1−0) or ¹²CO(1− 0) has at least two local maxima in the velocity interval of interest [8.25, 14.25] km s⁻¹, we then fix the number of velocity components to two for all three species. The denoised spectra ¹³CO(1− 0), C¹⁸O(1−0) and ¹²CO(1− 0) are analyzed iteratively in this order and as soon as two components are selected based on one of the three spectra, the velocities associated to the local maxima are used to initialize the estimations of the velocity Δ_V of each component for all three species. This ensures that the velocities of each component stay compatible for all three species during the fit. We analyze the denoised spectra in the above order because the ¹³CO(1− 0) line has both a good S/N and moderate opacities. The C¹⁸O(1−0) line delivers a good information on the underlying velocity structure because it is most often optically thin, but its limited S/N may hamper the Δ_V initialization. The saturation that happens for the ¹²CO(1−0) line also makes the determination of the Δ_V initializations inaccurate. Finally, when all the three denoised spectra have only one maximum, one component is independently fitted per species.

When initializing the parameters before maximizing the likelihood, we use two different assumptions to help the algorithm to converge towards reasonable solutions. First, when two components are detected, we use the same Δ_V and σ_V initializations for all the species so that the estimated parameters for each component remain correctly paired among the three species, as explainedabove. The white contours in Fig. 13 delimit the regions where two local maxima have been detected. Only 23% of the field of view requires two velocity components. As the estimations of T_ex, N, and σ_V are highly correlated (see Sect. 4.6), we systematically search in the 3D grid described in Sect. sectionlinking B.2 to initialize them. We stress that this is only during the initialization process that we use the same values Δ_V and σ_V for ¹³CO and C¹⁸O. The maximization of the log-likelihood is done independently on each species, ensuring that the estimations of Δ_V and σ_V may take a different value for each species.

Second, a single line of ¹²CO is available and it is quite optically thick. In our analysis of the CRB, we observed in Figs. 4 (d, e and g, h) that the estimation of logN is highly correlated with σ_V, when N is high. This may imply some degeneracy between the estimation of the velocity dispersion and the column density. To alleviate this issue, we first deal with ¹³CO(1−0), ¹³CO(2− 1), C¹⁸O(1−0), C¹⁸O(2−1), and we use the estimation of σ_V obtained on ¹³CO to fix σ_V for the estimation of the other parameters (T_ex, log N and Δ_V) in the analysis of the ¹²CO(1−0) line. If this assumption is false, the obtained estimations of the other parameters will be biased. While this procedure is not ideal, we empirically obtained estimations of T_ex and logN much closer to physical intuition: in particular, the estimations of the column density are 100 times too large when all four parameters are estimated. An analysis of the impact of a potential incorrect value of σ_V goes beyond the scope of the present paper.

5.3 Detailed analysis of two lines of sight

Figure 10 shows how the proposed estimator succeeds to fit the C¹⁸O, ¹³CO, and ¹²CO low J lines towards two lines of sight in the studied field of view (see red crosses Fig. 1), at offsets (803″, 473″) and (578″, −121″). The spectra on the left column are modeled with a single velocity component but they are asymmetric. This implies that our model is not perfectly adequate because it assumes that the line profile is symmetric. The issue is most problematic for ¹²CO, because the high opacity and the complex underlying velocity field imply a more complex profile with broad wings on each side of the line. A detailed solution for this issue is beyond the scope of the present paper.

The spectra on the right column are well-fitted with two different velocity components. The estimations for the C¹⁸O and ¹³CO species are physically relevant because the two velocity components are well separated in velocity. This is less obvious for ¹²CO, which presents a large velocity overlap of the two components.

5.4 Estimation of the quality of the fit and filtering out inaccurate estimations

In this section, we first compute the “energy” and the standard deviation of the fit residuals as two ways to quantify the quality of the fit. We then explain how we filter out inaccurate estimations from the physical analysis.

After a fit, we can define three different “energies” in the sense of the information theory.

• The “energy” of the measured signal is

(27)

• The “energy” of the estimated signal is

(28)

depending of the number of estimated components. The “energy” of the fit residual is (29) (30)

All the sums are computed on an interval of 10 km s⁻¹ around the maximum. The fit quality can then be quantified by comparing the “energy” in the residual with either the “energy” in the observed spectrum or the difference of “energy” between the observed and estimated signals . The former formula tells us the fraction of the observed “energy” that has not been fitted. The latter formula tells us whether the observed spectra has been under-fitted (positive value) or over-fitted (negative value). Measuring the residual “energies” is similar to computing a χ² in least-square fitting. Another way to quantify the quality of the fit is to compare the standard deviation of the residuals with the noise standard deviation (σ_b) on the observed spectrum. The fit is good when σ_r ~ σ_b.

We can, in addition, use the CRB framework to filter out pixels with inaccurate estimations. As explained in Appendix C.1, an estimation would be considered inaccurate when there is at least one estimation among the 3 × 3 neighboring pixels, for which at least one of the following conditions is satisfied: (31)

Fig. 10 Two examples of LTE fit of the CO isotopologues lines. In the titles, Tex is expressed in Kelvin, N in cm−2, ΔV in km s−1, and σV in km s−1. The plain lines show the data, and the dotted ones show the fit results. Values in red indicate estimations whose relative precision is larger than 20%. The associated lines of sight can be localized in Fig. 1 (see red crosses).

Fig. 11 Jointhistograms of the standard deviations of the residuals (σr) and of the noise (σb) for all the studied lines. Standard deviations are expressed in Kelvin. The dashed lines correspond to ratios 1 and 10.

5.5 Global analysis of the quality of the estimation

Figure 11 compares the standard deviation of the residuals (σ_r) with the standard deviation of the noise (σ_b) for all the lines studied here. If the fits were ideal, the joint histograms would only peak near a line of slope one. They thus suggest that the C¹⁸O lines are better fitted than the ¹³CO lines, and that the ¹²CO lines are the least well-fitted. We also checked that the residuals are larger (i.e., σ_r > σ_b) when the energy is large (not shown in the figure). As the S/N also increases with , this issue implies some mis-specification in the model. The best fit happens for the C¹⁸O(1–0) line that has the lowest opacity. In that case, the profiles Ψ are almost perfect Gaussian profiles. On the contrary the high opacity of the ¹²CO(1–0) line implies that the profiles are highly saturated, and any small kinematic perturbations create deviations from a Gaussian profile as shown in Fig. 10, where the asymmetry pattern could not be taken into account by the model. A better modeling could thus require additional velocity components especially for the ¹²CO(1− 0) line. Another limit of the model is that it does not encode self-absorption signatures that may happen at large opacity.

Figures 12 a, b show the spatial distributions of the “energy” of the observed spectra and of the fit residuals for the ¹³CO(1−0) line. Both images share the same color look-up table. On the left image, the yellow pixels correspond to bright molecular gas while the blue pixels corresponds to faint signal or noise associated with the IC 434 HII region. The “energy” of the residual still exhibits spatially coherent structures, but at a much smaller level than the “energy” of the measured spectra. Our estimator under-fits the observed spectra as shown in Fig. 12d. This may be related to the fact that our model does not fit asymmetric profiles. However, the fit quality is good as illustrated by the image of the energy ratio, which shows that the residual “energy” amounts to less than 1% of the signal (the dark blue color corresponds to 0.01% in image Fig. 12c) except in regions where the S/N becomes small. Figure D.1 shows the same quantities for all the lines modeled in this paper. Overall, the fitting method is able to recover all the emission with differences at the percent level or less for all lines.

Increasing the complexity of the line profile model to address the observed misspecifications could be hazardous. While the increase of the number of parameters in more complex models certainly allows one to decrease the difference between the observed spectrum and the model, it may also increase the variance of estimations in such proportions that obtained estimations may become useless. In other words, the simple model used here does not capture all the complexity of the physical processes, but it at least allows us to capture the processes that has already encoded. The fact that the residuals are smaller than 1% compared to the observed signal is sufficient to make the analysis of the excitation temperature, the column density, and the velocity dispersion pertinent for CO isotopologues. The main source of systematic errors in the column density determination results from the deviations from the local thermodynamic equilibrium, leading to a more complex partition function than the simple formula in Eq. (9). The effect is expected to be stronger in warm regions (T_ex ≥ 50−100 K) where many rotational levels are populated and contribute to the partition function. These warm regions occupy a small fraction of the total volume and therefore a bias would not affect the general conclusions. Non-LTE approaches will be developed to assess more quantitatively the magnitude of the effect and to provide recommendations on the best method depending on the molecular lines and the range of physical conditions that are studied.

Fig. 12 Spatial variations of the observed spectrum “energy” (a), of the residual “energy” (b), of their ratio in % (c), and of the ratio of “energy” in % that has not been modeled (d). The unit of the color look-up table is Kelvin2. White contours show the regions where two components have been detected.

6 Astrophysical implications

The proposed estimator (see Sect. 5.1) provides accurate column densities, excitation temperatures, and velocity dispersions in the framework of LTE excitation and radiative transfer. This allows us to carefully analyze the errors introduced by the simpler hypotheses that are commonly used for deriving CO isotopologues column densities.

6.1 Estimated parameters and associated uncertainties for C¹⁸O, ¹³CO, and ¹²CO

Figure 13 shows the spatial variations of the estimated parameters and associated uncertainties for the C¹⁸O, ¹³CO, and ¹²CO isotopologues. Table 2 lists the mean and standard deviation values of the excitation temperatures, column densities, velocity dispersions and opacities. These values are computed over the field of view that is observed for all the lines.

The peak-signal-to-noise ratio of the ¹²CO(1−0), ¹³CO(1− 0), and ¹³CO(2− 1) lines is large (> 20) over most of the studied field of view. The regions where the peak-signal-to-noise ratio of the C¹⁸O(1−0) and C¹⁸O(2−1) is larger than 20 still amounts to respectively 15 and 32% of the studied field of view. The CRBs are small for all estimated parameters (inside the red contours that delimit the regions where the relative precision is better than 20% for all estimated parameters) except near the regions of transitions between one and two velocity components (i.e., near the white contours). Even though the ¹²CO S/N is much larger than the C¹⁸O or the ¹³CO one, the ¹²CO estimations are more uncertain (see map of for ¹²CO in Fig. 13) because a single transition is available and the ¹²CO opacities are large (ranging from 5 to about 1000). However, the higher S/N for ¹²CO helps to derive the velocity field in regions where ¹³CO and C¹⁸O are not well detected (diffuse gas).

The largest variations are observed in the column density which varies from the detection limit near 10¹⁵ cm⁻² up to values larger than 10¹⁷ cm⁻² for ¹³CO. The comparison of the column density maps for ¹³CO and C¹⁸O suggests that the C¹⁸O molecules are confined to the high column density regions and avoid the cloud edges. The ¹²CO isotopologue shows a different behavior with emission extending over most of the imaged field of view and column densities ranging from ~10¹⁶ to ~ 10¹⁹ cm⁻². These behaviors are related to the difference of opacity of the isotopologues lines. Both ¹³CO lines have moderate opacities across the mapped region, with somewhat higher values for ¹³CO(2− 1) than for ¹³CO(1− 0). The ¹²CO is highly optically thick almost everywhere.

For all CO isotopologues, the excitation temperature presents coherent spatial variations with maximum values near the NGC 2023 star forming region. Even though the C¹⁸O emission is fainter and less extended than that of ¹³CO, both speciesprovide similar kinematics information. The velocity field is spatially regular with well-resolved gradients, for instance in the Horsehead nebula at the bottom right of the map. The velocity dispersion σ_V ranges from 0.4 to 0.8 km s⁻¹, with values around 0.3 km s⁻¹ in the Horsehead nebula (see also Hily-Blant et al. 2005), and somewhat narrower lines for C¹⁸O than for ¹³CO.

Fig. 13 Top: spatial variations of the estimated physical parameters. From left to right: centroid velocity, velocity dispersion, excitation temperature, column density, and line opacities. Bottom: spatial variations of the relative precisions. From left to right: relative precision on the centroid velocity, velocity dispersion, excitation temperature, column density, and line peak S/N. The black contours on the peak-signal-to-noise-ratio image delimit the regions where . On all images, red contours delimit the regions where the relative precision is better than 20% for all estimated parameters, and white contours delimit regions where two components have been estimated. In this lattercase, the images only show the estimation that is the closest (in terms of centroid velocity ΔV) to its neighboring pixels.

Fig. 14 Comparison between the map of the dust temperature and maps of temperature ratios. Inaccurate estimations are filtered out. The dust temperature is only presented in regions with an accurate estimation of parameters.

6.2 Excitation temperatures

Simplifying assumptions are often made when analyzing the CO rotational emission. The most usual one is that all isotopologues have the same excitation temperature. It is based on the similarity of the collisional cross sections. Because the opacity of a rotational transition scales with the molecular column density, the differences in abundances translate to different opacities. The main isotopologue (¹²CO) has optically thick lines while the rarer isotopologues (¹³CO and C¹⁸O) have lines either optically thin or with moderate opacities (see, e.g., Ripple et al. 2013). Another assumption for high density regions (n > 10⁵ cm⁻³ Goldsmith & Langer 1978) is the full thermalization of the lowest CO rotational levels at the temperature measured on dust (i.e., assuming the convergence of the dust and gas temperatures). Both hypotheses have been subject to scrutiny. Recently, in their clustering analysis of a one square degree map in the Orion B cloud, Bron et al. (2018) showed that the observed CO isotopologue line intensities and line ratios cannot be explained using these simple hypotheses and that differences in excitation temperatures should be taken into account.

As shown in Fig. 14, the excitation temperatures of the three isotopologues are different across the field of view. The ¹²CO isotopologue has the largest excitation temperature, followed by ¹³CO and C¹⁸O. Table 3 lists the typical ratios of excitation temperatures. They are COCO) ~ 1.7 and CO)∕T_ex(C¹⁸O) ~ 1.3. Figure 14 shows that these ratios vary significantly as a function of the total column density and dust temperature. Figure 15 suggests that higher dust temperature regions that trace higher UV illumination conditions tend to show larger differences in excitation temperatures between the CO isotopologues. The same trend is seen when the CO isotopologue excitation temperatures are compared to the dust temperature as in Fig. 14. A similar effect has been reported by Welty et al. (2018) for the diffuse-translucent cloud along the line of sight to HD 62542.

The excitation temperature of C¹⁸O, which traces the UV shielded regions, is on average lower than T_dust (the mean value of T_ex∕T_dust is 0.71, see Table 2). We find a similar situation for ¹³CO(1–0) (mean value T_ex ∕T_dust = 0.76), while most of the positions show ¹²CO(1–0) excitation temperatures larger than T_dust. This difference between the CO excitation temperature – which is a lower approximation of the gas kinetic temperature – and the dust temperature is indeed expected in photodissociation regions and UV-dominated regions where the gas kinetic temperature is larger than the dust temperature. This is different from the usual approximation that T_dust is a good approximation of the gas kinetic temperature in cold and shielded regions.

The difference in excitation temperatures between CO isotopologues can be explained by radiative trapping in the ¹²CO lines or by the presence of kinetic temperature gradients along the line of sight, especially near photodissociation regions, possibly combined with density gradients. Clear spatial patterns emerge in Fig. 14. The C¹⁸O and ¹³CO excitation temperatures get closer to T_dust in regions where T_dust is lower than 20 K. Because the dust emission strongly varies with its temperature (as T^(4+β) where β is the dust emissivity index and takes values between 1.5 and 2, Planck Collaboration XI 2014), the dust temperature derived from a single temperature fit of the spectral energy distribution in lines of sight combining a strongly UV illuminated region and a more shielded material does not represent the conditions in the UV shielded region well. The dust temperature can overestimate the temperature in the shielded gas that represents the bulk of the material. Somehow the illuminated region “overshines” as compared to the bulk of the matter.

The error introduced in the column density determination by using an incorrect excitation temperature can be estimated by examining the variation of the line column-to-intensity ratio, defined as the ratio of the column density of the species to the line integrated emission, as a function of the excitation temperature. Figure 16 shows the variation of the column-to-intensity ratio for the ¹³CO(1−0) and ¹³CO(2− 1) lines for different column densities of ¹³CO, assuming a velocity dispersion of 0.61 km s⁻¹. Column densities of correspond to optically thin lines, while the opacity becomes significant (i.e., τ ≥ 0.5) for 10¹⁶ −10¹⁷ cm⁻². In the optically thin case, the column-to-intensity ratio presents a shallow minimum which depends on the transition, rises rapidly at temperatures lower than the minimum and more slowly for excitation temperatures above the minimum. When the line opacity becomes significant, the shape of the column-to-intensity ratio curve changes andthe minimum shifts to higher excitation temperatures or possibly disappears. Therefore, using the ¹²CO excitation temperature for determining the column densities of ¹³CO and C¹⁸O leads to errors in the estimation of these column densities because the associated column-to-intensity ratio is inappropriate. For moderate column densities ( cm⁻²), the column densities can be underestimated in the case of low excitation temperatures, or overestimated when using a too high excitation temperature, depending on the position of T_ex relative to the minimum of the column-to-intensity curve. A CO excitation temperature near the minimum of the column-to-intensity curve minimizes the error because a small difference in T_ex in this region does not change the column-to-intensity. For high column densities (N ~ 10¹⁷ cm⁻²), the bias is significant at low excitation temperatures (T_ex < 20 K) because the column-to-intensity ratio (purple curve in Fig. 16) is rising fast at low excitation temperatures.

Fig. 15 Scatter plots between the CO isotopologue excitation temperatures. The color scale encodes the dust temperature Tdust. The ellipsesrepresent the interval of confidence for each estimation: Each ellipse is centered on the estimation of the excitation temperature for one pixel and its horizontal and vertical sizes are equal to the associated CRBs. Dashed ellipses correspond to pixels with two components. Inaccurate estimations are filtered out. Dashed red lines show the loci of ratios 1/2, 1, 2, and 4.

6.3 Abundances

Figure 17 presents maps of the CO isotopologue abundances relative to H₂ (see Sect. 2.2), and Fig. 18 shows scatter plots of the relationships between the CO isotopologue column densities and that of molecular hydrogen. In the mapped area, ¹³CO and ¹²CO can be fitted and analyzed for H₂ column densities larger than 10^21.5 cm⁻². The threshold for C¹⁸O is about twice higher at ~10^21.8 cm⁻². Indeed, as shown by Pety et al. (2017) and Orkisz et al. (2019), the threshold for the detection of C¹⁸O is A_V ~ 3 mag or N(H₂) = 10^21.5 cm⁻² in the Orion B molecular cloud, while the thresholds for ¹²CO and C¹³O are close to A_V = 1 mag.

Over the mapped area, the mean abundances are well defined at CO)∕N(H₂) = 10^{−5.6 ± 0.29} = 2.5 ± 1.5 × 10⁻⁶, and N(C¹⁸O)∕N(H₂) = 10^−6.7±0.15 = 2.0 ± 0.8 × 10⁻⁷ (see Table 2). These mean abundances are comparable to those of other molecular clouds in the solar neighborhood such as Taurus and Ophiuchus (Frerking et al. 1989). However, these values are about a factor of two lower than those predicted when no isotopic fractionation is assumed and when the non depleted gas phase carbon elemental abundances applicable to the Orion region are used, namely, C∕H = 1.4 × 10⁻⁴, ¹²C/¹³C = 57−67 and ¹⁶O/¹⁸O = 500−560 (Gerin et al. 2015; Langer & Penzias 1990; Wilson & Rood 1994), namely CO)∕N(H₂) = 4−5 × 10⁻⁶, and N(C¹⁸O)∕N(H₂) = 5−6 × 10⁻⁷. The difference is more pronounced for C¹⁸O because its abundance is affected by both photodissociation and freeze-out over a more significant fraction of the studied area.

Significant deviations from the mean values are present. The C¹⁸O abundance is not only lower near the photo-illuminated edges where molecules are photodissociated, but also in high column density and well shielded regions. In these latter regions, the dust temperature gets below the CO condensation temperature (T_dust < 17 K), and the CO molecules can rapidly freeze onto dust grains, lowering the gas phase abundances. This depletion effect is seen for C¹⁸O and ¹³CO supporting the explanation by a global freeze-out effect.

Because the ¹²CO lines are saturated over most of the region, and because a single line has been observed, the determination of its abundance is more uncertain. Nevertheless, fixing the value of the velocity dispersion of ¹²CO to the one estimated for ¹³CO (see Sect. 5.2) allows us to determine the ¹²CO column density for the pixels that have the least saturated line emission. Although fairly uncertain, the resulting values of the ¹²CO abundance relative to H₂ are close to 6 × 10⁻⁵, lower than the expected value using the carbon elemental abundance relative to H, 1.4 × 10⁻⁴ (Gerin et al. 2015), but similar to values obtained in Taurus by Pineda et al. (2010).

Fig. 16 Plots of the 13CO column densities per unit intensity of 13CO(1−0) (top) and 13CO (2−1) (bottom) as a function of the excitation temperature. This plot is done for four values of the 13CO column density, and σV = 0.61 km s−1.

Fig. 17 Comparison of the map of dust-traced H2 column density with maps of various column density ratios. Inaccurate estimations were filtered out. The dust-traced column density is only presented in regions with an accurate estimation of parameters.

Fig. 18 Scatter plot of the CO isotopologue column densities as a function of the dust-traced H2 column density. Inaccurate estimations are filtered out. The black dashed lines show the expected gas phase abundances relative to H2 with no depletion: C18O/H2 = 5 × 10−7, 13CO/H2 = 4 × 10−6, and CO/H2 = 1.4 × 10−4.

6.4 CO isotopologue column density ratios

Maps of the CO isotopologue column density ratios are shown in Fig. 17 and scatter plots are displayed in Fig. 19. With no isotopic fractionation and all carbon locked in CO, the expected CO isotopologue ratios are ¹²CO/¹³CO = 57−67 and ¹³CO/C¹⁸O = 7.5−9.8 using the elemental abundances given in the previous subsection. As shown in Fig. 19, the lower bound of the ratio of the ¹³CO and C¹⁸O column densities is indeed in the expected range at CO)∕N(C¹⁸O) = 8. The upper bound is close to 50 indicating that chemical effects play a significant role, by enhancing the ¹³CO abundance (fractionation) or destroying C¹⁸O (photodissociation).

Figure 20 suggests that the column density ratio is well correlated with the ratio of line intensities, but the column density ratio is a factor up to 1.75 smaller than the ratio of line intensities because of the difference in excitation temperatures and of the moderate opacity of the ¹³CO(1–0) line. Ratios of integrated intensities can therefore be used to estimate the column density ratios, but after checking with radiative transfer calculations for a possible multiplicative bias and introducing a correction if needed.

Although the derivations of ¹²CO column densities are uncertain, the column density ratio COCO) ranges between 10 and 60. In particular, the low values of the ratio remain even when the opacity of the ¹²CO line becomes small enough to accurately derive the column density. Such low values are found in diffuse and translucent gas as a consequence of efficient fractionation in ¹³C due to the exchange reaction between ¹³C⁺ and ¹²CO, that enhances the ¹³CO abundance (Liszt & Pety 2012). The physical conditions in the translucent envelope of Orion B seem to favor fractionation, which is not restricted to a small subset of the mapped area but it seen over wide areas. As discussed by Bron et al. (2018) the presence of chemical fractionation over the whole region can be identified by comparing the ratio of integrated intensities of the CO isotopologues, and looking at the data in the W(¹³CO)/W(¹²CO) versus W(¹³CO)/W(C¹⁸O) plane. The existence of this chemical fractionation implies that using a single value for the abundance ratios of CO isotopologues can introduce significant errors when attempting to correct for the CO opacity in computing its column densities as done by Barnes et al. (2018) for instance. This further increases the bias introduced by using the same excitation temperature for ¹²CO and ¹³CO ground state transitions. In addition to possibly biasing the results, using such simplifying assumptionsis also expected to increase the dispersion and affect the overall determination of the X_CO = N(H₂)∕W(CO) conversion factor.

Fig. 19 Scatter plots between CO isotopologue column densities. The color scale encodes the ratio of the excitation temperatures of the considered species in each panel. The ellipses represent the interval of confidence for each estimation: Each ellipseis centered on the estimation of the column density for one pixel and its horizontal and vertical sizes are equal to the associated CRBs. Dashed ellipses correspond to pixels with two components. Inaccurate estimations are filtered out. The dashed red lines show the loci of ratios 5, 20, and 60.

Fig. 20 Joint histogram of the logarithm of the estimated column density ratios and of the observed integrated intensity ratios. All the ratios are computed for the 13CO(1− 0) line over the C18O(1−0) one. The dashed red lines correspond respectively to ratios of one and two.

6.5 Velocity dispersions

Maps of the velocity dispersions for ¹³CO, C¹⁸O, and ¹²CO are shown in Fig. 13. The mean values are listed in Table 2 and the ratios for the different CO isotopologues are listed in Table 3. Our formalism includes the line broadening due to opacity (see, e.g., Phillips et al. 1979), which is very significant for ¹²CO. This implies that the actual velocity dispersion derived from high opacity lines like those of ¹²CO is smaller than the apparent line width. Because we fixed the ¹²CO velocity dispersion to the value obtained with ¹³CO in our estimation, the velocity dispersions of ¹²CO and ¹³CO are identical (see Sect. 5.2). The velocity dispersions of C¹⁸O and ¹³CO are howeverfitted independently. Both species show similar velocity dispersions but ¹³CO has consistently broader line profiles than C¹⁸O. The ratio between the velocity dispersions of ¹³CO and C¹⁸O is 1.1 (see Table 3). The ¹³CO emission is produced by a more extended volume along the line of sight than the C¹⁸O emission as shown by the lower threshold in N(H₂) where ¹³CO is detected as discussed above. When analyzing the filamentary structure of the Orion B molecular cloud, Orkisz et al. (2019) showed that the gas velocity dispersion determined from C¹⁸O reaches a minimum value in the filament ridges, and is always lower than the velocity dispersion determined by ¹³CO. The refined analysis presented here, which takes the opacity broadening effect into account, confirms the presence of gradients in velocity dispersion across the spatial extent of the cloud and along the line of sight, which are captured by the difference between ¹³CO and C¹⁸O. Inspecting the spatial distribution of the velocity in Fig. 13 suggests that the small excess of velocity dispersion for ¹³CO relative to C¹⁸O is more prominent in the regions with relatively low T_dust. This supports the hypothesis that this variation of velocity dispersion is tracing the starting point of the dissipation of turbulence when entering the dense filamentary skeleton of the molecular cloud.

7 Conclusion

This paper presents an analysis of the precision of the estimation of physical parameters (Δ_V, σ_V, N, T_ex) when trying to fit spectra of low J transitions for the most common CO isotopologues using the LTE radiative transfer model. This analysis was based on the Cramer Raobound (CRB) computation. We applied this analysis to the region of the Orion B molecular cloud that contains the Horsehead pillar, as well as the NGC 2023 and IC 434 HII regions with the following astrophysical results.

• The mean abundances of the CO isotopologues are consistent with previous determinations in other regions: X(¹²CO) ~ 6 × 10⁻⁵, X(¹³CO) = 2.5 ± 1.5 × 10⁻⁶, and X(C¹⁸O) = 2.0 ± 0.8 × 10⁻⁷.

• The excitation temperatures T_ex are different among the CO isotopologues. ¹²CO presents the highest T_ex, followed by ¹³CO and C¹⁸O. For ¹³CO and C¹⁸O, the excitation temperatures are lower than the dust temperature on average, while they are higher for ¹²CO. This systematic variation can be understood as resulting from gradients of physical conditions along the line of sight together with the increased effect of radiative trapping for the more abundant isotopologues.

• These differences in T_ex imply that the ratio of ¹³CO(1−0) and C¹⁸O(1−0) integrated intensities is not a direct measurement of the column density ratio , with a difference of up to a factor two. Moreover, this column density ratio exhibits regular spatial variations across the mapped region, with high values in the UV illuminated regions and low values in shielded regions. These low values are consistent with the ratio of ¹³C and ¹⁸O elemental abundances (i.e., a factor of about 8).

• In this nearby molecular cloud, the elemental abundances are uniform and the variations in CO isotopologue relative abundances are solely due to chemical processes (fractionation, photodissociation, freezing). The interpretation of variations of line integrated intensity ratios should therefore be performed with caution, taking into account radiative transfer and chemical effects.

We obtained the following results from the methodological viewpoint.

• This analysis has shown that it is important to take the opacity broadening effect into account when fitting the line profiles, even for moderate opacities as first discussed by Phillips et al. (1979). The estimation of the column density is correlated with that of the velocity dispersion when the line is optically thick and this correlation between column density and velocity dispersion must be taken into account when estimating uncertainties on the fitted parameters even for moderate line opacities. When the line is optically thin, the estimation of the column density is correlated with that of the excitation temperature except in a small interval where the ratio of the column density of the species to the integrated intensity of the line reaches a minimum (around T_ex = 10 K for the 1–0 transitions of CO isotopologues). This means that a small variation of the estimation of the excitation temperature or velocity dispersion leads to large errors on the estimation of the column density.

• This analysis also allows us to quantify the benefit of a simultaneous analysis of two rotational lines of the same species compared to the analysis of a single line. In particular, it alleviates the degeneracy described above. This is a rigorous demonstration of intuitive results. It is an argument in favor of the installation of dual-band receiver systems for telescopes like the IRAM-30 m or NOEMA.

• In order to derive the precision achieved on these parameters when trying to fit actual observations of the CO isotopologue lines towards the Orion B molecular cloud, we first showed that a (simple) maximum likelihood estimator is unbiased and efficient when the relative precision given by the CRB is better than 20%, and that it is possible to detect pixels for which the estimation of the parameters in LTE conditions is accurate (i.e., better than 20%).

• The residuals of the fit of the CO isotopologue lines amount to less than 1% of the original signal and the relative precision on the physical parameters is better than 20% for 63, 82, and 40% of the field of view for the ¹²CO, ¹³CO, and C¹⁸O species, respectively. The presence of structured residuals nevertheless indicates that the model remains sometimes too simple. In particular, asymmetric line profiles or the presence of line wings are incorrectly modeled. The ¹²CO line profiles are the most affected. Addressing the possibility of catching the complex shape of this spectrum is a motivating perspective that would generalize the approach initiated in this paper.

In the transition between regions where the number of required velocity components changes, some ambiguity on the velocities of the components occurs, and this impacts the estimations of all the other parameters. Fixing a priori Δ_V based on spatial processing (e.g. extending the ROHSA pre-processing) and thus applying only a gradient on σ_V, T_ex and N could improve the robustness of our estimations. This would be useful, in particular, for low S/N pixels. Another perspective is to use a spatial regularization criterion in the fit to improve all the estimations. This will be the subject of another paper (Vono et al., in prep.). Finally, trying to estimate the above physical parameters in regions that are more diffuse than on the studied field of view or for other species that have higher dipole moments (HCO⁺ or CS), requires the use of non-LTE models. Such non-LTE models would also be interesting to identify possible systematic effects coming from the use of the LTE approximation. An upcoming paper will study the large velocity gradient approximation of the radiative transfer using a similar CRB approach.

Acknowledgements

This work is based on observations carried out under project numbers 019-13, 022-14, 145-14, 122-15, 018-16, and finally the large program number 124-16 with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research also used data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). We thank CIAS for their hospitality during the many workshops devoted to the ORION-B project. This projecthas received financial support from the CNRS through the MITI interdisciplinary programs. This work was supported in part by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP, co-funded by CEA and CNES. JRG thanks Spanish MICI for funding support under grant AYA2017-85111-P. Finally, we thank the anonymous referee for helpful comments on the manuscript.

Appendix A: Gradient calculation

The spectrum at frequency ν is: (A.1)

and thus it isa function of the unknown parameters θ = [T_ex, logN, Δ_V, σ_V]. The following gradients are useful to derive the Fisher matrix (see Eq. (14)).

A.1 ∂s (ν)/∂T_ex

(A.2)

A.2

(A.3)

Thus (A.4)

and, finally, (A.5)

A.3

(A.6)

Thus (A.7)

A.4

(A.8)

Thus, (A.9)

where is numerically computed.

A.5 ∂s (ν)∕∂L_N

Let us introduce L_N = logN. (A.10)

A.6

Since N = 10^L_N and using Eq. (A.8), we get (A.11)

where ln is the natural logarithm.

A.7 ∂Ψ(ν)∕∂L_N

From Eq. (A.6), we get (A.12)

Then using Eq. (A.11), one gets (A.13)

A.8 ∂s (ν)∕∂Δ_V

(A.14)

A.9 ∂Ψ(ν)∕∂Δ_V

With and , we have (A.15)

Thus, (A.16)

Since , we get (A.17)

where (A.18)

Thus (A.19)

A.10 ∂s (ν) /∂σ_V

(A.20)

A.11 ∂Ψ(ν)∕∂σ_V

From Eq. (A.15), we get (A.21)

and with , we get (A.22)

Then, (A.23)

Thus (A.24)

which can also be written (A.25)

Appendix B: Maximum likelihood estimator

B.1 Definition for two lines of the same species and a single velocity component

We begin with the assumption that we are estimating the physical parameters (θ) of the LTE radiative transfer for two lines of the same species, and a single velocity component. Starting from Eq. (13), we write out the notation of the 2 K samples x_n,l as χ. The amount of noise for each line is fixed to σ_b,1, and σ_b,2. In this case,the maximum likelihood estimator (MLE) is (Garthwaite et al. 1995): (B.1)

Thus, is the argument that maximizes the likelihood for the observed sample χ. With two lines, the likelihood can be written as: (B.2)

And the log-likelihood is: (B.3)

B.2 Initialization of the unknown parameters

The log-likelihood function L(θ;χ) can have many local maxima. It is thus crucial to initialize the gradient near the global maximum. As explained in Sect. 4.1, the vector of unknown parameters has four components in the case of a single velocity component. A simple initial estimation of the typical velocity along the line of sight (Δ_V) is given by the velocity where the spectrum intensity is maximum. As the three other unknown parameters (T_ex, N, and σ_V) are highly correlated (see Sect. 4.6), we systematically search in a 3D grid defined as follows.

• We sample N logarithmically between 10¹² and 10¹⁸ cm⁻² with a step of 0.1 (i.e., 61 values).

• We sample T_ex logarithmically between 3 and 100 K with a step of 0.05 (i.e., 31 values).

• Finally, we first sample σ_V with 0.2, 0.3,..., 0.6, 1.2,..., 3.8 km s⁻¹ (i.e., 10 values) before refining the search with a step of 0.025 km s⁻¹ for σ_V ≤ 0.6 km s⁻¹, and of 0.05 km s⁻¹, otherwise.

These values have been fixed empirically based on simulations for which we tried to find a tradeoff between accuracy and computation time.

B.3 Maximization of the likelihood function through a scoring algorithm

The likelihood function is maximized here using Fisher’s scoring algorithm (Garthwaite et al. 1995). It is an iterative algorithm, (B.4)

where i is the ith iteration, p_i is a constant, I_F(θ) is the Fisher matrix (see Sect. 4) seen as a function of θ, and ∇_θ (θ) is the gradient, (B.5)

where j is the index of the unknown parameter.

In practice, at each iteration i, the algorithm tries p_i = 0.1, 0.4 and 0.8, and it makes a quadratic fit to get an estimation of p_i that minimizes the log-likelihood in the interval [0.1, 0.8]. Moreover, the inversion of I_F is made with the pseudo-inverse when I_F becomes singular (in the sense that the ratio between the largest and the smallest eigenvalues of I_F is larger than 10⁸). Finally, the iteration loop stops when the log-likelihood verifies or when the number of iterations reaches 1000.

B.4 Generalization to two velocity components

When two velocity components are needed, the log-likelihood of Eq. (B.3) becomes: (B.6)

where s_n,l(θ_m) is the spectrum corresponding to the component m ∈{1, 2}. The grid used in the initialization step now has six dimensions: . One solution is to consider the total Fisher matrix of size 8 × 8 (see Sect. 4.2), but it can lead to singularities. We empirically observe that iteratively applying the gradient to each component separately (θ₁, and then θ₂) actually leads to better results than using a gradient on the enlarged vector . Such a coordinate descent only requires inversions of 4 × 4 Fisher matrices.

B.5 Computing load and optimization

From the computational viewpoint, the estimation may be carried out many times either because it may be applied to many different lines of sight or because it may be used in Monte Carlo simulations. It is thus useful to actually compute in advance a 5D set of per line. In this 5D set, T_ex, N and σ_V are sampled as described in Appendix B.2. Moreover, we consider ten different values of Δ_V with a step of 0.05 km s⁻¹ and we restrict the range of explored channels to the velocity range where the lines appear, that is, an interval of 26.5 km s⁻¹ around the initially estimated Δ_V. This last point substantially decreases the computation time. For a single species, the computation of the 5D set of takes around 13 s in our Matlab implementation on a standard 2016 laptop. Each subsequent estimation (initialization plus gradient) takes around 0.05 or 1.0 second when estimating one or two velocity components, respectively.

For the considered ORION-B data, along with the initialization proposed in Appendix B.2, the median number of iterations required to reach convergence is 40 or 200 when estimating one or two velocity components, respectively.

Appendix C: Performance of the maximum likelihood estimator

Monte Carlosimulations with P independent realizations of the estimator are used toanalyze its performance. We simulate data for the ¹³CO(1− 0) and ¹³CO(2− 1) lines using the different values of T_ex and N already used in Sect. 4.3. We here compute P = 200 simulations with random noise for each pair of (T_ex, N) values.

After showing that MLE estimates of the parameters can be injected in the computation of the CRB to detect accurate estimations, we show that the proposed estimator is unbiased and efficient.

C.1 Detection of (in)accurate estimations

In order to physically interpret the estimations computed on observed data, it is crucial to remove inaccurate estimations. We thus need a way to quickly detect inaccurate estimations. Figure C.1 shows the fraction of the Monte Carlo realizations for each pair of (T_ex, N) values that deliver “accurate” estimations of all parameters. Here, “accurate” indicates that all the following conditions are simultaneously satisfied: (C.1)

In these equations, ρ is the relative precision required for all the parameters. At first sight, Eqs. (C.1) and (25) seem identical. However, we here use estimation of the parameters , while we used values of θ used to simulate the data in Sect. 4.7. The value R(ρ) is the fraction of “accurate” estimations detected without a priori information on the parameters. The contours in Fig. C.1 correspond to the frontiers where the relative precision on the actual values θ of the four parameters is ρ = 20 or 10%. Figure C.1 thus clearly suggests that these frontiers are close to the pixels for which R(ρ) ≃ 0.5. The test R(ρ) ≥ 0.5 thus gives a fair detection of accurate (and inaccurate) estimations. Furthermore, these maps show that, most of the time, all P estimations are either inaccurate (R(ρ) = 0 in blue) or accurate (R(ρ) = 1, in yellow). In other words, a single estimation is often sufficient to detect whether it is accurate or not.

Fig. C.1 Variations of the relative number of accurate estimations as a function of the column density and excitation temperature. The contour corresponds to the frontier where the relative precision on the actual values of the four parameters is ρ = 20% (solid red contour on the left) and 10% (dashed red contour on the right). The 13CO(1−0) and 13CO(2− 1) lines are simulated with σb = 100 mK, ΔV = 1.1 km s−1, and σV = 0.61 km s−1.

Fig. C.2 Variations of the bias (left) and efficiency (right) of the maximum likelihood estimator as a function of the column density and excitation temperature. The solid and dashed red contours correspondto the frontiers where the relative precisions on the estimations on the four parameters are ρ = 20, and 10%, respectively, as defined in Fig. C.1. The 13CO(1−0) and 13CO(2− 1) lines are simulated with σb = 100 mK, ΔV = 1.1 km s−1, and σV = 0.61 km s−1.

More precisely, the light blue area in Fig. C.1.a corresponds to a value R(0.2) = 0.15. This means that we still have a 15% chance of considering an estimation as accurate when it is in fact inaccurate, when N ≃ 10¹⁵ cm⁻² and T_ex > 50 K. It is possible to improve this on observed data because adjacent pixels on the sky have physical parameters that are partially correlated. We can thus assume that accurate and inaccurate estimations are spatially grouped, and the detection can be improved by computing Eq. (C.1) in a sliding window of size 3 × 3 pixels. In practice, we remove estimations of pixels for which one of the neighbors in the 3 × 3 pattern violates Eq. (C.1).

C.2 Bias and variance of the estimator

The bias and variance of the maximum likelihood estimator can be estimated with (Garthwaite et al. 1995): (C.2)

where P is the number of simulations in the Monte Carlo analysis, θ_i are the actual values of the parameters, and are the estimated values. Figure C.2 shows the variations of the ratios and as a function of T_ex and N.

It is evident that the bias of the proposed estimator is negligible compared to its standard deviation (i.e., ), and that its variance reaches the Cramer Rao bound (i.e., ) in the region where the CRBs are sufficiently small to get accurate estimations.

Appendix D: Additional figures

In Sect. 5.5, Fig. 12 shows the residuals for only one line. Figure D.1 is a generalization with regard to the other lines.

Fig. D.1 Spatial variations of the observed spectrum “energy” (first column), of the residual “energy” (second column), of their ratio (third column), and of the ratio of “energy” that has not been modeled. The unit of the color look-up table is Kelvin2 or % depending on the column. White contours show the regions where two components have been detected.

References

All Tables

All Figures

	Fig. 1 Spatial distributions of the integrated intensity (in K km s−1) of the considered lines. The maps have been rotated counterclockwise by 14 degrees from the RA/Dec J2000 reference frame. The spatial offsets are given in arcsecond from the projection center located at 05h 40m 54.270s, −02°28′00.00′′. Red crosses stand for two particular lines of sight, which are analyzed in Fig. 10.
In the text

Fig. 2 Variations of the square root of the Cramer-Rao bound (CRB) of the centroid velocity (ΔV, left panel) and velocity dispersion (σV, middle panel) in km s−1 as afunction of the column density and the excitation temperature. Right panel: variations of a function of the product of the number of channels (K) and the signal-to-noise ratio . In all three cases, the data are simulated assuming that 13CO(1−0) and 13CO(2− 1) are measured and the unit of the image contours are km s−1. In this simulation, σb = 100 mK, ΔV = 1.1 km s−1 and σV = 0.61 km s−1.

In the text

Fig. 3 Top: variations of the square root of the CRB of Tex in Kelvin as a function of the column density and the excitation temperature. Bottom: functions of the line opacities. Left: only 13CO (1−0) is analyzed. Middle: only 13CO(2−1) is analyzed. Right: both 13CO(1−0) and 13CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to Tex and N values whichlead to singular Fisher matrices. For this analysis, σb = 100 mK, ΔV = 1.1 km s−1 and σV = 0.61 km s−1. The constant σV,0 = 1 km s−1 is introduced to have expressions in (d–f) that depend on σV, but remain homogeneous to a temperature.

In the text

Fig. 4 Top row: variations of the square root of the CRB of logN as a function of the column density and the excitation temperature. Second and third row: correlation coefficients between efficient estimators of (Tex, logN), and (logN, σV) in the second and third rows. respectively (defined in Eqs. (23) and (24)). Bottom row: variations of functions of the opacities. Left: only 13CO(1− 0) is analyzed. Middle: only 13CO(2−1) is analyzed. Right: both 13CO(1−0) and 13CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to Tex and N values whichlead to singular Fisher matrices. For this analysis, σb = 100 mK, ΔV = 1.1 km s−1 and σV = 0.61 km s−1.

In the text

Fig. 5 Illustration of the correlation between N and σV estimations when a single line (13CO(1−0)) is available. The blue points in the scatter plots show the estimations of N and σV obtained with a Monte Carlo simulation of individual spectra that share the same physical parameters and different realizations of a white Gaussian noise with standard deviation σb = 100 mK. The parameters are Tex = 18 K, N = 1017.5 cm−2, ΔV = 1.1 km s−1, and σV = 0.61 km s−1.

In the text

Fig. 6 Noise standard deviation σb,ρ in mK whichensures that relative precisions are better than ρ% (for details see Eq. (25)). Top: analysis of a single line. Bottom: analysis of two lines 13CO(1− 0) and 13CO(2− 1). The contours for 10 and 100 mK are highlighted because these σb values bracket the values reached during typical observations at the IRAM-30 m. For this analysis, ΔV is fixed at 1.1km s−1, but the computations are done for four different values of σV (0.3, 0.6, 1.3,and 2.0 km s−1) and then projected on the (N, Tex) plane (see text for details).

In the text

	Fig. 7 Same as Fig. 6, except that the contours show the variations of the peak-signal-to-noise ratio, , required to reach a given relative accuracy (ρ%).
In the text

	Fig. 8 Noise standard deviation σb,ρ in mK which ensures relative precisions better than 20% when a single line of 12CO is analyzed. Other details are identical to Fig. 6.
In the text

	Fig. 9 Noise standard deviation σb,ρ in mK which ensures relative precisions better than 20% when a couple of 12CO lines are observed: J = 1−0 and a higher J line. Other details are identical to Fig. 6.
In the text

	Fig. 10 Two examples of LTE fit of the CO isotopologues lines. In the titles, Tex is expressed in Kelvin, N in cm−2, ΔV in km s−1, and σV in km s−1. The plain lines show the data, and the dotted ones show the fit results. Values in red indicate estimations whose relative precision is larger than 20%. The associated lines of sight can be localized in Fig. 1 (see red crosses).
In the text

	Fig. 11 Jointhistograms of the standard deviations of the residuals (σr) and of the noise (σb) for all the studied lines. Standard deviations are expressed in Kelvin. The dashed lines correspond to ratios 1 and 10.
In the text

	Fig. 12 Spatial variations of the observed spectrum “energy” (a), of the residual “energy” (b), of their ratio in % (c), and of the ratio of “energy” in % that has not been modeled (d). The unit of the color look-up table is Kelvin2. White contours show the regions where two components have been detected.
In the text

Fig. 13 Top: spatial variations of the estimated physical parameters. From left to right: centroid velocity, velocity dispersion, excitation temperature, column density, and line opacities. Bottom: spatial variations of the relative precisions. From left to right: relative precision on the centroid velocity, velocity dispersion, excitation temperature, column density, and line peak S/N. The black contours on the peak-signal-to-noise-ratio image delimit the regions where . On all images, red contours delimit the regions where the relative precision is better than 20% for all estimated parameters, and white contours delimit regions where two components have been estimated. In this lattercase, the images only show the estimation that is the closest (in terms of centroid velocity ΔV) to its neighboring pixels.

In the text

	Fig. 14 Comparison between the map of the dust temperature and maps of temperature ratios. Inaccurate estimations are filtered out. The dust temperature is only presented in regions with an accurate estimation of parameters.
In the text

Fig. 15 Scatter plots between the CO isotopologue excitation temperatures. The color scale encodes the dust temperature Tdust. The ellipsesrepresent the interval of confidence for each estimation: Each ellipse is centered on the estimation of the excitation temperature for one pixel and its horizontal and vertical sizes are equal to the associated CRBs. Dashed ellipses correspond to pixels with two components. Inaccurate estimations are filtered out. Dashed red lines show the loci of ratios 1/2, 1, 2, and 4.

In the text

	Fig. 16 Plots of the 13CO column densities per unit intensity of 13CO(1−0) (top) and 13CO (2−1) (bottom) as a function of the excitation temperature. This plot is done for four values of the 13CO column density, and σV = 0.61 km s−1.
In the text

	Fig. 17 Comparison of the map of dust-traced H2 column density with maps of various column density ratios. Inaccurate estimations were filtered out. The dust-traced column density is only presented in regions with an accurate estimation of parameters.
In the text

	Fig. 18 Scatter plot of the CO isotopologue column densities as a function of the dust-traced H2 column density. Inaccurate estimations are filtered out. The black dashed lines show the expected gas phase abundances relative to H2 with no depletion: C18O/H2 = 5 × 10−7, 13CO/H2 = 4 × 10−6, and CO/H2 = 1.4 × 10−4.
In the text

Fig. 19 Scatter plots between CO isotopologue column densities. The color scale encodes the ratio of the excitation temperatures of the considered species in each panel. The ellipses represent the interval of confidence for each estimation: Each ellipseis centered on the estimation of the column density for one pixel and its horizontal and vertical sizes are equal to the associated CRBs. Dashed ellipses correspond to pixels with two components. Inaccurate estimations are filtered out. The dashed red lines show the loci of ratios 5, 20, and 60.

In the text

	Fig. 20 Joint histogram of the logarithm of the estimated column density ratios and of the observed integrated intensity ratios. All the ratios are computed for the 13CO(1− 0) line over the C18O(1−0) one. The dashed red lines correspond respectively to ratios of one and two.
In the text

Fig. C.1 Variations of the relative number of accurate estimations as a function of the column density and excitation temperature. The contour corresponds to the frontier where the relative precision on the actual values of the four parameters is ρ = 20% (solid red contour on the left) and 10% (dashed red contour on the right). The 13CO(1−0) and 13CO(2− 1) lines are simulated with σb = 100 mK, ΔV = 1.1 km s−1, and σV = 0.61 km s−1.

In the text

Fig. C.2 Variations of the bias (left) and efficiency (right) of the maximum likelihood estimator as a function of the column density and excitation temperature. The solid and dashed red contours correspondto the frontiers where the relative precisions on the estimations on the four parameters are ρ = 20, and 10%, respectively, as defined in Fig. C.1. The 13CO(1−0) and 13CO(2− 1) lines are simulated with σb = 100 mK, ΔV = 1.1 km s−1, and σV = 0.61 km s−1.

In the text

	Fig. D.1 Spatial variations of the observed spectrum “energy” (first column), of the residual “energy” (second column), of their ratio (third column), and of the ratio of “energy” that has not been modeled. The unit of the color look-up table is Kelvin2 or % depending on the column. White contours show the regions where two components have been detected.
In the text