C18O, 13CO, and 12CO abundances and excitation temperatures in the Orion B molecular cloud: An analysis of the precision achievable when modeling spectral line within the Local Thermodynamic Equilibrium approximation

CO isotopologue transitions are routinely observed in molecular clouds to probe the column density of the gas, the elemental ratios of carbon and oxygen, and to trace the kinematics of the environment. We aim 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. We use the Cramer Rao Bound (CRB) technique to analyze and estimate the precision of the physical parameters in the framework of local-thermodynamic-equilibrium excitation and radiative transfer with an additive 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 efficient for a common range of excitation temperatures and column densities (Tex>6 K, N>1e14 - 1e15 cm-2). Contrary to the general assumptions, the different 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 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 C18O emission is 10% smaller than that of 13CO. The substantial gain resulting from the simultaneous analysis of two different rotational transitions of the same species is rigorously quantified. The CRB technique is a promising avenue for analyzing the estimation of physical parameters from the fit of spectral lines.


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 well suited for the emission in the low energy rotational lines of molecules such as carbon monoxide, which are accessible at millimeter wavelengths. The advent of sensitive broadband heterodyne receivers provides homogeneous data sets of various CO isotopologues and other species with high signal-to-noise ratios over large fields of view. The ORION-B IRAM-30m large program (Outstanding Radio-Imaging of OrioN B, co-PIs: J. Pety and M. Gerin) aims at imaging 5 square degrees towards the southern part of the Orion B molecular cloud over most of the 3 mm atmospheric window. Carbon monoxide is especially interesting 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) show 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) and Gratier et al. (submitted as a companion paper) show qualitatively and quantitatively that 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/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 on 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 involve 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 is important. It has been studied in de-tail through laboratory experiments (e.g. Stark et al. 2014) and in models of photo-dissociation regions (e.g. Visser et al. 2009). In observations, it is clearly seen as an offset between the threshold for the apparition of 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) show, in an 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, 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 ORION-B require new statistical analysis methods using the information in an optimal way, and a derivation of the physical parameters and their associated errors with a rigorous methodology. For instance, the emission of the lowest rotational transitions of the three major isotopologues of carbon monoxide, 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 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 will here apply this technique in the simplest possible model framework, i.e., the emission of lines in Local Thermodynamic Equilibrium (LTE), which can be fully expressed with analytical equations.
The level populations of interstellar molecules result from the balance of collisional (and possibly radiative) excitation and radiative & collisional de-excitation. Therefore the level populations often deviate from LTE conditions because the collisions are not efficient enough to populate all energy levels according to a Boltzmann distribution. With its low dipole moment (0.1 Debye) and high abundance relative to H 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, i.e., 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 theoretically studied. For instance, using non local, non LTE radiative transfer models of a uniform (constant density and temperature) spherical cloud, Bernes (1979) shows 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 is concluded that the LTE model is mostly valid for the ground state transition and deviations from this approximation increase with the quantum number of the upper level (van der Tak et al. 2007).
With a wide range of physical conditions, from bright far-UV illuminated regions to cold and shielded regions through diffuse and translucent gas irradiated by a moderate radiation field, the Orion B molecular cloud is an ideal place to probe to which extent fractionation and selective photodissociation can modify the elemental abundance ratio. It is also a good region to probe 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 computes and analyzes 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 discusses the performance of this estimator.

Description of the data
We will try 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 will compare our results with the dust-traced H 2 column density and dust temperature. This section describes the associated data sets. The 3 mm data were obtained with the IRAM-30m as part of the ORION-B large program. Pety et al. (2017) present in detail the acquisition and reduction of the dataset used in this study. In short, the used data were acquired at the IRAM-30m telescope using the EMIR receiver and Fourier transform spectrometer from August 2013 to November 2014. The frequency range from 84 to 116 GHz was completely sampled at 200 kHz spectral resolution. The J = 1 − 0 lines of the CO isotopologues analyzed here are observed in a single receiver tuning. These lines are thus well inter-calibrated. The absolute flux calibration at 3 mm for the IRAM-30m telescope is estimated to be better than 5%.

