ATLASGAL-selected massive clumps in the inner Galaxy. X. Observations of atomic carbon at 492 GHz

(Abridged) In this paper, we investigate the physical conditions of [CI]-traced gas in high-mass star-forming regions by analyzing APEX [CI] 492 GHz single-pointing observations of the ATLASGAL Top100 sources along with other multi-wavelength data. Our 98 sources are clearly detected in [CI] 492 GHz emission, and the observed integrated intensities and line widths tend to increase toward evolved stages of star formation. In addition to these"main"components that are associated with the Top100 sample, 41 emission and two absorption features are identified by their velocities toward 28 and two lines of sight respectively as"secondary"components. The secondary components have systematically smaller integrated intensities and line widths than the main components. We found that [CI] 492 GHz and 13CO(2-1) are well correlated with the 13CO(2-1)-to-[CI] 492 GHz integrated intensity ratio varying from 0.2 to 5.3. In addition, we derived the H2-to-[CI] conversion factor, X(CI), by dividing 870 micron-based H2 column densities by the observed [CI] 492 GHz integrated intensities and found that X(CI) ranges from 2.3e20 to 1.3e22 with a median of 1.7e21. In contrast to the strong correlation with 13CO(2-1), [CI] 492 GHz has a scattered relation with the 870 micron-traced molecular gas. Finally, we performed LTE and non-LTE analyses of the [CI] 492 GHz and 809 GHz data for a subset of the Top100 sample and inferred that [CI] emission likely originates from warm (kinetic temperature>60 K), optically thin (opacity<0.5), and highly pressurized (thermal pressure ~ e5 to e8 K/cm3) regions.


Introduction
High-mass stars play a key role in the evolution of galaxies (e.g., Kennicutt & Evans 2012). For instance, massive stars inject a substantial amount of radiative and mechanical energy into the surrounding gas via ultraviolet (UV) radiation fields, stellar outflows/winds, and supernova explosions. This stellar feedback has been argued to be critical for regulating star formation in galaxies (e.g., Cox 1981;Silk 1997;Ostriker et al. 2010;Hopkins et al. 2014). In addition, high-mass stars enrich the interstellar medium with heavy elements throughout their lives, driving the chemical evolution of galaxies (e.g., Matteucci 2021).
Despite their importance, how high-mass stars form and interact with the surrounding gas is not yet fully understood (e.g., Motte et al. 2018). One of the challenges has been that massive star-forming regions tend to be widely separated and located far awary from us due to the rarity of high-mass stars. Furthermore, high-mass stars reach the main sequence while deeply embedded in their parental molecular clouds, preventing us from prob-ing the early stage of star formation with traditional optical and near-infrared (NIR) observations. To overcome these difficulties and provide a comprehensive census of high-mass star-forming regions, a number of Galactic plane surveys at far-infrared (FIR) and submillimeter (submm) wavelengths have been performed in recent years, including the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009). ATLASGAL is a 870 µm survey of the inner Galaxy (covering |l| < 60 • with |b| < 1.5 • and 280 • < l < 300 • with −2 • < b < 1 • ) conducted with the Atacama Pathfinder EXperiment (APEX) 12 m submillimeter telescope in Chile. This survey identified ∼10 4 massive clumps (e.g., Contreras et al. 2013;Csengeri et al. 2014;Urquhart et al. 2014), and the brightest ∼100 of them (Top100 hereafter) have been a subject of various follow-up studies, e.g., Giannetti et al. (2014) (investigating CO depletion and isotopic ratios), Csengeri et al. (2016) (probing SiO-traced shocked gas), Kim et al. (2017) (searching for H ii regions), and König et al. (2017) (characterizing physical properties based on dust continuum emission). More specifically, König et al. (2017) constructed dust spectral energy distributions (SEDs) from 8 µm to Article number, page 1 of 27 A&A proofs: manuscript no. draft_040622 870 µm and measured dust temperatures and fluxes. Along with distance information, the authors then estimated clump masses, luminosities, and column densities on ∼19 ′′ scales. These dustbased analyses showed that the Top100 sample could be divided into four groups, i.e., 70 µm weak sources (70w), mid-infrared (MIR) weak sources (IRw), MIR-bright sources (IRb), and MIRbright sources associated with radio continuum emission (H ii regions), and the dust temperature and bolometric luminosity increase along this sequence, implying that the Top100 sample represents different evolutionary stages of high-mass star formation.
We have been conducting a multi-telescope campaign for observing carbon species in the Top100 clumps to address the following questions: (1) What are the physical conditions of the interstellar medium leading to high-mass star formation? (2) How do the energetics of gas surrounding massive young stellar objects and protostellar objects vary as the sources evolve? Carbon is a critical ingredient in the evolution of the interstellar medium (e.g., Henning & Salama 1998). In its different forms, i.e., C + , C 0 , and CO, it is one of the primary coolants and provides a powerful diagnostic tool for examining the physical conditions and energetics of gas. For example, [C ii], [C i], and low-J CO transitions (upper J 4) can be used to probe the properties of photodissociation regions, while mid-and high-J CO lines can trace shock-heated gas (e.g., Pon et al. 2014;Lee et al. 2019). Among these carbon phases, atomic carbon has received relatively little attention, mainly because [C i] finestructure transitions are typically fainter than low-J CO emission and are heavily affected by atmospheric absorption. Theoretically, atomic carbon is expected to be abundant at the surfaces of molecular clouds (e.g., Langer 1976), but subsequent observations have shown that [C i] emission is widespread throughout the clouds, invoking an interest in [C i] emission as a tracer of total gas mass (e.g., Frerking et al. 1989;Schilke et al. 1995;Kramer et al. 2008;Oka et al. 2001;Shimajiri et al. 2013). This is the first of a series of papers presenting results from our [C ii] and [C i] campaign for the Top100 sources. As a pilot study, we here focus on examining the physical conditions of [C i]-traced gas based on APEX [C i] 492 GHz single-pointing observations, and the two science questions described above will be addressed with more data in forthcoming papers. The organization of this paper is as follows. In Sect. 2, we present APEX [C i] and 13 CO(2-1) observations of the Top100 sources and describe the determination of line parameters (including central velocities and line widths) via Gaussian fitting. The observed properties of [C i] 492 GHz are probed across different evolutionary stages of high-mass star formation (Sect. 3) and are compared to molecular gas tracers (e.g., 13 CO(2-1) and dust-based H 2 ; Sect. 4). In addition, the [C i] 492 GHz and 809 GHz data for a small sub-set of the Top100 sources are analyzed to derive the physical conditions of [C i]-traced gas (Sect. 5). Finally, we discuss our results and summarize main findings in Sect. 6.

