Issue 
A&A
Volume 645, January 2021



Article Number  A26  
Number of page(s)  27  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202037776  
Published online  23 December 2020 
C^{18}O, ^{13}CO, and ^{12}CO abundances and excitation temperatures in the Orion B molecular cloud
Analysis of the achievable precision in modeling spectral lines within the approximation of the local thermodynamic equilibrium
^{1}
Aix Marseille Univ., CNRS, Centrale Marseille, Institut Fresnel,
Marseille,
France
email: antoine.roueff@fresnel.fr
^{2}
LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités,
75014
Paris, France
^{3}
Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allee Geoffroy SaintHilaire,
33615
Pessac, France
^{4}
Laboratoire de Physique de l’Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, Sorbonne Paris Cité,
Paris, France
^{5}
IRAM, 300 rue de la Piscine,
38406
Saint Martin d’Hères,
France
^{6}
Instituto de Física Fundamental (CSIC). Calle Serrano 121,
28006
Madrid, Spain
^{7}
Chalmers University of Technology, Department of Space,
Earth and Environment,
412 93 Gothenburg, Sweden
^{8}
University of Toulouse, IRIT/INPENSEEIHT, CNRS,
2 rue Charles Camichel, BP 7122,
31071
Toulouse Cedex 7, France
^{9}
LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités,
92190
Meudon, France
^{10}
Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, GIPSALab,
38000
Grenoble, France
^{11}
Univ. Lille, CNRS, Centrale Lille, UMR 9189 – CRIStAL,
59651
Villeneuve d’Ascq, France
^{13}
Instituto de Astrofísica, Pontificia Universidad Católica de Chile,
Av. Vicuña Mackenna 4860,
7820436
Macul,
Santiago, Chile
^{13}
Institut de Recherche en Astrophysique et Planétologie (IRAP), Université Paul Sabatier,
Toulouse Cedex 4, France
^{14}
National Radio Astronomy Observatory,
520 Edgemont Road,
Charlottesville,
VA
22903, USA
^{15}
Canadian Institute for Theoretical Astrophysics, University of Toronto,
60 Saint George Street, 14th floor,
Toronto,
ON,
M5S 3H8, Canada
^{16}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot,
Sorbonne Paris Cité,
91191
GifsurYvette,
France
^{17}
School of Physics and Astronomy, Cardiff University,
Queen’s buildings,
Cardiff
CF24 3AA, UK
Received:
19
February
2020
Accepted:
28
April
2020
Context. CO isotopologue transitions are routinely observed in molecular clouds for the purpose of probing the column density of the gas and the elemental ratios of carbon and oxygen, in addition to tracing the kinematics of the environment.
Aims. Our study is aimed at estimating the abundances, excitation temperatures, velocity field, and velocity dispersions of the three main CO isotopologues towards a subset of the Orion B molecular cloud, which includes IC 434, NGC 2023, and the Horsehead pillar.
Methods. We used the Cramer Rao bound (CRB) technique to analyze and estimate the precision of the physical parameters in the framework of localthermodynamicequilibrium (LTE) excitation and radiative transfer with added white Gaussian noise. We propose a maximum likelihood estimator to infer the physical conditions from the 1–0 and 2–1 transitions of CO isotopologues. Simulations show that this estimator is unbiased and proves efficient for a common range of excitation temperatures and column densities (T_{ex} > 6 K, N > 10^{14}−10^{15} cm^{−2}).
Results. Contrary to general assumptions, the various CO isotopologues have distinct excitation temperatures and the line intensity ratios between different isotopologues do not accurately reflect the column density ratios. We find mean fractional abundances that are consistent with previous determinations towards other molecular clouds. However, significant local deviations are inferred, not only in regions exposed to the UV radiation field, but also in shielded regions. These deviations result from the competition between selective photodissociation, chemical fractionation, and depletion on grain surfaces. We observe that the velocity dispersion of the C^{18}O emission is 10% smaller than that of ^{13}CO. The substantial gain resulting from the simultaneous analysis of two different rotational transitions of the same species is rigorously quantified.
Conclusions. The CRB technique is a promising avenue for analyzing the estimation of physical parameters from the fit of spectral lines. Future works will generalize its application to nonLTE excitation and radiative transfer methods.
Key words: ISM: molecules / ISM: clouds / radiative transfer / methods: data analysis / methods: statistical
© A. Roueff et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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^{2}−10^{5} cm^{−3}, Draine 2011) are wellsuited for the emission in the lowenergy 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 signaltonoise ratios (S/N) over large fields of view. The Outstanding RadioImaging of OrioN B (ORIONB IRAM; coPIs: J. Pety and M. Gerin) 30m 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 ^{12}CO(1− 0), ^{13}CO(1− 0), and C^{18}O (1−0) lines are indeed tracing the molecular gas well. Their quantitative comparisons show that the H_{2} column density deduced from the dust emission can be accurately estimated from the ^{12}CO(1− 0), ^{13}CO(1− 0), and C^{18}O (1−0) lines in the column density range from 10^{21} to ≳ 10^{22} cm^{−2}.
Assuming identical excitation temperatures, the opacity of the ground state transitions is expected to be smaller for ^{13}CO than for ^{12}CO, and even smaller for C^{18}O, because of the difference in elemental abundances, ^{12}C∕^{13}C ~ 60 and ^{16}O ∕^{18}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 ^{13}C^{+} and ^{12}CO leads to an enhancement of the ^{13}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 ^{13}CO emitting region and brings it closer to that of ^{12}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 ^{13}C elemental abundance from CO observations only. Up to now, the most reliable determinations of the ^{12}C∕^{13}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^{18}O and the doubly isotopic species ^{13}C^{18}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 photodissociation regions (e.g., Visser et al. 2009). In observations, it is clearly seen as an offset between the threshold for the apparition of ^{12}CO (near A_{V} = 0.5 mag) and C^{18}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 ^{13}CO abundance is enhanced through fractionation in the same regions where the C^{18}O abundance decreases due to selective photodissociation. This leads to a broad range of the ^{13}CO/C^{18}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^{18}O reaches a minimum value in the filament ridges and that it is always lower than the velocity dispersion determined by ^{13}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 ORIONB 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, ^{12}CO, ^{13}CO, and C^{18}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 deexcitation. 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_{2}, the low energy rotational lines of carbon monoxide are bright and easily thermalized in collisions with H_{2}, 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 nonlocal, nonLTE radiative transfer models of a uniform (constant density and temperature) spherical cloud, Bernes (1979) has shown that the excitation temperatures of the ^{12}CO(1− 0) and ^{12}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 farUV 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 ^{12}CO, ^{13}CO, and C^{18}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 ^{13}CO(1−0), ^{13}CO(2− 1), C^{18}O (1−0), C^{18}O (2−1) and ^{12}CO(1− 0) lines towards parts of the Orion B molecular cloud. We compared our results with the dusttraced H_{2} column density and dust temperature. This section describes the associated data sets.
2.1 IRAM30 m observations
2.1.1 3 mm CO lines from the ORIONB large program
The 3 mm data were obtained with the IRAM30 m as part of the ORIONB 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 IRAM30m 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 intercalibrated. The absolute flux calibration at 3 mm for the IRAM30 m telescope is estimated to be better than 5%.
2.1.2 1 mm CO lines
The ^{13}CO(2−1) and C^{18}O (2−1) lines were also observed at the IRAM30 m in 2006 (PI: N. Peretto) using the ABCD generation of receivers and the VESPA autocorrelator. The two lines were observed simultaneously ensuring an excellent intercalibration.
Data reduction was carried out using the GILDAS^{1}/CLASS software. The contribution of the atmosphere was first removed (ON–OFF procedure) and the data were calibrated to the scale using the standard chopperwheel method (Penzias & Burrus 1973). The data were then converted to mainbeam temperatures using the standard forward (0.94) and mainbeam (0.62) efficiencies for the ABCD receiver around 220 GHz^{2}. The resulting absolute flux calibration is estimated to be better than 10%. We subtracted a firstorder baseline from every spectrum, excluding the velocity range from 5 to 15 km s^{−1} 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 IRAM30 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^{4} τ_{850} mag, and the visual extinction into H_{2} column density using N(H_{2})∕A_{V} = 0.9 × 10^{21}H cm^{−2} mag^{−1}.
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^{h}40^{m}54.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 IRAM30 m angular resolution ranges from 11.5″ at 220 GHz to 23.5″ at 110 GHz. The positionpositionvelocity 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 positionpositionvelocity cubes^{3} of 129 × 170 × 80 pixels, each pixel covering 9″ × 9″ × 0.5 km s^{−1} (Nyquist sampling at 3 mm). Figure 1 shows the maps of the intensity integrated between 0 and 20 km s^{−1} 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 05^{h} 40^{m} 54.270^{s}, −02°28′00.00^{′′}. Red crosses stand for two particular lines of sight, which are analyzed in Fig. 10. 
Properties of the observed lines.
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 multiresolution 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 wellknown (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 ONOFF calibration and spectrum baselining to subtract the slowlyvarying continuum residual from the receiver and the atmosphere) delivers a noise b that is centered (i.e., with zeromean) 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 lowvelocity 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 CramerRao bound analysis
In this paper, we aim to estimate the physical parameters of the LTE model presented in Sect. 3 based on the ORIONB 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 CramerRao 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 (^{13}CO(1−0)) or two (^{13}CO(1−0) and ^{13}CO(2− 1)) lines are available. Second, we consider which ^{12}CO line should be observed to improve the precision achieved when only ^{12}CO(1− 0) observations are available.
Fig. 2 Variations of the square root of the CramerRao 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 signaltonoise ratio . In all three cases, the data are simulated assuming that ^{13}CO(1−0) and ^{13}CO(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^{4} 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 ORIONB project. To generate Figs. 2–4, we assume that the two measured lines are ^{13}CO(1− 0) and ^{13}CO(2− 1). The two corresponding opacities are noted τ_{1} and τ_{2}. The number of samples is K = 80 for each line at a spectral resolution of 0.5 km s^{−1}. 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 CramerRao 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^{13} cm^{−2} and 10^{19} cm^{−2}, and the other two parameters are kept constant at Δ_{V} = 1.1 km s^{−1} and σ_{V} = 0.61 km s^{−1} (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 T_{ex} in Kelvin as a function of the column density and the excitation temperature. Bottom: functions of the line opacities. Left: only ^{13}CO (1−0) is analyzed. Middle: only ^{13}CO(2−1) is analyzed. Right: both ^{13}CO(1−0) and ^{13}CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to T_{ex} 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^{16} cm^{−2} and T_{ex} ≥ 12 K, the square root of both CRBs are smaller than 0.01 km s^{−1}. 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 signaltonoise 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^{−1}. 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 ^{13}CO(1−0) or ^{13}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^{16} cm^{−2}. 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 ^{13}CO(1− 0) line, it shows a different behaviour for the ^{13}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 ^{13}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^{16}−10^{17.5}] cm^{−2}. 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^{14}−10^{16.5}] cm^{−2}. 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 ^{13}CO(1−0) or ^{13}CO(2− 1). In this equation, σ_{V,0} = 1 km s^{−1} 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 (T_{ex}, 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 ^{13}CO(1− 0) is analyzed. Middle: only ^{13}CO(2−1) is analyzed. Right: both ^{13}CO(1−0) and ^{13}CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to T_{ex} 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 nondiagonal 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) anticorrelated. Starting from a modeled spectrum with log (N cm^{−2}) = 17.5, σ_{V} = 0.61 km s^{−1}, 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 anticorrelated. However, the standard deviation on the log (N cm^{−2}) and σ_{V} estimations are 0.017 and 0.005 km s^{−1}, 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^{16}−10^{17}] cm^{−2} (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^{15}−10^{17}] cm^{−2} (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^{17} cm^{−2}, 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 τ_{1} and τ_{2} are close to 1. With two lines (see Fig. 4l), the best precision is when τ_{1} < 1 and τ_{2} > 1.
Fig. 5 Illustration of the correlation between N and σ_{V} estimations when a single line (^{13}CO(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 T_{ex} = 18 K, N = 10^{17.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 ^{13}CO(1− 0) and ^{13}CO(2− 1). The contours for 10 and 100 mK are highlighted because these σ_{b} values bracket the values reached during typical observations at the IRAM30 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, T_{ex}) 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^{5} 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^{−1}. 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 IRAM30 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 ^{13}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 peaksignaltonoise ratio defined for one transition, l, by: (26)
and for two transitions by . Figure 7 shows the variations of minimum peaksignaltonoise ratio for similar conditions as in Fig. 6. Analyzing only the ^{13}CO(1−0) line requires at least a S/N of 100 to get a relative precision of 20%. Adding the ^{13}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 peaksignaltonoise ratio, , required to reach a given relative accuracy (ρ%). 
4.8 Complementing ^{12}CO(1–0) observations
All the previous analyses were done for the ^{13}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 ^{12}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 ^{12}CO. A global pattern is seen, especially for the higher energy transitions ^{12}CO(4−3), ^{12}CO(5− 4), ^{12}CO(6− 5), For instance, if N lies in the interval [10^{17}, 10^{19}] cm^{−2} and T_{ex} > 24 K, the ^{12}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 groundbased 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 ^{12}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^{−2}. 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 ^{12}CO transitions (van der Tak et al. 2007) because of the higher critical densities of the higherJ transitions. Nevertheless, the usefulness of mildly excited lines remains valid. NonLTE approaches will be developed in the future that will provide a quantitative assessment of the diagnostic power of these lines.
5 Application to the ORIONB 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 ORIONB 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 ^{12}CO 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 ^{12}CO 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 ^{13}CO(1− 0), C^{18}O(1−0), and ^{12}CO(1− 0) lines. If one of the denoised spectra for ^{13}CO(1−0), C^{18}O(1−0) or ^{12}CO(1− 0) has at least two local maxima in the velocity interval of interest [8.25, 14.25] km s^{−1}, we then fix the number of velocity components to two for all three species. The denoised spectra ^{13}CO(1− 0), C^{18}O(1−0) and ^{12}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 ^{13}CO(1− 0) line has both a good S/N and moderate opacities. The C^{18}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 ^{12}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 ^{13}CO and C^{18}O. The maximization of the loglikelihood 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 ^{12}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 ^{13}CO(1−0), ^{13}CO(2− 1), C^{18}O(1−0), C^{18}O(2−1), and we use the estimation of σ_{V} obtained on ^{13}CO to fix σ_{V} for the estimation of the other parameters (T_{ex}, log N and Δ_{V}) in the analysis of the ^{12}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^{18}O, ^{13}CO, and ^{12}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 ^{12}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 wellfitted with two different velocity components. The estimations for the C^{18}O and ^{13}CO species are physically relevant because the two velocity components are well separated in velocity. This is less obvious for ^{12}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

The “energy” of the estimated signal is
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^{−1} 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 underfitted (positive value) or overfitted (negative value). Measuring the residual “energies” is similar to computing a χ^{2} in leastsquare 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, T_{ex} 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^{18}O lines are better fitted than the ^{13}CO lines, and that the ^{12}CO lines are the least wellfitted. 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 misspecification in the model. The best fit happens for the C^{18}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 ^{12}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 ^{12}CO(1− 0) line. Another limit of the model is that it does not encode selfabsorption 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 ^{13}CO(1−0) line. Both images share the same color lookup 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 underfits 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. NonLTE 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 lookup table is Kelvin^{2}. White contours show the regions where two components have been detected. 
Statistics of the estimated parameters over all the pixels for the three CO isotopologues.
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^{18}O, ^{13}CO, and ^{12}CO
Figure 13 shows the spatial variations of the estimated parameters and associated uncertainties for the C^{18}O, ^{13}CO, and ^{12}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 peaksignaltonoise ratio of the ^{12}CO(1−0), ^{13}CO(1− 0), and ^{13}CO(2− 1) lines is large (> 20) over most of the studied field of view. The regions where the peaksignaltonoise ratio of the C^{18}O(1−0) and C^{18}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 ^{12}CO S/N is much larger than the C^{18}O or the ^{13}CO one, the ^{12}CO estimations are more uncertain (see map of for ^{12}CO in Fig. 13) because a single transition is available and the ^{12}CO opacities are large (ranging from 5 to about 1000). However, the higher S/N for ^{12}CO helps to derive the velocity field in regions where ^{13}CO and C^{18}O are not well detected (diffuse gas).
The largest variations are observed in the column density which varies from the detection limit near 10^{15} cm^{−2} up to values larger than 10^{17} cm^{−2} for ^{13}CO. The comparison of the column density maps for ^{13}CO and C^{18}O suggests that the C^{18}O molecules are confined to the high column density regions and avoid the cloud edges. The ^{12}CO isotopologue shows a different behavior with emission extending over most of the imaged field of view and column densities ranging from ~10^{16} to ~ 10^{19} cm^{−2}. These behaviors are related to the difference of opacity of the isotopologues lines. Both ^{13}CO lines have moderate opacities across the mapped region, with somewhat higher values for ^{13}CO(2− 1) than for ^{13}CO(1− 0). The ^{12}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^{18}O emission is fainter and less extended than that of ^{13}CO, both speciesprovide similar kinematics information. The velocity field is spatially regular with wellresolved 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^{−1}, with values around 0.3 km s^{−1} in the Horsehead nebula (see also HilyBlant et al. 2005), and somewhat narrower lines for C^{18}O than for ^{13}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 peaksignaltonoiseratio 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 (^{12}CO) has optically thick lines while the rarer isotopologues (^{13}CO and C^{18}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^{5} cm^{−3} 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 ^{12}CO isotopologue has the largest excitation temperature, followed by ^{13}CO and C^{18}O. Table 3 lists the typical ratios of excitation temperatures. They are COCO) ~ 1.7 and CO)∕T_{ex}(C^{18}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 diffusetranslucent cloud along the line of sight to HD 62542.
The excitation temperature of C^{18}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 ^{13}CO(1–0) (mean value T_{ex} ∕T_{dust} = 0.76), while most of the positions show ^{12}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 UVdominated 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 ^{12}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^{18}O and ^{13}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 columntointensity 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 columntointensity ratio for the ^{13}CO(1−0) and ^{13}CO(2− 1) lines for different column densities of ^{13}CO, assuming a velocity dispersion of 0.61 km s^{−1}. Column densities of correspond to optically thin lines, while the opacity becomes significant (i.e., τ ≥ 0.5) for 10^{16} −10^{17} cm^{−2}. In the optically thin case, the columntointensity 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 columntointensity ratio curve changes andthe minimum shifts to higher excitation temperatures or possibly disappears. Therefore, using the ^{12}CO excitation temperature for determining the column densities of ^{13}CO and C^{18}O leads to errors in the estimation of these column densities because the associated columntointensity ratio is inappropriate. For moderate column densities ( cm^{−2}), 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 columntointensity curve. A CO excitation temperature near the minimum of the columntointensity curve minimizes the error because a small difference in T_{ex} in this region does not change the columntointensity. For high column densities (N ~ 10^{17} cm^{−2}), the bias is significant at low excitation temperatures (T_{ex} < 20 K) because the columntointensity ratio (purple curve in Fig. 16) is rising fast at low excitation temperatures.
Statistics of ratios of the estimated parameters over all the pixels for the three CO isotopologues.
Fig. 15 Scatter plots between the CO isotopologue excitation temperatures. The color scale encodes the dust temperature T_{dust}. 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_{2} (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, ^{13}CO and ^{12}CO can be fitted and analyzed for H_{2} column densities larger than 10^{21.5} cm^{−2}. The threshold for C^{18}O is about twice higher at ~10^{21.8} cm^{−2}. Indeed, as shown by Pety et al. (2017) and Orkisz et al. (2019), the threshold for the detection of C^{18}O is A_{V} ~ 3 mag or N(H_{2}) = 10^{21.5} cm^{−2} in the Orion B molecular cloud, while the thresholds for ^{12}CO and C^{13}O are close to A_{V} = 1 mag.
Over the mapped area, the mean abundances are well defined at CO)∕N(H_{2}) = 10^{−5.6 ± 0.29} = 2.5 ± 1.5 × 10^{−6}, and N(C^{18}O)∕N(H_{2}) = 10^{−6.7±0.15} = 2.0 ± 0.8 × 10^{−7} (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^{−4}, ^{12}C/^{13}C = 57−67 and ^{16}O/^{18}O = 500−560 (Gerin et al. 2015; Langer & Penzias 1990; Wilson & Rood 1994), namely CO)∕N(H_{2}) = 4−5 × 10^{−6}, and N(C^{18}O)∕N(H_{2}) = 5−6 × 10^{−7}. The difference is more pronounced for C^{18}O because its abundance is affected by both photodissociation and freezeout over a more significant fraction of the studied area.
Significant deviations from the mean values are present. The C^{18}O abundance is not only lower near the photoilluminated 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^{18}O and ^{13}CO supporting the explanation by a global freezeout effect.
Because the ^{12}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 ^{12}CO to the one estimated for ^{13}CO (see Sect. 5.2) allows us to determine the ^{12}CO column density for the pixels that have the least saturated line emission. Although fairly uncertain, the resulting values of the ^{12}CO abundance relative to H_{2} are close to 6 × 10^{−5}, lower than the expected value using the carbon elemental abundance relative to H, 1.4 × 10^{−4} (Gerin et al. 2015), but similar to values obtained in Taurus by Pineda et al. (2010).
Fig. 16 Plots of the ^{13}CO column densities per unit intensity of ^{13}CO(1−0) (top) and ^{13}CO (2−1) (bottom) as a function of the excitation temperature. This plot is done for four values of the ^{13}CO column density, and σ_{V} = 0.61 km s^{−1}. 
Fig. 17 Comparison of the map of dusttraced H_{2} column density with maps of various column density ratios. Inaccurate estimations were filtered out. The dusttraced 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 dusttraced H_{2} column density. Inaccurate estimations are filtered out. The black dashed lines show the expected gas phase abundances relative to H_{2} with no depletion: C^{18}O/H_{2} = 5 × 10^{−7}, ^{13}CO/H_{2} = 4 × 10^{−6}, and CO/H_{2} = 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 ^{12}CO/^{13}CO = 57−67 and ^{13}CO/C^{18}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 ^{13}CO and C^{18}O column densities is indeed in the expected range at CO)∕N(C^{18}O) = 8. The upper bound is close to 50 indicating that chemical effects play a significant role, by enhancing the ^{13}CO abundance (fractionation) or destroying C^{18}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 ^{13}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 ^{12}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 ^{12}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 ^{13}C due to the exchange reaction between ^{13}C^{+} and ^{12}CO, that enhances the ^{13}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(^{13}CO)/W(^{12}CO) versus W(^{13}CO)/W(C^{18}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 ^{12}CO and ^{13}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_{2})∕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 ^{13}CO(1− 0) line over the C^{18}O(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 ^{13}CO, C^{18}O, and ^{12}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 ^{12}CO. This implies that the actual velocity dispersion derived from high opacity lines like those of ^{12}CO is smaller than the apparent line width. Because we fixed the ^{12}CO velocity dispersion to the value obtained with ^{13}CO in our estimation, the velocity dispersions of ^{12}CO and ^{13}CO are identical (see Sect. 5.2). The velocity dispersions of C^{18}O and ^{13}CO are howeverfitted independently. Both species show similar velocity dispersions but ^{13}CO has consistently broader line profiles than C^{18}O. The ratio between the velocity dispersions of ^{13}CO and C^{18}O is 1.1 (see Table 3). The ^{13}CO emission is produced by a more extended volume along the line of sight than the C^{18}O emission as shown by the lower threshold in N(H_{2}) where ^{13}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^{18}O reaches a minimum value in the filament ridges, and is always lower than the velocity dispersion determined by ^{13}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 ^{13}CO and C^{18}O. Inspecting the spatial distribution of the velocity in Fig. 13 suggests that the small excess of velocity dispersion for ^{13}CO relative to C^{18}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(^{12}CO) ~ 6 × 10^{−5}, X(^{13}CO) = 2.5 ± 1.5 × 10^{−6}, and X(C^{18}O) = 2.0 ± 0.8 × 10^{−7}.

The excitation temperatures T_{ex} are different among the CO isotopologues. ^{12}CO presents the highest T_{ex}, followed by ^{13}CO and C^{18}O. For ^{13}CO and C^{18}O, the excitation temperatures are lower than the dust temperature on average, while they are higher for ^{12}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 ^{13}CO(1−0) and C^{18}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 ^{13}C and ^{18}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 dualband receiver systems for telescopes like the IRAM30 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 ^{12}CO, ^{13}CO, and C^{18}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 ^{12}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 preprocessing) 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 nonLTE models. Such nonLTE 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 01913, 02214, 14514, 12215, 01816, and finally the large program number 12416 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://gouldbeltherschel.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, INAFIFSI Rome and INAFArcetri, 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 ORIONB 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, cofunded by CEA and CNES. JRG thanks Spanish MICI for funding support under grant AYA201785111P. 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.3
A.4
where is numerically computed.
A.5 ∂s (ν)∕∂L_{N}
Let us introduce L_{N} = logN. (A.10)
A.6
Since N = 10^{LN} 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.9 ∂Ψ(ν)∕∂Δ_{V}
A.10 ∂s (ν) /∂σ_{V}
A.11 ∂Ψ(ν)∕∂σ_{V}
From Eq. (A.15), we get (A.21)
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 loglikelihood is: (B.3)
B.2 Initialization of the unknown parameters
The loglikelihood 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^{12} and 10^{18} cm^{−2} 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^{−1} (i.e., 10 values) before refining the search with a step of 0.025 km s^{−1} for σ_{V} ≤ 0.6 km s^{−1}, and of 0.05 km s^{−1}, 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 loglikelihood in the interval [0.1, 0.8]. Moreover, the inversion of I_{F} is made with the pseudoinverse 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^{8}). Finally, the iteration loop stops when the loglikelihood verifies or when the number of iterations reaches 1000.
B.4 Generalization to two velocity components
When two velocity components are needed, the loglikelihood 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 (θ_{1}, and then θ_{2}) 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^{−1} 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^{−1} 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 ORIONB 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 ^{13}CO(1− 0) and ^{13}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 ^{13}CO(1−0) and ^{13}CO(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 ^{13}CO(1−0) and ^{13}CO(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^{15} cm^{−2} 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 lookup table is Kelvin^{2} or % depending on the column. White contours show the regions where two components have been detected. 
References
 André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bernes, C. 1979, A&A, 73, 67 [NASA ADS] [Google Scholar]
 Barnes, P. J., Hernandez, A. K., Muller, E., & Pitts, R. L. 2018, ApJ, 866, 19 [CrossRef] [Google Scholar]
 Bonaca, A., & Hogg, D. W. 2018, ApJ, 867, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Bron, E., Daudon, C., Pety, J., et al. 2018, A&A, 610, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cernicharo, J., & Guelin, M. 1987, A&A, 176, 299 [NASA ADS] [Google Scholar]
 Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press) [Google Scholar]
 Espinosa, S., Silva, J. F., Mendez, R. A., Lobos, R., & Orchard, M. 2018, A&A, 616, A95 [CrossRef] [EDP Sciences] [Google Scholar]
 Frerking, M. A., Keene, J., Blake, G. A., & Phillips, T. G. 1989, ApJ, 344, 311 [NASA ADS] [CrossRef] [Google Scholar]
 Garthwaite, P. H., Jolliffe, I. T., & Jones, B. 1995, Statistical Inference (London: Prentice Hall Europe) [Google Scholar]
 Gerin, M., Ruaud, M., Goicoechea, J. R., et al. 2015, A&A, 573, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Goldsmith, P. F., & Langer, W. D. 1978, ApJ, 222, 881 [NASA ADS] [CrossRef] [Google Scholar]
 Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gratier, P., Pety, J., Bron, E., et al. 2020, A&A, 645, A27 [CrossRef] [EDP Sciences] [Google Scholar]
 HilyBlant, P., Teyssier, D., Philipp, S., & Güsten, R. 2005, A&A, 440, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keene, J., Schilke, P., Kooi, J., et al. 1998, ApJ, 494, L107 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, W. D., & Penzias, A. A. 1990, ApJ, 357, 477 [NASA ADS] [CrossRef] [Google Scholar]
 le Chevalier, F. 1989, Principes de Traitement des Signaux radar et Sonar (Paris, Milan, Barcelone, Mexico: Masson) [Google Scholar]
 Leung, C. M.,& Liszt, H. S. 1976, ApJ, 208, 732 [NASA ADS] [CrossRef] [Google Scholar]
 Liszt, H. S. 2006, A&A, 458, 507 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liszt, H. S. 2017, ApJ, 835, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Liszt, H. S., & Pety, J. 2012, A&A, 541, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266 [Google Scholar]
 Marchal, A., MivilleDeschênes, M.A., Orieux, F., et al. 2019, A&A, 626, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martín, S., Muller, S., Henkel, C., et al. 2019, A&A, 624, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mather, J. C., Cheng, E. S., Cottingham, D. A., et al. 1994, ApJ, 420, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Menten, K. M., Reid, M. J., Forbrich, J., & Brunthaler, A. 2007, A&A, 474, 515 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Orkisz, J. H., Peretto, N., Pety, J., et al. 2019, A&A, 624, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ossenkopf, V., Röllig, M., Neufeld, D. A., et al. 2013, A&A, 550, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Penzias, A. A., & Burrus, C. A. 1973, ARA&A, 11, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Pety, J. 2005, SF2A2005: Semaine de l’ Astrophysique Francaise (Paris: EDP Sciences Conference Series), eds. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 721 [Google Scholar]
 Pety, J., Guzmán, V. V., Orkisz, J. H., et al. 2017, A&A, 599, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Phillips, T. G., Huggins, P. J., Wannier, P. G., & Scoville, N. Z. 1979, ApJ, 231, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998, J. Quant. Spectr. Rad. Transf., 60, 883 [Google Scholar]
 Pineda, J. L., Goldsmith, P. F., Chapman, N., et al. 2010, ApJ, 721, 686 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration I. 2011, A&A, 536, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ripple, F., Heyer, M. H., Gutermuth, R., Snell, R. L., & Brunt, C. M. 2013, MNRAS, 431, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, N., André, P., Könyves, V., et al. 2013, ApJ, 766, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Sliwa, K., Wilson, C. D., Aalto, S., & Privon, G. C. 2017, ApJ, 840, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Stark, G., Heays, A. N., Lyons, J. R., et al. 2014, ApJ, 788, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Stoica, P., & Moses, R. 2005, Spectral Analysis of Signals (New Jersey: Prentice Hall) [Google Scholar]
 van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Welty, D., Sonnentrucker, P. G., Rachford, B., Snow, T., & York, D. G. 2018, AAS Meeting Abstracts, 231, 247.28 [Google Scholar]
 Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, R. W., Jefferts, K. B., & Penzias, A. A. 1970, ApJ, 161, L43 [NASA ADS] [CrossRef] [Google Scholar]
See http://www.iram.fr/IRAMFR/GILDAS for more information about the GILDAS software (Pety 2005).
For details, see http://www.iram.es/IRAMES/mainWiki/Iram30mEfficiencies
The data products associated with this paper are available at http://www.iram.fr/~pety/ORIONB
All Tables
Statistics of the estimated parameters over all the pixels for the three CO isotopologues.
Statistics of ratios of the estimated parameters over all the pixels for the three CO isotopologues.
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 05^{h} 40^{m} 54.270^{s}, −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 CramerRao 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 signaltonoise ratio . In all three cases, the data are simulated assuming that ^{13}CO(1−0) and ^{13}CO(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 T_{ex} in Kelvin as a function of the column density and the excitation temperature. Bottom: functions of the line opacities. Left: only ^{13}CO (1−0) is analyzed. Middle: only ^{13}CO(2−1) is analyzed. Right: both ^{13}CO(1−0) and ^{13}CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to T_{ex} 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 (T_{ex}, 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 ^{13}CO(1− 0) is analyzed. Middle: only ^{13}CO(2−1) is analyzed. Right: both ^{13}CO(1−0) and ^{13}CO(2− 1) are analyzed simultaneously. In (a) and (b) pixels in grey correspond to T_{ex} 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 (^{13}CO(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 T_{ex} = 18 K, N = 10^{17.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 ^{13}CO(1− 0) and ^{13}CO(2− 1). The contours for 10 and 100 mK are highlighted because these σ_{b} values bracket the values reached during typical observations at the IRAM30 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, T_{ex}) plane (see text for details). 

In the text 
Fig. 7 Same as Fig. 6, except that the contours show the variations of the peaksignaltonoise 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 ^{12}CO 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 ^{12}CO 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, T_{ex} 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 lookup table is Kelvin^{2}. 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 peaksignaltonoiseratio 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 T_{dust}. 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 ^{13}CO column densities per unit intensity of ^{13}CO(1−0) (top) and ^{13}CO (2−1) (bottom) as a function of the excitation temperature. This plot is done for four values of the ^{13}CO column density, and σ_{V} = 0.61 km s^{−1}. 

In the text 
Fig. 17 Comparison of the map of dusttraced H_{2} column density with maps of various column density ratios. Inaccurate estimations were filtered out. The dusttraced 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 dusttraced H_{2} column density. Inaccurate estimations are filtered out. The black dashed lines show the expected gas phase abundances relative to H_{2} with no depletion: C^{18}O/H_{2} = 5 × 10^{−7}, ^{13}CO/H_{2} = 4 × 10^{−6}, and CO/H_{2} = 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 ^{13}CO(1− 0) line over the C^{18}O(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 ^{13}CO(1−0) and ^{13}CO(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 ^{13}CO(1−0) and ^{13}CO(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 lookup table is Kelvin^{2} or % depending on the column. White contours show the regions where two components have been detected. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.