1 mm CO lines
The 13 CO(2−1) and C 18 O(2−1) were also observed at the IRAM-30m in 2006 (PI: N. Peretto) using the ABCD generation of receivers and the VESPA auto-correlator. The two lines were observed simultaneously ensuring an excellent inter-calibration.
Data reduction was carried out using the GILDAS 1 /CLASS software. The contribution of the atmosphere was first removed (ON-OFF procedure) and the data were calibrated to the T A scale using the standard chopper-wheel method (Penzias & Burrus 1973). The data were then converted to main-beam temperatures using the standard forward (0.94) and main-beam (0.62) efficiencies for the ABCD receiver around 220 GHz 2 . The resulting absolute flux calibration is estimated to be better than 10%. We subtracted a first order baseline from every spectrum, excluding the velocity range from 5 to 15 km s −1 in the Local Standard of Rest (LST) frame. The spectra were finally gridded into a data cube through a convolution with a Gaussian kernel of full width at half maximum ∼ 1/3 of the IRAM-30m telescope beamwidth at the line rest frequency.

Herschel observations
In order to get independent constraints on the physical conditions in the Orion B cloud, we use 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 .

Field of view
We will jointly analyze the J = 1 − 0 and J = 2 − 1 lines of the CO isotopologues. We thus restrict 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  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 Figure 10.
in offsets (δx, δy) in arcseconds from this projection center. The IRAM-30m angular resolution ranges from 11.5 at 220 GHz to 23.5 at 110 GHz. The position-position-velocity cubes of each line were smoothed to a common angular resolution of 23.5 to avoid resolution effects during the comparison. At a distance of 400 pc (Menten et al. 2007), the sampled linear scales range from ∼ 0.045 pc to ∼ 3 pc.
The spectral and spatial axes were resampled in order to share the same spatial grid and velocity axis for all lines. The spectroscopic observations thus provide position-positionvelocity 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.

Noise
In this paper, the standard deviation of the noise σ b is estimated only on negative values of each spectrum through 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 it assumes that the baselining removed any intensity offset. Table 1 lists the median noise estimated after spectral resampling and angular smoothing.

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 will adapt our formalism to handle such cases, it is not obvious to devise a robust statistical test to deduce the best number of components that must be used. This is particularly true at transitions between regions where the number of required velocity components changes to get a good fit. To address this issue, we used the ROHSA (Marchal et al. 2019) algorithm that makes a Gaussian decomposition based on a multi-resolution process from coarse to fine grid. We only used here 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 The data products associated with this paper are available at http: //www.iram.fr/~pety/ORION-B

Radiative Transfer in Local Thermodynamic Equilibrium
Molecular line emission and absorption in the case of Local Thermodynamic Equilibrium (LTE) are well known (see, e.g., Mangum & Shirley 2015). In this Section, we just summarize the associated notations and equations so that we can 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 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 use in radioastronomy. The data reduction (atmospheric ON-OFF calibration and spectrum baselining to subtract the slowly-varying continuum residual from the receiver and the atmosphere) delivers a noise b that is centered (i.e., with zero-mean), and whose variance can be considered constant over each line profile. We assume that two lines (l ∈ {1, 2}) from the same species are 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 ν red l = ν l 1 − V c , where V is the velocity of the emitting cell of gas in the LSR frame and c is the speed of light. This equation is the radio low-velocity approximation of the Doppler effect. The Doppler effect due to the motion of the observer relative to the LSR is automatically taken into account in the data acquisition process. Therefore, each line of the dataset is analyzed in the LSR frame. In this frame each line is centered around a typical velocity, noted ∆ V . This velocity is related to the redshifted centroid frequency of the line, ν cent l , through a particular case of the previous equation 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 ν red l can be written as where T CMB is the known CMB temperature (T CMB = 2.73 K, Mather et al. 1994), T ex is the unknown excitation temperature along the line of sight, and J is a measure of intensity at a given temperature where B(T, ν) is the spectral distribution of the radiation of a black body at temperature T . The term 1 − exp(−Ψ(ν)) in Eq. (4) represents the emission/absorption by the emitting/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 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 will lead to larger values). The function φ is a Gaussian profile 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 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 its energy (in units of Kelvin), and N 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 The parameter α l is related to the line opacity τ l 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 n up n low = g up g low exp − hν l kT ex .
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 , ν cent l , Ψ l , s, and x will also depend on the sky position.