[C i] 492 GHz
Single-pointing observations of the [C i] 3 P 1 -3 P 0 transition at 492.16065 GHz were made toward 98 of the Top100 sources using the APEX FLASH + (First Light APEX Submillimeter Heterodyne instrument; Klein et al. 2014) receiver in 2016 (project IDs: M0039-97 andM0013-98). Two different off positions were observed in the total power mode to avoid contamination, while additional observations were performed in the wobbler mode to accurately measure continuum levels. The characteristic of the observed sources, including locations, distances, and classifications, are summarized in Table A.1. The obtained spectra were processed with the GILDAS CLASS 1 software and were converted into units of main-beam brightness temperature (T MB ) using a forward efficiency (F eff ) of 0.95 and a main-beam efficiency (B eff ) of 0.46 based on Mars observations. The final spectra on 13 ′′ scales were smoothed to a velocity resolution of 0.46 km s −1 and have a root-mean-square (rms) noise level of 0.1-0.6 K (median ∼ 0.3 K).

[C i] 809 GHz
Nine of the Top100 sources were mapped in the [C i] 3 P 2 -3 P 1 transition at 809.34197 GHz using the APEX CHAMP + (Carbon Heterodyne Array of the MPIfR; Güsten et al. 2008) instrument in 2013. The mapping was made in the On-The-Fly (OTF) mode with a coverage of 1.6 ′ × 1.6 ′ , and the observed sources are marked in Table A.1. The GILDAS CLASS software was employed for data reduction, and the processed data cubes were translated into the T MB units by adopting F eff = 0.95 and B eff = 0.38 based on Mars observations. The final data on 8 ′′ scales were smoothed to a velocity resolution of 0.70 km s −1 and have a median rms noise level of ∼1 K.

13 CO(2-1)
The 98 sources in our [C i] 492 GHz survey were also observed in the 13 CO(2-1) transition at 220.39868 GHz (single-pointing) using the APEX PI230 receiver in 2016 (project ID: M0024-97). Unlike the [C i] 492 GHz observations, only one off position was observed for each source in the total power mode. The obtained 13 CO(2-1) data were reduced using the GILDAS CLASS software and were converted into the T MB units by assuming F eff = 0.95 and B eff = 0.66 (Tang et al. 2018). The final science-ready spectra have a resolution of 30 ′′ and a rms noise level of 0.04-0.14 K (median ∼ 0.06 K) per 0.66 km s −1 channel.

Deriving line parameters
While the observed [C i] 809 GHz spectra can be approximated as single Gaussians, the [C i] 492 GHz and 13 CO(2-1) spectra often show complex shapes with multiple components (Appendix B). For our analyses, we identified individual components by fitting Gaussian functions to the final spectra and estimated line parameters such as central velocities and full width at half maximum (FWHMs). When two components are close to each other, we examined the difference between their central velocities and considered them as one component if the difference is smaller than the sum of their FWHMs. In addition, we calculated the integrated intensity (I([C i]) or I( 13 CO)) of each identified component by integrating the emission over a velocity range over which the component is clearly detected. For [C i] 809 GHz, we smoothed the cubes and used the spectra that were extracted at a resolution of 13 ′′ to compare with the [C i] 492 GHz data (Sect. 5). The final [C i] and 13 CO(2-1) spectra with the fitted Gaussians are shown in Fig. B.1, and the derived line parameters and integrated intensities are presented in Table A.2. GHz spectra of the four evolutionary groups (70w, IRw, IRb, and H ii region groups in blue, green, yellow, and red). To derive these spectra, sources with contaminated reference positions were excluded.

Observational results
In this section, we mainly discuss the observed properties of [C i] 492 GHz (integrated intensities and FWHMs) and how they vary in a diversity of environments.

[C i] 492 GHz detection rate
We regard components whose peak-to-rms ratios are equal to or higher than five as detections. With this threshold, we found that all 98 sources are detected in [C i] 492 GHz emission. The exact breakdown of the sources is as follows: 14, 31, 31, and 22 for the 70w, IRw, IRb, and H ii region groups respectively. These different evolutionary groups have comparable median distances (∼4 kpc) and [C i] 492 GHz sensitivities (∼0.3 K), making comparisons of the groups straightforward. In addition to the main components that are directly associated with the Top100 sources, secondary components are detected in both emission and absorption (Fig. B.1): 41 emission components toward 28 lines of sight and two absorption components toward two lines of sight. Throughout this paper, we refer to these components as "diffuse clouds", since they likely trace colder and less dense gas than the main 98 components.

Average [C i] 492 GHz spectra
To examine how [C i] 492 GHz changes across different environments, we first derived average spectra of the four evolutionary groups by shifting the observed spectra to a central velocity of 0 km s −1 based on the source velocities and stacking them (Fig. 1). For this, 94 sources with clean reference positions were considered (12, 30, 30, and 22 for the 70w, IRw, IRb, and H ii region groups respectively). Fig. 1 shows that the characteristics of [C i] 492 GHz emission systematically change across the four evolutionary stages in terms of brightness and line width. Specifically, the peak temperature increases from 5.8 K to 7.8 K for the IRb and H ii region groups, while the 70w spectrum has a slightly higher peak than the IRw spectrum (4.7 K versus 3.8 K). On the other hand, the FWHM increases toward more evolved stages: 4.3, 5.6, 5.9, and 8.4 km s −1 for the 70w, IRw, IRb, and H ii region groups. The peak temperature and FWHM values here were estimated from Gaussian fitting of the average spectra.

Integrated intensities and line widths of the main components
Next we probed how statistically different [C i] 492 GHz properties are between the evolutionary stages by constructing cumulative distribution functions (CDFs) of the integrated intensity and FWHM for the main components (Fig. 2). Like in Sect. 3.2, we used the 94 sources with clean reference positions only. The top panel of Fig. 2 suggests that the 70w, IRw, and IRb groups have comparable distributions of the integrated intensity, while the H ii region group has systematically higher values (median I(C i) ∼ 23, 21, 26, and 67 K km s −1 for the 70w, IRw, IRb, and H ii region group respectively). This result is supported by two-sided Kolmogorov-Smirnov (K-S) tests 2 at a significance level of 5%, where only the H ii region group is found to be statistically distinct from the rest. Similarly, we found that the 70w and H ii region groups have the smallest and largest FWHMs, while the IRw and IRb groups have intermediate FWHMs (Fig.  2 bottom): median FWHM ∼ 3.7, 4.6, 4.4, and 7.0 km s −1 for the 70w, IRw, IRb, and H ii region group respectively. This conclusion again agrees with K-S tests where the 70w group is found to be statistically different from the IRw and H ii region groups and the H ii region group stands out against the rest.

Comparison between the main and diffuse cloud components
As noted in Sect. 3.1, 41 components are detected along 28 lines of sight as secondary emission components. In comparison to the main components that are associated with the Top100 sources, these secondaries have systematically smaller values of the [C i] 492 GHz integrated intensity and line width (Fig. 3). For example, the median I(C i) is 29 K km s −1 and 4 K km s −1 for the main and diffuse cloud components respectively, while the median FWHM is 5 km s −1 and 2 km s −1 . For these estimates and the histograms in Fig. 3, spectra with clean reference positions were used (94 and 38 for the main and diffuse cloud components respectively). In summary, our APEX survey suggests that [C i] 492 GHz is ubiquitous and traces not only massive dense clumps, but also non-star-forming relatively diffuse gas. The observed integrated intensities and FWHMs systematically vary across different environments, which is likely driven by local conditions of the interstellar medium.