Cramer-Rao Bound analysis
In this article, we aim at estimating the physical parameters of the LTE model presented in section 3 on the ORION-B data (see section 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) ( θ − θ) 2 , 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 MSE is large, one does not know whether it is due to the choice of the estimator or to a lack of information in the data. In estimation theory, the Fisher matrix allows to quantify the amount of information in the considered problem. It provides a reference precision, named Cramer-Rao Bound (CRB), which does not depend on the choice of a specific estimator of the searched quantity, but only on the physical model and the statistical properties of the noise (see, e.g., Bonaca & Hogg 2018;Espinosa et al. 2018).
Mathematically speaking, the CRB noted B(θ) is simply a lower bound on the MSE of unbiased estimators (see Eq. (15)). Indeed, the MSE is equal to the estimation variance ( θ − θ ) 2 for an unbiased estimator because the estimation MSE is in general equal to the sum of its variance and its bias ( θ − θ) squared, i.e., Therefore, a high CRB value implies that any unbiased estimator θ will necessarily have a high dispersion around the true value θ.
A high CRB can be understood as a lack of information in the underlying model with respect to the considered level of noise.
When it occurs, one solution can be to introduce additional a priori knowledge or to make another measurement with a better signal-to-noise ratio. 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 section 5.1 we check with Monte Carlo simulations that an efficient estimator (i.e., whose variance reaches the CRB) exists.
We here 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. First, what is the maximum noise tolerable to get a given relative precision on these parameters? We here 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, which 12 CO line should be observed to improve the precision achieved when only 12 CO(1−0) observations are available?  (2) can be written over K discrete frequency channels as ∀n ∈ {1, ..., K} x n,l = s n,l + b n,l .
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 line l, can simply be computed as (Stoica & Moses 2005) where [ A] i j stands for the term (i, j) of the matrix A. In our case, the physical model for s will be the LTE radiative transfer introduced in Sect. 3 and the vector of unknown parameters is θ = [T ex , log N, ∆ V , σ V ] T , which we also write θ = [θ 1 , θ 2 , θ 3 , θ 4 ] T 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 θ i (here T ex , log N, Each diagonal term of the inverse of the Fisher matrix B(θ i ) = [I −1 F ] ii is called the Cramer Rao bound of the corresponding parameter. We will note them B(T ex ), B(log N), B(∆ V ) and B(σ V ). These CRBs do not depend on the choice of the estimation algorithm θ and are usually asymptotically reached by the maximum likelihood estimator (Garthwaite et al. 1995). The CRB can thus be considered as a reference precision of the estimation problem.
The calculation of the gradients ∂s n,l ∂θ i i=1,2,3,4 is detailed in appendix A.

Generalization to two lines and two velocity components
We will use the CRB analysis on the case where we observe two different lines (l ∈ {1, 2}) of the same species. We will assume that these lines are well separated in frequency so that their frequency supports are disjoint The Fisher matrix of the set x n,l 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 will also use the CRB framework in the case where each observed line is emitted from two independent velocity components, i.e., from two gas components characterized by different values of the unknown parameters (θ m with m ∈ {1, 2}). Equation 16 that encodes the spectrum for line l can then be written as where S n,l = s n,l (θ 1 ) + s n,l (θ 2 ).
This means that the composite line profile is considered as the simple sum of two velocity components that do not radiatively interact. This assumption is only correct if the two velocity components are sufficiently separated in velocity. This case with two components is the most complex model we will study in this paper. In this case, the number of unknown parameters is 8 (instead of 4) and thus the size of the Fisher Matrix is 8 × 8 (instead of 4 × 4).

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 thus empirically analyze their evolution as a function of the physical properties of the analyzed medium in a particular case taken from the ORION-B project. To generate figures 2, 3, and 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 13 CO(1−0) 13 CO(2−1) 13 CO(1−0) and 13 CO(2−1)  (d) The values of the Cramer-Rao Bounds of the unknown parameters (i.e., T ex , log N, ∆ V and σ V ) are then computed as a function of the values of T ex and N. The Fisher matrices are computed following Eq. (14), and then numerically inverted to obtain The excitation temperature T ex is sampled logarithmically between 3 K 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 to figures 2, 3, and 4, where N and T ex varies horizontally, and vertically, respectively. In these figures, the variations of B 1/2 (θ i ) are shown instead of the variations of B(θ i ) 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 signal-to-noise ratio, but also the amount of information measured by the Fisher matrix because of the non linearity in the radiative transfer equation.

Precision of the estimation of the centroid velocity ∆ V
and the associated velocity dispersion σ V Figure 2(a-b) shows variations of B 1/2 (∆ V ) and B 1/2 (σ V ). 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 will have a small dispersion around the actual values. This can be written ∆ V = 1.10 ± 0.01 km s −1 and σ V = 0.61 ± 0.01 km s −1 . Figure 2(c) shows a function of KR, where R is the signalto-noise ratio defined by This expression is used in signal processing to quantify the signal-to-noise ratio on the "energy" of the signal. In our case, we empirically find that, as a rule of thumb 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 signal-to-noise ratio and on σ 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). 13 CO(2−1) 13 CO(1−0) and 13 CO(2−1) (a) B 1/2 (log N) (j) log τ 1 (k) log τ 2 (l) log τ 1 log τ 2 we first mention 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 monotonous 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 Fig. 3 (a-b)) will give a reference precision on T ex of 0.1 K ≤ B 1/2 (T ex ) < 10 K 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 [1 − exp(−Ψ(ν cent l ))] in Eq (4) approaches unity. The CRB almost linearly depends on log T ex above 6 K. The analysis of two lines (see Figure 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, i.e., 0.1 K ≤ B 1/2 (T ex ) < 10 K for typically N > [10 14 − 10 16.5 ] cm −2 .
Here again the precision almost linearly depends on log T ex above 6 K.
The second row of Figure 3 shows functions of the opacities. Trying for several values of σ V , we empirically obtain 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 Hence, according to Eq. (21) and Eq. (22), when opacities are close to one, the gain in precision (in standard deviation) is around 20 when one observes 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.

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 log N estimations is given by where B(T ex , log N) = [I −1 F ] 12 is the non-diagonal element of the inverse Fisher matrix, see Eq. (14). We also introduce the correlation coefficient between log N estimations and σ V estimations where B(log N,  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 log N) or (log N and σ V ), only when values of |γ| = 1. In this case, the variance of these estimations becomes infinite. Figure 5 shows a simulated example where log(N) and σ V can be accurately estimated even though they are highly (but not completely) anti-correlated. Starting from a modeled spectrum with log(N/ cm −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 log N and T ex and even more ambiguities between estimations of log N and σ V (in particular for large values of N). When a single line is observed (Fig. 4 d e, g, h), |γ(T ex , log N)| and |γ(log N, σ V )| are mostly larger than 0.9 for small N (in Fig. 4 d-e yellow pixels correspond to γ > 0.99 and in g-h to γ > 0.9). For high N, the ambiguity with T ex decreases, but not the one with σ V (in Fig. 4 d-e light blue values are −0.5 < γ < 0 while in Fig. 4 g-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 γ. Figure 4 f and 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 Fig. 4 a-b), B 1/2 (log N) < 0.1 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. 4 c), the situation greatly improves: B 1/2 (log N) < 0.1 for N > [10 15 − 10 17 ] cm −2 (depending on T ex ) and T ex ≥ 6 K. Figures 4 a, b, c also show local minima of the B 1/2 (log N) 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 B 1/2 (log N) shows that the smallest values of B 1/2 (log N) (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.

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 where ρ is a fixed threshold. In other words, instead of analyzing the precision for a given amount of noise, one can also ana-lyze the tolerable level of noise to ensure an intended precision (herein described by ρ). Such an analysis will be useful to design an observation program and optimize the telescope time needed to reach the scientific goal. Up to this point of the paper, we 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-30m time estimator for the EMIR receivers 5 indicates that we can achieve a sensitivity of 100 mK in 30 to 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 1 to 3 hours depending on the frequency and the observing mode (frequency or position switching). The colored contours thus correspond to the "fast/slow" acquisition mode at the IRAM-30m. 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, one can also analyze the peak-signal-to-noise ratio defined for one transition l by  Figure 6, except that the contours show the variations of the peak-signal-to-noise ratio, P ρ , required to reach a given relative accuracy (ρ%). and for two transitions by P = max l=1,2 P l . Figure 7 shows the variations of minimum peak-signal-to-noise ratio P ρ for similar conditions as in Fig. 6. Analyzing only the 13 CO(1−0) line requires at least a signal-to-noise ratio of 100 to get a relative precision of 20%. Adding the 13 CO(2−1) line in the analysis reduces the minimum signal-to-noise ratio by a factor up to 10 to reach the same relative precision. The required signal-to-noise ratio increases at high column densities because the lines become optically thick, and at a combination of low column density and high excitation.
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). We here ask two questions. First, what is the best J line to observe 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 a CRB point of view) because it allows to reach 20% accuracy over this broad range of parameters for a noise level of 100 mK. However this line is not easy to access from ground based tele-scopes because of the limited atmospheric transmission at the line frequency of 690 GHz.
Second, what is the best J line to be observed 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 keep a good 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 higher-J transitions. Nevertheless, the usefulness of mildly excited lines remains valid. Non-LTE approaches will be developed in the future that will provide a quantitative assessment of the diagnostic power of these lines.

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