Relation between [C i] 492 GHz and molecular gas tracers
From which part of a molecular cloud [C i] 492 GHz emission arises has long been a subject of debate. Considering the small difference between the ionization threshold of C 0 (11.26 eV) and the dissociation threshold of CO (11.09 eV), C 0 is expected to form a narrow layer between C + and CO, limiting its usefulness as a tracer of total gas mass. However, a significant body of research has revealed that this expectation could be wrong and C 0 could exist deep inside molecular clouds (e.g., Frerking et al. 1989;Schilke et al. 1995;Kramer et al. 2008). For example, the proportionality between the [C i] 492 GHz and low-J 13 CO integrated intensities has often been interpreted as evidence of the widespread [C i] 492 GHz emission throughout molecular clouds (e.g., Keene et al. 1997;Oka et al. 2001;Shimajiri et al. 2013;Izumi et al. 2021 On the other hand, 45 secondary components are newly detected in 13 CO(2-1) only, which is in contrast to a handful number of components seen in [C i] 492 GHz only (5). This striking difference partially results from the higher sensitivity of the 13 CO(2-1) data (Sect. 2).

[C i] 492 GHz as a tracer of "CO-dark" molecular gas
The [C i] 492 GHz components without corresponding 13 CO(2-1) detections can be candidates for "CO-dark" molecular gas and require an in-depth look. In Fig. 4, we compare [C i] 492 GHz, 13 CO(2-1), H i, and 12 CO(1-0) for the five sources toward which [C i] components without association with 13 CO(2-1) emission are identified (indicated as gray dashed lines). For G08.71−0.41, Article number, page 4 of 27   4 shows that G08.71−0.41 and G14.49−0.14 likely suffer from bad reference positions in 13 CO(2-1) at the velocities of our interest (−9 km s −1 and 55 km s −1 respectively), which means that 13 CO(2-1) could have been detected. On the other hand, 13 CO(2-1) is simply too faint at the velocities of our interest for G336.96−0.23, G337.92−0.48, and G351.57+0.77 (−101 km s −1 , −19 km s −1 , and −16 km s −1 respectively), implying the I( 13 CO)-to-I([C i]) ratios lower than 0.7, 0.9, and 0.5 respectively. For these estimates, 13 CO(2-1) 5σ values for the velocity ranges where [C i] emission is integrated are considered. In general, Fig. 4 shows that H i is abundant toward the five sources. In particular, for G337.92−0.48, the [C i] component without associaiton with 13 CO(2-1) emission corresponds to the strongest absorption feature, implying a substantial contribution from the cold H i . On the contrary, 12 CO(1-0) is not detected at the velocities of our interest for G336.96−0.23 and G337.92−0.48 (rms of ∼1-2 K). These results suggest that the [C i] 492 GHz secondary components without corresponding 13 CO(2-1) detections have a large amount of atomic gas and could trace "CO-dark" molecular gas. This conclusion, however, is drawn from a couple lines of sight and needs to be confirmed by extensive observations of Article number, page 5 of 27 is the velocity resolution of the [C i] spectra. A similar calculation was done for the Type 2, but using the number of 13 CO channels over which [C i] is detected, 13 CO 5σ, and velocity resolution of the 13 CO spectra instead. These derived upper limits are indicated as the arrows. Left: log scale. The measured minimum, median, and maxium I( 13 CO)-to-I([C i]) ratios (0.2, 1.8, and 5.3) are indicated as the dashed lines in different shades of gray color. Right: linear scale. A linear one-to-one relation is shown as the black dashed line. [C i] 492 GHz, 13 CO(2-1), and 12 CO(1-0) with matching resolutions and sensitivities.

Correlation between the [C i] and 13 CO integrated intensities
In Fig. 5, we show a comparison between the [C i] 492 GHz and 13 CO(2-1) integrated intensities for all three types. To provide different perspectives, both linear and log versions are presented. For this comparison, we excluded 41 detections that suffer from bad reference positions in [C i] 492 GHz or 13 CO(2-1) and used arrows to mark upper limits for the Types 1 and 2.
In general, we found that the [C i] 492 GHz and 13 CO(2-1) integrated intensities are well correlated. For instance, the Type 0 detections have a strong correlation between [C i] 492 GHz and 13 CO(2-1) with a Spearman rank correlation coefficient of 0.95. This suggests that [C i] 492 GHz and 13 CO(2-1) arise from similar regions of the interstellar medium, which could be supported by comparable line widths of the two transitions ( Fig. 6): the median of the absolute difference between the [C i] 492 GHz and 13 CO(2-1) line widths is 0.8 km s −1 . Considering a factor of two larger resolution of the 13 CO(2-1) data (13 ′′ versus 30 ′′ ), this difference in FWHM is surprisingly small. Finally, both [C i] 492 GHz and 13 CO(2-1) are likely not too optically thick (Sect. 5.2; Urquhart et al. 2021), contributing to the good correlation between the two transitions.

Ratio of 13 CO to [C i]
While [C i] 492 GHz and 13 CO(2-1) are strongly correlated, the ratio of 13 CO(2-1) to [C i] 492 GHz does vary, e.g., the I( 13 CO)to-I([C i]) ratio ranges from 0.2 to 5.3 with a median of 1.8 (Fig.  5 left). These variations can be seen in Fig. 5 (right) as well, where the data points deviate from the linear one-to-one relation more significantly as the [C i] and 13 CO integrated intensities increase beyond ∼50 K km s −1 . Considering that I([C i]) increases toward more evolved groups (Sect. 3.3), the systematic deviation from the linear one-to-one relation implies that the ratio of 13 CO(2-1) to [C i] 492 GHz depends on the classification of the sources.  To examine how the ratio of 13 CO(2-1) to [C i] 492 GHz depends on the classification of the sources, we constructed CDFs of I( 13 CO)/I([C i]) for the 70w, IRw, IRb, H ii region, and diffuse cloud groups and found that the diffuse cloud and 70w groups have systematically lower ratios (Fig. 7). For example, the fraction of the sources with I( 13 CO)/I([C i]) < 2 are 87%, 80%, 48%, 57%, and 41% for the diffuse cloud, 70w, IRw, IRb, and H ii region groups respectively. This conclusion is also supported by two-sided K-S tests, where the diffuse cloud and 70w groups are found to be drawn from the same distribution, while the null hypothesis is rejected at a significance level of 5% when the diffuse cloud components are compared against the IRw/IRb/H ii region groups.
The observed variations in the 13 CO(2-1)-to-[C i] 492 GHz integrated intensity ratio most likely result from local conditions of gas and implies that the Top100 sources are in different stages of carbon phase transition (C + -, C 0 -, or CO-dominated). The exact stages of carbon phase transition, as well as the underlying physical and energetic conditions, will be examined in forthcoming papers where observations of important cooling lines ([C ii], [C i], [O i], and CO transitions) for the Top100 sources are analyzed using state-of-the-art models of photodissociation regions and shocks (e.g., Le Petit et al. 2006;Flower & Pineau des Forêts 2015).

[C i] 492 GHz and 870 µm-based
In addition to 13 CO(2-1), we examined a relation between [C i] 492 GHz and the total H 2 column density by calculating the ratio of N(H 2 ) to I([C i]) (X(C i) factor). For this calculation, we extracted the 870 µm-based N(H 2 ) estimates from König et al. (2017) for the 94 main components whose [C i] 492 GHz spectra have clean reference positions and scaled them by (13/19) 2 to compensate the difference in resolution (13 ′′ and 19 ′′ for [C i] 492 GHz and H 2 respectively; this scaling assumes that H 2 uniformly fills the 19 ′′ beam). The derived X(C i) ranges from 2.3 × 10 20 to 1.3 × 10 22 with a median of 1.7 × 10 21 (in units of cm −2 (K km s −1 ) −1 ; Fig. 8 left), and this median value is consistent with recent observations of the massive star-forming region RCW 38 (1.4 × 10 21 ; Izumi et al. 2021), as well as predictions from numerical simulations (∼10 21 ; Offner et al. 2014;Glover et al. 2015). While showing a Gaussian distribution with the well-defined peak of 1.7 × 10 21 , X(C i) has large variations (a factor of 55). To examine if there is any trend in these variations, we constructed CDFs for the four evolutionary groups (Fig. 8 right). The constructed CDFs show that the H ii region group has more sources with high X(C i) values. For example, sources with X(C i) > 1.7 × 10 21 constitute 77% of the H ii region group, while accounting for 25%, 57%, and 33% of the 70w, IRw, and IRb groups. When examined at a significance level of 5% for two-sided K-S tests, the H ii region group was found to be statistically different from the 70w and IRb groups.
Sect. 4.2.1 showed that N(H 2 ) for a given I([C i]) is not constant and varies across the evolutionary groups. To better understand the origin of these variations, we examined how X(C i) changes with N(H 2 ) (Fig. 9) and found that X(C i) has a weak dependence on N(H 2 ) (Spearman's rank correlation coefficient of 0.6; we did not attempt to perform least-squares fitting, since proper fitting is impossible without uncertainties in the H 2 column density). The same dependence of X(C i) on the H 2 column density was also found by Offner et al. (2014) and Izumi et al. (2021). In particu-Article number, page 7 of 27 A&A proofs: manuscript no. draft_040622 lar, Offner et al. (2014) performed hydrodynamic simulations of Galactic molecular clouds with two different UV radiation fields (G UV = 1 and 10 Draine fields) and showed that X(C i) increases with N(H 2 ) 10 22 cm −2 (gray lines in Fig. 9). This increase of X(C i) with N(H 2 ) was interpreted as a result of [C i] 492 GHz becoming optically thick and less sensitive to high H 2 column densities. In our case, however, [C i] 492 GHz is likely not too optically thick (Sect. 5). Our result could rather indicate the flattening of C 0 abundance at high H 2 column densities, which can be verified by a comprehensive excitation study based on extensive [C i] 492 GHz and 809 GHz observations.

Relation between I([C i]) and N(H 2 )
The observed variations in X(C i) reflect a scattered relation between I([C i]) and N(H 2 ) (Fig. 10), which is in sharp contrast to the tight correlation between I([C i]) and I( 13 CO) (Fig.  5 right). While bearing in mind possible sources for the scatter (e.g., variations in the assumed properties of dust and gas, such as the dust opacity at 870 µm and dust-to-gas mass ratio), we take our result at face value and interpret it as a result of [C i] 492 GHz being a tracer of moderately dense molecular gas that is more closely associated with 13 CO(2-1). Indeed, Duarte-Cabral et al. (2021) recently performed a cross-match between the ATLASGAL 870 µm and SEDIGISM 13 CO(2-1) surveys and found that 4824 ATLASGAL clumps are enveloped within 1709 SEDIGISM clouds, implying that 13 CO(2-1) probes the molecular gas that surrounds denser and colder clumps traced by 870 µm. Similarly, Plume et al. (2000) and Shimajiri et al. (2013) found that [C i] 492 GHz is poorly correlated with CS(1-0) and C 18 O(1-0) (dense molecular gas tracers) in Orion A. In general, our finding suggests that [C i] 492 GHz would not be a good tracer of total gas mass if the cloud mass is dominated by dense and cold molecular gas.

Physical conditions traced by [C i]
Two fine-structure transitions of [C i] at 492 GHz and 809 GHz are invaluable to constrain the physical properties of gas. For instance, in the optically thin regime, their flux ratio constrains the excitation temperature (T ex ) as the integrated intensity is proportional to the upper state column density. In more general situations, including high optical depth cases, the physical properties of gas can be derived through escape probability radiative trans- fer models such as RADEX (van der Tak et al. 2007). In this section, we examine the physical conditions of [C i]-traced gas in a small sub-set of the Top100 sources by analyzing the [C i] 492 GHz and 809 GHz observations in the commonly adopted approximation of local thermodynamic equilibrium (LTE). Considering that what kind of the interstellar medium [C i] traces is not entirely clear, we perform a non-LTE calculation as well to evaluate the impact of sub-critical densities.  Table 2.

LTE calculation
Given that the observed sources are dense clumps, we could assume LTE conditions and derive physical conditions including the excitation temperature, opacities at 492 GHz and 809 GHz (τ 1−0 and τ 2−1 ), and atomic carbon column density (N(C)) (Appendix C). But before the LTE calculation, we first explored how the ratio of the 809 GHz to 492 GHz peak main-beam brightenss temperature (T 2−1 /T 1−0 ) changes with the excitation temperature in two extreme cases: optically thin and thick regimes. When τ ≪ 1, Eqs. (C.1) and (C.5) suggest that T 2−1 /T 1−0 can be written as On the other hand, when τ ≫ 1, T 2−1 /T 1−0 can be written as Eqs.
(1) and (2) are plotted as a function of the excitation temperature in Fig. 11.   In the case of G338.78+0.48, the measured T peak,1−0 is most likely a lower limit, since its [C i] 492 GHz spectrum is affected by a contaminated reference position. would be about half as bright when the gas is cool (T ex ∼ 10-20 K) and optically thick. The observed T 2−1 /T 1−0 ranges from ∼0.9 to ∼1.7 ( Table 2), suggesting that the [C i]-emitting gas is likely warm and optically thin.
To estimate T ex , τ 1−0 , τ 2−1 , N(C), and their associated 1σ uncertainties, we performed Monte Carlo (MC) simulations where the [C i] spectra are perturbed 1000 times based on the rms values (Gaussian errors are assumed). Each perturbed spectrum is fitted with a Gaussian function to measure the peak temperature, and Eqs. (C.2), (C.3), and (C.8) are applied to the measured temperatures to derive T ex , τ 1−0 , τ 2−1 , and N(C). For our calculation of N(C), τ 1−0 is used since it has a higher sensitivity than τ 2−1 . Finally, the distribution of each parameter is examined to estimate a median value, as well as ±34% boundaries from the median (so that the data points within the boundaries are 68% of the total number). These median values and 1σ boundaries are provided in Table 2, and the normalized histograms of the derived parameters for G345.49+0.32 and G340.00−0.23 are shown in Fig. 12 to demonstrate our calculation. We note that G345.00−0.23 has asymmetric distributions, leading to a large uncertainty in T ex in particular. This is because the measured temperature ratio of 1.7 is on the asymptotic part of the T 2−1 /T 1−0 Article number, page 9 of 27 A&A proofs: manuscript no. draft_040622 Fig. 13. Predicted 809 GHz-to-492 GHz temperature ratios for N(C) = 10 16 , 10 17 , and 10 18 cm −2 (left, middle, and right respectively). In each plot, the overlaied contours range from 10% to 90% of the peak temperature ratio with 20% steps. versus T ex curve (Fig. 11) and therefore a small change in the ratio can result in a large variation in T ex . For the five sources where [C i] 809 GHz is not detected, the lower and upper limits on the physical parameters are provided in Table 2. Zmuidzinas et al. (1988) observed Galactic star-forming regions (OMC 1, NGC 2024, S140, DR 21(OH), W3, M17, and W51) in [C i] 492 GHz and 809 GHz and performed similar excitation analyses. Compared to Zmuidzinas et al. (1988) results (32-77 K), our excitation temperatures (49-186 K) are slightly higher. In particular, we found a tendency of higher temperatures for more evolved sources (94 K and 186 K for the two H ii region sources), though this needs to be verified by extensive [C i] 809 GHz observations in the future. Our optical depth values are, on the other hand, lower than Zmuidzinas et al. (1988), i.e., τ 1−0 ∼ 0.01-0.32 versus 0.22-1.25.

Non-LTE consideration
While there is certainly a significant amount of dense gas in the observed massive clumps, whether the [C i] lines arise from this dense gas or from the lower density gas that surrounds the dense clumps is currently not clear. Since our previous calculations are based on the LTE assumption, it is therefore important to examine the impact of sub-critical densities on the derived parameters.
To this end, we employed RADEX (van der Tak et al. 2007), where the radiative transfer equation for isothermal and homogeneous gas is solved based on the escape probability method to compute the intensities of atomic and molecular transitions. To cover a wide range of physical conditions, we created a RADEX model grid with the following parameters: kinetic temperature T k = 10-500 K, H 2 density n(H 2 ) = 10 2 -10 6 cm −3 , and C 0 column density N(C) = 10 16 -10 18 cm −2 (uniformly sampled in log space with ten, seventeen, and nine points). In addition to these key parameters, the line width of 5 km s −1 (median of the measured FWHM 1−0 ; Table 2) and the background radiation of 2.73 K were assumed. To illustrate the RADEX modeling results, we present the predicted T 2−1 /T 1−0 for N(C) = 10 16 , 10 17 , and 10 18 cm −2 in Fig. 13. Fig. 13 illustrates a well-known degeneracy between the temperature and density: high density/low temperature and low density/high temperature combinations produce the same line ratio.
To reproduce the observed T 2−1 /T 1−0 of 0.9-1.7, T k 60 K and n(H 2 ) 10 3 cm −3 are required, mostly ensuring LTE (densities lower than ∼10 4 cm −3 mainly result in a difference between the kinetic and excitation temperatures). To explore the RADEX predictions more closely, we then selected the N(C) values that are closest to our LTE-based estimates and extracted the physical parameters that reproduce the observed line ratios within 1σ uncertainties (Table 3). Our results show that the H 2 density, kinetic temperature, and thermal pressure are poorly constrained, while the optical depths are relatively well constrained and not too different from our LTE values. The large variations in T k and n(H 2 ) are mostly driven by the density-temperature degeneracy and hinder us from inferring any systematic variation across the different evolutionary groups. On the other hand, the predicted brightness temperatures are comparable to the observed values. All these results imply that our LTE-based optical depths and C 0 column densities are not heavily affected by sub-critical densities and [C i] emission arises from warm (T k 60 K), optically thin (τ < 0.5), and highly pressurized regions (P/k B ∼ a few (10 5 -10 8 ) K cm −3 ) in several of the Top100 sources.

Summary
In this paper, we present APEX observations of [C i] 492 GHz toward ATLASGAL-selected massive clumps in the inner Galaxy (Top100 sample). Our target sources have been extensively examined in dust continuum and spectral lines and constitute a representative sample of high-mass star-forming regions covering a significant range of masses (∼20-10 4 M ⊙ ), bolometric luminosities (∼60-10 6 L ⊙ ), and evolutionary phases (70 µm weak, IR-weak, IR-bright, and H ii regions). To examine the physical conditions of [C i]-traced gas in the Top100 sources, we combined our APEX [C i] 492 GHz spectra with multi-wavelength observations of atomic and molecular gas and obtained the following results.  Notes. (a) The observed temperature ratio and its 1σ uncertainty (in parentheses). (b) The input N(C) for RADEX that is closest to the LTE-based N(C). (c) The RADEX parameters that are found to reproduce the observed temperature ratios within 1σ uncertainties (P/k B = thermal pressure). The ranges, median values, and 1σ uncertainties are presented in rows.
GHz is prevalent in the inner Galaxy and traces not only massive clumps, but also non-star-forming relatively diffuse gas. GHz diffuse cloud components without corresponding 13 CO(2-1) emission contain a significant amount of atomic gas and could trace "CO-dark" molecular gas. This conclusion, however, is based on few lines of sight and needs to be verified by sensitive [C i] 492 GHz, 13 CO(2-1) and 12 CO(1-0) observations in the future. 5. We derived the X(C i) factors by dividing the 870 µm-based H 2 column densities by the [C i] 492 GHz integrated inten-sities and found that X(C i) (in units of cm −2 (K km s −1 ) −1 ) ranges from 2.3 × 10 20 to 1.3 × 10 22 with a median of 1.7 × 10 21 . Our median value is in good agreement with recent observations of the massive star-forming region RCW 38 (1.4 × 10 21 ), as well as predictions from numerical simulations of Galactic molecular clouds (∼10 21 ). 6. In contrast to the tight correlation with 13 CO(2-1), GHz shows a much more scattered relation with the 870 µmtraced molecular gas. We interpreted this as [C i] 492 GHz and 13 CO(2-1) probing warm molecular gas that surrounds denser and colder clumps traced by 870 µm and cautioned not to use [C i] 492 GHz as a tracer of total gas mass if dense and cold molecular gas dominates the cloud mass. 7. our LTE and non-LTE analyses showed that [C i] emission arises from warm (T k 60 K), optically thin (τ < 0.5), and highly pressurized regions (P/k B ∼ a few (10 5 -10 8 ) K cm −3 ) in several of the Top100 sources.
Based on a large sample of sensitive [C i] 492 GHz spectra, this paper investigated the physical properties of [C i]-traced gas in a wide range of environments and shedded light on the origin of [C i] 492 GHz emission by probing the relation with molecular gas tracers. One of our key results is the systematic variation in the ratio of 13 CO(2-1) to [C i] 492 GHz, which implies that the Top100 sources are in different stages of carbon phase transition (C + -, C 0 -, or CO-dominated). We will examine the exact stages of carbon phase transition and the underlying physical and energetic conditions in forthcoming papers Article number, page 11 of 27 A&A proofs: manuscript no. draft_040622 where [C ii], [C i], [O i], and CO observations from our multitelescope campaign are analyzed using state-of-the-art models of photodissociation regions and shocks (e.g., Le Petit et al. 2006;Flower & Pineau des Forêts 2015).  Notes. The physical characteristics are from König et al. (2017), and the columns are as follows. Name: ATLASGAL source name; l: Galactic longtidue; b: Galactic latitude; d: distance; src : source velocity; log 10 (N(H 2 )): dust-based N(H 2 ) in log; Class: source classification; [C i] 809 GHz?: presence of the [C i] 809 GHz mapping data  Notes. These are the data that were used for the analyses in Sect. 4.1. Among all fitted components, only those with clean reference positions were selected. When one of the two lines is not detected, 5σ of the undetected transition is integrated over the velocity range where the other transition is seen, providing a upper limit on the integrated intensity (indicated as the downward arrows). The columns are as follows. Name: ATLASGAL source name; c : representative velocity of the component (determined based on Gaussian fitting); l : lower velocity for emission integration; u : upper velocity for emission integration; I: integrated intensity; σ(I): 1σ on the integrated intensity.

Appendix C: LTE calculation
To calculate the physical properties of the [C i]-traced main components (Sect. 5.1), we combined the APEX [C i] 492 GHz and 809 GHz observations and followed the LTE-based procedure as below. First, we started with the definition of optical depth: where c is the speed of light, ν is the frequency of the transition, h is the Planck constant, k B is the Boltzmann constant, A ul is Einstein's spontaneous emission coefficient, φ is the normalized spectral distribution of the transition as a function of velocity, and N u is the column density of atoms in the upper transition state. With this definition of optical depth, the τ 2−1 -to-τ 1−0 ratio can be written as where g u and E u are the statistical weight and energy of the upper transition state. The partition function Z is defined as