Proposed estimator
The maximum likelihood estimator (MLE) is a good candidate because under mild conditions, it reaches the CRB asymptotically, i.e., 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 (i.e., 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.

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 com-ponent will 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 signal-to-noise ratio, 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 signal-to-noise ratio 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 explained above. 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. 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 log-likelihood is done independently on each species, ensuring that the estimations of ∆ V and σ V may take a different value for each species.
Second, a single line of 12 CO is available and it is quite optically thick. In our analysis of the CRB, we observed in Figure 4 (d-e and g-h) that the estimation of log N 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 12 CO(1−0) and 12 CO(2−1) 12 CO(1−0) and 12 CO(3−2) 12 CO(1−0) and 12 CO(4−3) 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 log N 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. 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.

Detailed analysis of two lines of sight
The spectra on the right column are well fitted 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.

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 will 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 E s l = n s 2 n,l ( θ) or E s l = n s 2 n,l ( θ 1 ) + s 2 n,l ( θ 2 ), depending of the number of estimated components. -The "energy" of the fit residual is where r n,l = x n,l − s n,l ( θ), or r n,l = x n,l − s n,l ( θ 1 ) − s n,l ( θ 2 ).
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 (E r l /E x l ) or the difference of "energy" between the observed and estimated signals ({E x l − E s l }/E x l ). The former formula tells us the fraction of the observed "energy" that has One velocity component (δx = 803 , δy = 473 ) Two velocity components (δx = 578 , δy = −121 ) 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).
not been fitted. The latter formula tells us whether the observed spectra has been under-fitted (positive value) or over-fitted (negative value). Measuring the residual "energies" is similar to computing a χ 2 in least-square fitting. Another way to quantify the quality of the fit is to compare the standard deviation of the residuals σ r = 1 K−1 n r 2 n,l 1/2 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 Sect. C.1, an estimation will 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 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 well fitted. We also checked that the residuals are larger (i.e., σ r > σ b ) when the energy E x is large (not shown in the figure). As the signal-to-noise ratio also increases with E x , this issue implies some misspecification of 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 will 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 self-absorption signatures that may happen at large opacity. Figures 12 a and 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 look-up table. On the left image, the yellow pixels correspond to bright molecular gas while the blue pixels corresponds to faint signal or noise associated with the IC 434 Hii region. The "energy" of the residual still exhibits spatially coherent structures, but at a much smaller level than the "energy" of the measured spectra. Our estimator under-fits the observed spectra as shown in Fig. 12d. This may be related to the fact that our model does not fit asymmetric profiles. However, the fit quality is good as illustrated by the image of the energy ratio, which shows that the residual "energy" amounts to less than 1% of the signal (the dark blue color corresponds to 0.01% in image 12c) except in regions where the signal-to-noise ratio 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 it already encodes. The fact that the residuals are smaller than 1% compared to the observed signal is sufficient to make the analysis of the excitation temperature, the column density, and the velocity dispersion pertinent for CO isotopologues. The main source of systematic errors in the column density determination results from the deviations from the local thermodynamic equilibrium, leading to a more complex partition function than the simple formula in Eq. (9). The effect is expected to be stronger in warm regions (T ex ≥ 50 − 100 K) where many rotational levels are populated and contribute to the partition function. These warm regions occupy a small fraction of the total volume and therefore a bias would not affect the general conclusions. Non-LTE approaches will be developed to assess more quantitatively the magnitude of the effect and to provide recommendations on the best method depending on the molecular lines and the range of physical conditions that are studied.

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  Relative precision on the centroid velocity, velocity dispersion, excitation temperature, column density, and line peak signal-to-noise ratio. The black contours on the peak signal-to-noise-ratio image delimit the regions where P ≥ 3. 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 latter case, the images only show the estimation that is the closest (in terms of centroid velocity ∆ V ) to its neighboring pixels.
dispersions and opacities. These values are computed over the field of view that is observed for all the lines.
The peak-signal-to-noise ratio of the 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 peak-signal-to-noise 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 T ex ( 13 CO) 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. Fig. 15. Scatter plots between the CO isotopologue excitation temperatures. The color scale encodes the dust temperature T dust . The ellipses represent 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. 12 CO signal-to-noise ratio is much larger than the C 18 O or the 13 CO one, the 12 CO estimations are more uncertain (see map of B 1/2 (log N) for 12 CO in Figure 13) because a single transition is available and the 12 CO opacities are large (ranging from 5 to about 1 000). However, the higher signal-to-noise ratio 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 species provide similar kinematics information. The velocity field is spatially regular with well-resolved gradients, for instance in the Horsehead nebula at the bottom right of the map. The velocity Table 3. Statistics of ratios of the estimated parameters over all the pixels for the three CO isotopologues. We stress that σ V ( 12 CO) is fixed to σ V ( 13 CO). 13 CO/C 18 O 12 CO/ 13 CO T ex /T ex 1.3 ± 0.33 1.7 ± 0.39 log(N/N ) 1.2 ± 0.17 1.4 ± 0.33 σ V /σ V 1.1 ± 0.25 Fixed to 1 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 Hily-Blant et al. 2005), and somewhat narrower lines for C 18 O than for 13 CO.

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 questioned. Recently, in their clustering analysis of a one square degree map in the Orion B cloud, Bron et al. (2018) have shown 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 present. 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 T ex ( 12 CO)/T ex ( 13 CO) ∼ 1.7 and T ex ( 13 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 diffuse/translucent cloud along the line of sight to HD62542.
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 photo-dissociation regions and UV-dominated regions where the gas kinetic temperature is larger than the dust temperature. This is different from the usual approximation that T dust is a good approximation of the gas kinetic temperature in cold and shielded regions.
The difference in excitation temperatures between CO isotopologues can be explained by radiative trapping in the 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 column-to-intensity ratio, defined as the ratio of the column density of the species to the line integrated emission, as a function of the excitation temperature. Figure 16 shows the variation of the column-tointensity 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 N13 CO ∼ 10 14 − 10 15 cm −2 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 column-to-intensity ratio presents a shallow minimum which depends on the transition, rises rapidly at temperatures lower than the minimum and more slowly for excitation temperatures above the minimum. When the line opacity becomes significant, the shape of the column-to-intensity ratio curve changes and the 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 column-to-intensity ratio is inappropriate. For moderate column densities (N12 CO < 10 16 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 column-to-intensity curve. A CO excitation temperature near the minimum of the column-to-intensity curve will minimize the error because a small difference in T ex in this region will not change the column-to-intensity. For high column densities (N ∼ 10 17 cm −2 ), the bias is significant at low excitation temperatures (T ex < 20 K) because the column-to-intensity ratio (purple curve in Fig. 16) is rising fast at low excitation temperatures. 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 N( 13 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). But these values are about a factor of two lower than those predicted when one assumes no isotopic fractionation, and uses the non depleted gas phase carbon elemental abundances applicable to the Orion region, 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 N( 13 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 freeze-out 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 photo-illuminated edges where molecules are photodissociated, but also in high column density and well shielded regions. In these latter regions, the dust temperature gets below the CO condensation temperature (T dust < 17 K), and the CO molecules can rapidly freeze onto dust grains, lowering the gas phase abundances. This depletion effect is seen for C 18 O and 13 CO supporting the explanation by a global freeze-out 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).

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 N( 13 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) and/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 N( 12 CO)/N( 13 CO) 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 will further increase 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 assumptions is also expected to increase the dispersion and affect the overall determination of the X CO = N(H 2 )/W(CO) conversion factor.

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 ap-parent 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 however fitted 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.

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 Rao bound (CRB) computation. We applied this analysis to the region of the Orion B molecular cloud that contains the Horsehead pillar and 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 N( 13 CO)/N(C 18 O), 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 will lead to large errors on the estimation of the column density. -This analysis also allows us to quantify the benefit of a simultaneous analysis of two rotational lines of the same species compared to the analysis of a single line. In particular, it alleviates the degeneracy described above. This is a rigorous demonstration of intuitive results. It is an argument in favor of the installation of dual-band receiver systems for telescopes like the IRAM-30m 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 pre-processing) and thus applying only a gradient on σ V , T ex and N could improve the robustness of our estimations. This would be useful, in particular, for low signal-to-noise ratio 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 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.
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 1. Spatial variations of the observed spectrum "energy" (first column), of the residual "energy" (second column), of their ratio (third column), and of the ratio of "energy" that has not been modeled. The unit of the color look-up table is Kelvin 2 or % depending on the column. White contours show the regions where two components have been detected.