C i and CO in nearby spiral galaxies I. Line ratio and abundance variations at ∼ 200 pc scales

We present new neutral atomic carbon [C i ]( 3 P 1 → 3 P 0 ) mapping observations within the inner ∼ 7kpc and ∼ 4kpc of the disks of NGC3627 and NGC4321 at a spatial resolution of 190pc and 270pc, respectively, using the Atacama Large Millimeter / Submillimeter Array (ALMA) Atacama Compact Array (ACA). We combine these with the CO(2 − 1) data from PHANGS-ALMA, and literature [C i ] and CO data for two other starburst and / or active galactic nucleus (AGN) galaxies (NGC1808, NGC7469) with the aim of studying: (a) the spatial distributions of C i and CO emission; (b) the observed line ratio R C i / CO = I [C i ](1 − 0) / I CO(2 − 1) as a function of various galactic properties; and (c) the abundance ratio of [C i / CO]. We ﬁnd excellent spatial correspondence between C i and CO emission and nearly uniform R C i / CO ∼ 0 . 1 across the majority of the star-forming disks of NGC3627 and NGC4321. However, R C i / CO strongly varies from ∼ 0 . 05 at the center of NGC4321 to > 0 . 2 − 0 . 5 in NGC1808’s starbursting center and NGC7469’s center with an X-ray-luminous AGN. Meanwhile, R C i / CO does not obviously vary with (cid:104) U (cid:105) , which is in line with predictions from photodissociation-dominated region (PDR) models. We also ﬁnd a mildly decreasing R C i / CO value with an increasing metallicity over 0 . 7 − 0 . 85 Z (cid:12) , which is consistent with the literature. Assuming various typical interstellar medium (ISM) conditions representing giant molecular clouds, active star-forming regions, and strong starbursting environments, we calculated the (non)local-thermodynamic-equilibrium radiative transfer and estimated the [C i / CO] abundance ratio to be ∼ 0 . 1 across the disks of NGC3627 and NGC4321, similar to previous large-scale ﬁndings in Galactic studies. However, this abundance ratio likely experiences a substantial increase, up to ∼ 1 and (cid:38) 1 − 5 in NGC1808’s starburst and NGC7469’s strong AGN environments, respectively. This result is in line with the expectations for cosmic-ray dominated region (CRDR) and X-ray dominated region (XDR) chemistry. Finally, we do not ﬁnd robust evidence for a generally CO-dark-and-C i -bright gas in the disk areas we probed.


Neutral atomic carbon in the interstellar medium
Neutral atomic carbon (C 0 or C i) is an important phase of carbon in the interstellar medium (ISM) besides ionized carbon (C + or C ii) and carbon monoxide (CO).Its 3 P state is split into three fine-structure levels.The 3 P 1 → 3 P 0 transition line at 492.16065 GHz, hereafter [C i] (1-0), has an upper energy level of E u = 23.620K and a critical density of n crit ∼ 1.0 × 10 3 cm −3 (at a temperature of 50 K) 1 , and the 3 P 2 → 3 P 1 transition line at 809.34197 GHz has E u = 62.462K and n crit ∼ 3.4 × 10 3 cm −3 .Therefore, they are easily excited in the cold ISM environment.Given the high abundance of carbon in the ISM, C i, CO and C ii are the most powerful cold ISM tracers, and are the most widely used tools to probe cold ISM properties and hence galaxy evolution at cosmological distances.
The interstellar C i is produced from the photodissociation of CO molecules by UV photons and the recombination of C ii, both Deceased. 1 n crit ≡ A/γ, where A is the Einstein A coefficient for the transition and γ the rate of collisional de-excitation with hydrogen molecules, with values from Schroder et al. (1991) as compiled in the Leiden Atomic and Molecular Database (Schöier et al. 2005), https://home.strw.leidenuniv.nl/~moldata/happening in photodissociation regions (PDRs;Langer 1976;de Jong et al. 1980;Tielens & Hollenbach 1985a,b;van Dishoeck & Black 1986, 1988;Sternberg & Dalgarno 1989, 1995;Hollenbach et al. 1991;Hollenbach & Tielens 1999;Kaufman et al. 1999;Wolfire et al. 2010;Madden et al. 2020;Bisbas et al. 2021, also see reviews by Genzel & Stutzki 1989;Jaffe et al. 1985;Hollenbach & Tielens 1997;Wolfire et al. 2022).In the simple plane-parallel (i.e., 1-D) PDR model, C i exists in a layer where UV photons can penetrate through the molecular gas but cannot maintain a high carbon ionizing rate.Such a layer is suggested to have an intermediate gas surface density (N H 2 ∼ 1-4 × 10 20 cm −2 ), fairly cold temperature (10 to a few tens K), and moderate visual extinction (A V ∼ 1−4).CO molecules dominate the more shielded interior of this layer and C ii dominates the exterior.Hydrogen molecules are gradually photodissociated to atoms across this layer (A V ∼ 0.1 − 4).This plane-parallel scenario qualitatively explains the layered structure of the Orion Bar PDR in our Galaxy on sub-pc scales (e.g., Genzel & Stutzki 1989;Tielens et al. 1993;Tauber et al. 1994;Hogerheijde et al. 1995;Hollenbach & Tielens 1999;Goicoechea et al. 2016), and the ρ Ophiuchi A PDR in a recent study with an unprecedented resolution of 360 AU (Yamagishi et al. 2021), where the interplay between individual young, massive (O, B) stars and the ISM can be directly observed.
However, such simple models cannot fully explain molecular clouds.Plenty of observations in Galactic molecular clouds, H ii regions and star-formation complexes have revealed a widespread C i distribution over most areas of molecular clouds on scales of a few to ten pc, and a relatively high C i fractional abundance even at a large A V depth > 10 mag (Phillips et al. 1980;Phillips & Huggins 1981;Wootten et al. 1982;Keene et al. 1985Keene et al. , 1987;;Jaffe et al. 1985;Zmuidzinas et al. 1986Zmuidzinas et al. , 1988;;Genzel et al. 1988;Frerking et al. 1989;Plume et al. 1994Plume et al. , 1999;;Gerin et al. 1998;Tatematsu et al. 1999;Ikeda et al. 1999Ikeda et al. , 2002;;Yamamoto et al. 2001;Zhang et al. 2001;Kamegai et al. 2003;Oka et al. 2004;Izumi et al. 2021).
This observational evidence actually requires PDRs to be clumpy (Stutzki et al. 1988;Genzel et al. 1988;Burton et al. 1990;Meixner & Tielens 1993;Spaans & van Dishoeck 1997;Kramer et al. 2004Kramer et al. , 2008;;Pineda et al. 2008;Sun et al. 2008), with a volume filling factor much less than unity (e.g., ∼ 0.1 − 0.3; Stutzki et al. 1988), so that UV photons can penetrate through most of the cloud, except for the densest clumps.Therefore, the spatial co-existence of C i and low rotational transition (low-J) CO lines can be well explained by such an inter-clump medium inside the molecular clouds.Alternatively, cosmic rays (CRs) have also been proposed to be the reason for the C i-CO association as they can penetrate much deeper than UV photons into clouds, and thus more uniformly dissociate CO molecules into C i inside clouds (Field et al. 1969;Padovani et al. 2009;Papadopoulos 2010;Ivlev et al. 2015;Papadopoulos et al. 2018).
At sub-galactic scales, when individual clouds are smoothed out, R Ci10/CO10 is found to lie in a narrower range of ∼ 0.1-0.3.For example, across the Galactic plane, the ratio is found to change only slightly within a galactocentric radius of 3-7 kpc (∼ 0.08-0.12,overall mean = 0.105 ± 0.004; Oka et al. 2005); whereas in the starburst galaxies like M 82 and NGC 253, R Ci10/CO10 or the [C i]/CO abundance ratio is a factor of 3-5 higher than the Galactic value (Schilke et al. 1993;White et al. 1994;Krips et al. 2016).Recent highresolution (hundred pc scales) [C i] (1-0) mapping studies in the nearby spiral galaxy M 83 (Miyamoto et al. 2021) and starburst 2 The line intensity I has brightness temperature units K km s −1 .
Meanwhile, in the central 100 pc of NGC 7469 close to the X-ray-luminous active galactic nucleus (AGN), R Ci/CO is found to be substantially enhanced by a factor of > 10 (reaching R Ci10/CO10 ≈ 1; Izumi et al. 2020; see also the Circinus Galaxy's AGN with R Ci10/CO32 ≈ 0.9 in Izumi et al. 2018).
Furthermore, in the low-metallicity Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC) environments, a factor of 1.5 − 3× higher R Ci/CO values have been found (Bolatto et al. 2000b,a; see also other low-metallicity galaxies : Bolatto et al. 2000c;Hunt et al. 2017).
These studies indicate that R Ci/CO is sensitive to ISM conditions among and within galaxies.However, previous extragalactic studies are mainly focused on starburst galaxies and bright galaxy centres.How R Ci/CO behaves across the disks and spiral arms of typical star-forming galaxies is much less explored.Whether C i can trace the so-called CO-dark molecular gas at a low gas density or metallicity is also still in question.Understanding R Ci/CO variation and quantifying the C i and CO excitation in extragalactic environments is now particularly important given the rapidly growing numbers and diverse galactic environments of C i line observations at cosmological distances (e.g., Weiß et al. 2003Weiß et al. , 2005;;Pety et al. 2004;Wagg et al. 2006;Walter et al. 2011;Danielson et al. 2011;Alaghband-Zadeh et al. 2013;Gullberg et al. 2016;Bothwell et al. 2017;Popping et al. 2017;Yang et al. 2017;Emonts et al. 2018;Valentino et al. 2018Valentino et al. , 2020;;Bourne et al. 2019;Nesvadba et al. 2019;Dannerbauer et al. 2019;Jin et al. 2019;Brisbin et al. 2019;Boogaard et al. 2020;Cortzen et al. 2020;Harrington et al. 2021;Dunne et al. 2021;Lee et al. 2021).
1.3.What drives R Ci/CO ?-the need for high-resolution mapping across wide galactic environments High-resolution mapping of C i and CO from disks and spiral arms of nearby galaxies is indispensable for understanding and quantifying the ISM C i and CO physics.In particular, the following questions will be answered only with such observations: Do C i and low-J CO have the same spatial distributions across galaxy disks?Does the R Ci/CO line ratio vary significantly across galaxy disks?How does the R Ci/CO line ratio vary among different local galaxies?How can the observed R Ci/CO line ratio be reliably converted to the C i and CO column density ratio (i.e., solving their radiation transfer equations and level population statistics and obtaining optical depths and excitation temperature)?
Currently, high-resolution (hundreds of parsec scale) imaging/mapping of C i in nearby galaxies relies on radio telescopes at high altitudes, e.g., the Atacama Large Millimeter/submillimeter Array (ALMA), Atacama Submillimeter Telescope Experiment (ASTE, 10-m single-dish), Atacama Pathfinder Experiment (APEX, 12-m single-dish), and is limited to the available weather conditions corresponding to the high frequency of C i lines.Therefore, only a handful of galaxies have been imaged/mapped.These include NGC 6240 by Cicone et al. (2018), Circinus Galaxy by Izumi et al. (2018), NGC 613 by Miyamoto et al. (2018), NGC 1808 by Salak et al. (2019), NGC 7469 by Izumi et al. (2020), IRAS F18293-3413 by Saito et al. (2020), NGC 6052 by Michiyama et al. (2020, C i undetected), 36 (U)LIRGs by Michiyama et al. (2021), M 83 by Miyamoto et al. (2021, mapped the centre and northern part with ASTE), Arp 220 by Ueda et al. (2022), and NGC 1068 by Saito et al. (2022a,b).Nevertheless, all of these observations, except for M 83, targeted the centres of starburst, IR-luminous and/or AGN host galaxies.M 83 is the only typical star-forming main sequence galaxy having a high-quality C i mapping.However, its map is still limited to a quarter of its inner ∼3 kpc disk.
In this work, we present two new C i mapping observations with the ALMA Atacama Compact Array (or Morita Array; 7-m dish, hereafter ACA or 7m) in nearby spiral galaxies NGC 3627 and NGC 4321.Our observations aimed at detecting C i in the disks, and thus are much larger-area and deeper than previous C i maps.By further adding a starburst and an AGN-host galaxy from the literature into our analysis, we present a comprehensive study of the R Ci10/CO21 line ratio (hereafter R Ci/CO ; ratios of other C i/CO transitions will be explicitly written otherwise) in nearby star-forming disk and starburst and AGN environments.In order to quantify the C i and CO abundance ratio, we have further calculated the local thermodynamic equilibrium (LTE) and non-LTE radiative transfer of C i and CO under several representative ISM conditions.The results from this study could thus be one of the most comprehensive local benchmarks for understanding the C i and CO line ratio and abundance variations.
The paper, as the first one of our series of papers, is organized as follows.Sect. 2 describes the sample selection, observations and data reduction.Sect. 3 presents the spatial distribution and variation of R Ci/CO .Sect. 4 presents the analysis of C i and CO excitation and radiative transfer and thus the abundance (species column density) ratio.In Sect. 5 we discuss various topics relating to the [C i/CO] abundance variation, CRDR, XDR, and COdark gas.Finally we conclude in Sect.6.Our companion paper II (D. Liu et al. in prep.)presents the detailed radiative transfer calculation and C i and CO conversion factors, and is related to the abundance calculation in this paper.

Target selection
Our targets, NGC 3627 and NGC 4321, are selected from a joint sample of the PHANGS-ALMA CO(2-1) (Leroy et al. 2021b), PHANGS-HST (Lee et al. 2022) The PHANGS-ALMA sample consists of 90 nearby starforming galaxies initially selected at distances of 2 d/Mpc 23, with a stellar mass log(M /M ) > 9.75, and which are not too inclined and visible to ALMA.They represent the typical star-forming main sequence galaxies in the local Universe.The PHANGS-HST sample (Lee et al. 2022) is a subsample of 38 PHANGS-ALMA galaxies, providing high-resolution stellar properties.The PHANGS-MUSE large program (Emsellem et al. 2022) targeted a subsample of 19 PHANGS-ALMA galaxies suitable for VLT/MUSE optical integral field unit (IFU) observations, achieving similar spatial coverage as in PHANGS-ALMA and providing rich nebular emission lines, attenuation, stellar age and metallicity and H ii region information (e.g., Kreckel et al. 2019Kreckel et al. , 2020;;Pessa et al. 2021;Santoro et al. 2022;Williams et al. 2022).
For the legacy value, we selected our sample from the joint ALMA+MUSE+HST sample pool, considered CO surface brightness, mapping area, C i line frequency and the transmission in the ALMA Band 8.In this work, we present the new C i data of NGC 3627 and NGC 4321.Because of the ALMA Band 8 sensitivity and the expected fainter C i line intensity than CO, the C i mapping areas are smaller than their CO(2-1) maps, but they are still the largest (∼ 1 − 2.35 ) and deepest (rms ∼ 0.04 K per ∼ 5 km s −1 ) C i mapping with ALMA.
Furthermore, in order to study the C i and CO in different galactic environments, we selected two starbursting and/or AGN-host galaxies from the literature whose ALMA [C i] (1-0) and CO  data are available in the archive: NGC 1808 (Salak et al. 2019) and NGC 7469 (Izumi et al. 2020).NGC 1808 has ALMA mosaic observations covering its ∼ 1 kpc starburst ring, and NGC 7469 has ALMA single-pointing observations covering its inner ∼ 3 kpc area.

New ALMA ACA Band 8 observations
Our C i observations for NGC 3627 were taken under the ALMA project 2018.1.01290.S (PI: D. Liu) during June 3rd to Sept. 29th, 2019 with ACA-only at Band 8.The total on-source integration time was 30.8 hours.A mosaic with 149 pointings was used to map an area of 94 × 136 at a position angle of −27 • , which covers the galaxy centre, bar and the majority of the spiral arms, and is about 1/3 of the full PHANGS-ALMA CO(2-1) map area.The achieved line sensitivity per 5 km s −1 is 45-55 mK (107-130 mJy/beam) across the [C i] (1-0) line frequencies with a beam of 3.46 .
The C i mapping for NGC 4321 were taken under the ALMA project 2019.1.01635.S (PI: D. Liu; NGC 1365 is also partially observed for this project) during Dec. 3rd, 2019 to July 1st, 2021, also with ACA-only at Band 8.The on-source integration time reached 16.1 hours for the mapping of an area of 55 × 42 with 23 pointings using the 7m array.The achieved line sensitivity is similar to that of NGC 3627.More detailed information is provided in Table 1.
No total power observations were requested in our presented observation, therefore we matched the uv ranges of the raw CO(2-1) and [C i] (1-0) visibility data during our data reduction.We verified that missing flux will not substantially (∼ 30%) bias our analysis using three methods as detailed in Appendix B. This includes a CASA simulation of visibilities mimicking our ACA mosaic observations in NGC 3627, from which we find a missing flux as small as 7% in C i-bright pixels and up to ∼ 50% in the faintest pixels (and mostly < 30%).
Our data reduction follows the PHANGS-ALMA imaging and post-processing pipeline (Leroy et al. 2021a), including imaging and deconvolution, mosaicking, broad and strict/signal mask generation, and moment map creation.The additional steps in this study are: (a) uv-clipping in the CO(2-1) data to match the uv sampling range of [C i] (1-0); and (b) making a joint mask of the individual CO(2-1) and [C i] (1-0) masks then extracting the moment maps (in the joint broad mask) and line ratios maps (in the joint strict mask).More technical details are provided in Appendix A.
Our [C i] (1-0) and CO(2-1) line intensity maps, the line ratio maps and HST images for the two PHANGS galaxies, NGC 3627 and NGC 4321, are presented in Figs. 1 and 2, respectively.Our data cubes and moment maps are also made publicly available at the PHANGS-ALMA CADC public repository3 .

ALMA observations of the additional targets from the literature
The NGC 1808 ALMA [C i] (1-0), CO (1-0) and CO (2-1) observations are presented in Salak et al. (2019) and all are interferometry plus total power, i.e., 12m+7m+tp, thus having no missing flux issue.To study the C i/CO line ratio consistently with our main targets, we only used the [C i] (1-0) and CO (2-1) data.We re-reduced the raw data under ALMA project 2017.1.00984.S (PI: D. Salak) with the observatory calibration pipeline, then imaged, primary-beam-corrected and short-spacing-corrected using the PHANGS-ALMA pipeline in the same way as for our main targets (Appendix A).Note that in this step we re-reduced and imaged the total power data with the PHANGS-ALMA singledish pipeline (Herrera et al. 2020;Leroy et al. 2021a) and did the short-spacing correction with the PHANGS-ALMA pipeline right after the primary beam correction.The properties of the final CO-C i beam-matched data cubes are listed in Table 1.
The NGC 7469 ALMA data are presented in Izumi et al. (2020) and we used the archival raw data of [C i] (1-0) and CO (2-1) under ALMA project 2017.1.00078.S (PI: T. Izumi).These are single pointings toward the centre.[C i] (1-0) is observed with ACA 7m and CO (2-1) with the 12m array, both without total power.We re-reduced the data following the same procedure as for NGC 3627 and NGC 4321 with the PHANGS-ALMA pipeline.
Our re-reduced [C i] (1-0) and CO(2-1) line intensity and ratio maps for these additional galaxies are shown in Fig. 3.

Spatial distributions
In Figs. 1, 2 and 3, we present the beam-matched [C i] (1-0) and CO (2-1) integrated line intensity maps and their ratio maps.The C i/CO line ratio map is computed only for pixels within the joint strict mask area of the two lines, as described in Appendix A, representing high-confidence signals.
In NGC 3627 (Fig. 1), we visually marked 23 regions based on the C i and CO commonly-detected (S/N > 3), strict-mask areas for later analysis.Region IDs are sorted by their [C i] (1-0) brightness.Region 1 corresponds to the galaxy centre, Regions 2 and 3 are the southern and northern bar ends, Region 4 is an offset peak next to the galaxy centre, and other regions are mostly on the spiral arms.Most of these regions show a similar R Ci/CO ∼ 0.1, except for, e.g., Regions 20 and 22 which exhibit a somewhat higher ratio of R Ci/CO ∼ 0.5.The edge of Region 4 also shows an enhancement of C i yet it is hard to tell whether this is physical or simply caused by edge effects, e.g., as seen in our simulations (Appendix C).
We also examined the integrated spectra of each region in Appendix D. The C i and CO line profiles are not always wellmatched in width and shape.In the galaxy centre (Region 1) and at bright bar-ends (Regions 2, 3, 5), C i and CO have similar line widths but C i has a lower intensity than CO, whereas in some fainter regions (9,10,13,14), C i is narrower than CO and relatively weaker.In Regions 4, 20 and 22 where R Ci/CO is high, C i and CO have similar line widths but the C i intensity is enhanced.
In NGC 4321 (Fig. 2), we mapped a smaller area due to its fainter CO brightness than NGC 3627.There are no regions with CO (2-1) integrated intensities > 100 K km s −1 outside our selected area in NGC 4321, whereas in comparison the NGC 3627 bar ends have CO (2-1) integrated intensities 300 K km s −1 .The C i mapping in NGC 4321 has a similar depth as in NGC 3627.
We marked 7 regions in NGC 4321 to guide the visual inspection.A clear R Ci/CO deficit can be seen in the galaxy centre (Region 1), with R Ci/CO ∼ 0.06, nearly a factor of two smaller than its inner disk which has similar R Ci/CO ∼ 0.1 as seen in the NGC 3627 disk.Some C i enhancement can be found at the edges or slightly outside these regions, where CO and C i are both faint (I CO (2-1) ∼ 30 − 100 K km s −1 ) and the edge effect of interferometric missing flux might play a role (see Appendix C).
The R Ci/CO spatial distributions of the other two starburst/AGN galaxies (Fig. 3) show in general more enhanced R Ci/CO than in the NGC 3627 and NGC 4321 disks.Large variations from the centres to outer areas can be seen, especially at the the S/N threshold of ∼ 3-5 (see contours).
To quantify the spatial co-localization of the C i and CO emission, we use the JACoP tool (Bolte & Cordelières 2006) 4  in the ImageJ software (Schneider et al. 2012) to measure the commonly used co-localization indicators, namely the Pearson's coefficient (PC; equal to 0 and 1 for null and 100% localization, respectively) and Manders' coefficients (M1 and M2; which measure the overlap between two images and are equal to one for perfect overlap).We turn on the option of using Costes' automatic threshold (Costes et al. 2004) which automatically set a pixel value threshold for each image to minimize the contribution of noise for the comparison (Bolte & Cordelières 2006).For the C i and CO emission in our galaxies, we find PC/M1/M2 values all very close to 1, that is, nearly 100% colocalization (PC/M1/M2 = 0.934/0.985/0.990for NGC 3627, 0.816/0.986/0.964for NGC 4321, 0.978/0.890/0.906for NGC 1808, and 0.935/0.999/0.998for NGC 7469).
Overall, the spatial comparison of the C i and CO intensities demonstrates that the two species are well-correlated when observed at hundred parsec scales.Meanwhile, we do see certain spectral differences in the line profiles as shown in Appendix D. This may suggest that variations exist at smaller scales than our resolution, which could be the sign of different evolutionary stages (e.g.Kruijssen & Longmore 2014) and would be interesting for higher-resolution follow-ups, e.g., with the main ALMA 12m array.

R Ci/CO versus CO surface brightness and gas surface density
In Fig. 4 we show the scatter points of C i and CO line brightnesses and the R Ci/CO line ratios.Systematic differences in R Ci/CO can be found among these galaxies and within their galactic environments.We discuss these environments separately as follows.

Star-forming and starburst disk environments
The majority of molecular gas in the 3-7 kpc star-forming disks of NGC 3627 and NGC 4321 exhibits a tight distribution of R Ci/CO around 0.1 (within about ±50%).The distribution is fairly flat over an order of magnitude in I CO (2-1) .We denote this as the "SF disk" regime.The ratio rises to about 0.2 within a galactocentric radius of few hundred parsec to ∼ 1 kpc starburst areas in NGC 1808 and NGC 7469, which we define as the "SB disk" regime.
The nearly flat "SB disk" distribution is consistent with earlier studies.Salak et al. (2019) found that in NGC 1808 CO  .They also studied CO(1-0), finding that CO (1-0) , i.e., with a steeper slope.Saito et al. (2020) found a similar slope between CO (1-0) in IRAS F18293-3413 (whose CO(2-1) data is insufficient for such a study as previously mentioned).This correlation slope is likely to change with different CO transitions due to the CO excitation's dependence on the CO luminosity or star formation rate themselves (e.g., Daddi et al. 2015;Liu et al. 2021).For higher-J CO transitions, Michiyama et al. (2021) studied C i and CO(4-3) lines in 36 (U)LIRGs, finding CO (4-3) .The upward tails seen at the CO-faint-end in the right panel of Fig. 4 are mostly due to the detection limits in each galaxy.The color shading indicates the 5-σ limit for each galaxy, assuming a typical line width of 10 km s −1 .Region 22 of NGC 3627 and the outer pixels in the two starburst/AGN galaxies make the tips of the blue and gray faint-end tails, respectively.But most of these features should still be spurious.We performed CASA simulations following the real observations of NGC 3627 to verify what we would get if R Ci/CO were constant across the galaxies (Appendix C) We find that noise can indeed boost a simulated R Ci/CO of 0.2 to a measurement as high as 0.5 in a NGC 3627-like large mosaic, or even higher to > 1 in single-pointing ALMA data without total power or sufficient short spacing.However, the faint-end tail from the simulation is not very similar to that of NGC 3627, raising some questions about the nature of the high-R Ci/CO in Region 22.
Furthermore, we also did stacking experiments using the CO line profile as the prior to measure the C i line fluxes in bins of CO line brightnesses or galactocentric radii.We find no systematically higher R Ci/CO at fainter CO pixels or at the outer radii.

Galaxy centre environments
Most interestingly, the 250 pc centres of all four galaxies show dramatically different R Ci/CO .It varies by a factor of ten from ∼ 0.05 in the centre of NGC 4321 to ∼ 0.1 in the centre of NGC 3627, then to ∼ 0.15-0.2 at the centre of NGC 1808, and to ∼ 0.2-0.5 in the nucleus of NGC 7469.This "nuclei" trend has a slope of about 0.8, i.e., CO (2-1) , and corresponds to a Σ mol range of about 10 3 to 10 4 M pc −2 if adopting an α CO = 4.3 M (K km s −1 pc 2 ) −1 (Bolatto et al. 2013) and a CO excitation R 21 = 0.65 (Leroy et al. 2022).
The reasons for the nuclei showing this trend likely include both the effect of varying [C i/CO] abundance ratio and the excitation and radiative transfer.They are determined by the actual ISM conditions in each environment, i.e., the gas kinetic temperature, C i and CO column densities (N Ci and N CO ), line widths (∆v), and H 2 volume density (n H 2 ), etc.They will be analyzed in Sect. 4.

Histogram and curve-of-growth of the line ratio distributions
Fig. 5 shows the histogram and cumulative distribution functions of R Ci/CO (left panel) and the mean R Ci/CO and the 16-th to 84-th percentiles in bins of CO line intensity (right panel) in our galaxies.The histogram is computed as the sum of CO line intensity in each bin of R Ci/CO , representing the amount of gas at each R Ci/CO .The NGC 4321 histogram is lower than others, showing that it has the least total molecular gas mass among the four galaxies within mapped areas.The NGC 3627 histogram has a peak about 0.5 dex higher than that of NGC 4321, but their peak locations are consistently at log R Ci/CO ∼ −1, same conclusion as our previous scatter plots show.The two SB galaxies NGC 1808 and NGC 7469, despite having much smaller observed areas (∼ 1-2 kpc), have more than an order of magnitude higher peaks hence total molecular gas masses.The locations of their histogram peaks also shift to a systematically higher R Ci/CO .It is clearly seen that both SF and SB disks lack a high-R Ci/CO (log R Ci/CO > −0.35) gas component which dominates the gas in the strong AGN host galaxy NGC 7469.
The curve-of-growth in the left panel of Fig. 5 further shows the cumulative CO brightness in bins of R Ci/CO , and reflects the fact that the majority of gas in our galactic disks has a narrow distribution of R Ci/CO (see also last few rows in Table 1).In NGC 3627, NGC 4321 and NGC 1808, the higher-and lower-  R Ci/CO gas outside ±0.1 dex around the mean R Ci/CO only sums up to < 5-10% of their total gas masses.
In the right panel of Fig. 5, we illustrate how the mean R Ci/CO changes with different gas surface density, as well as the scatter of the data point distribution.The 2D pixels are binned by CO brightness, and the mean and 16%/50%/84%-th percentiles of the R Ci/CO are computed for each bin.It sketches out the mean trends in Fig. 4 (right panel), and also reveals that a) the scatter of R Ci/CO is the largest (±0.20 dex) in the starburst disk of NGC 7469 compared to other galaxies; b) the scatter is the smallest (±0.06 dex) in the star-forming disks of NGC 4321 and NGC 3627 where Σ mol ∼ 100-500 M pc −2 ; c) the scatters are very similarly (±0.10 dex) at a relatively high Σ mol ∼ 1000 M pc −2 in all our galaxies except for NGC 4321's centre.

R Ci/CO versus ISRF
Here we study how R Ci/CO correlates with the interstellar radiation field (ISRF) intensity (U; in units of Milky Way mean ISRF; Habing 1968), which is a measure of the UV radiation field strength impacting the ISM.Since C i originates from PDRs, its abundance should be highly related to the PDR properties, thus whether or how R Ci/CO correlates with the ISRF that illuminates the PDR is a key question.
The ISRF intensity U is usually measured by the re-radiated dust emission, as the original UV photons are highly attenuated by dust in the ISM.Draine & Li (2007, hereafter DL07) developed a series of dust grain models and synthesized SED templates that can describe nearby galaxies' observed dust SEDs.A number of studies (e.g., Draine et al. 2007;Galliano et al. 2011;Dale et al. 2012;Aniano et al. 2012Aniano et al. , 2020;;Daddi et al. 2015;Chastenet et al. 2017;Liu et al. 2021;Chastenet et al. 2021) have used the DL07 models to fit the dust SEDs of galaxies near and far.For our study, we adopt the mass-averaged ISRF U maps from an ongoing effort by J. Chastenet et al. (in prep.) which is a continuation of the work published for other galaxies by Chastenet et al. (2017Chastenet et al. ( , 2021) ) with the same method.In this method, the Herschel far-infrared, Spitzer and/or WISE near-infrared images 5 are first PSF-matched to a common resolution ∼ 18 (Herschel SPIRE 250µm), then SED fitting is performed using the DL07 models at each pixel.Each dust SED is fitted by two dust components, one in the ambient ISM with a constant ISRF intensity U min , and the other in regions with higher U with a power law distributed ISRF from U min up to U max = 10 7 .The fitting then gives the best-fit U min and the power-law index γ of the U distribution.The mean ISRF intensity U is calculated as the dust mass-weighted mean of the ambient and power-law ISRF 5 For the two (U)LIRGs, the WISE maps are saturated, so their inferred ISRF intensity U should be considered with a large uncertainty.
intensities (in DL07, U PDR ≥ 10 2 is considered to be associated with the PDR), and better represents the overall ISRF than U min .
Fig. 6 presents the U maps and the R Ci/CO versus U scatter plot.Since the spatial resolution of the U maps is 18 , we calculate a weighted R Ci/CO in each circular aperture shown in Fig. 6 whose diameter equals the U map resolution, with weighting following the 2D Gaussian beam.In addition, we also show the central pixel value for each region in Fig. 6.We restrict the analysis to regions where both U and R Ci/CO are measured.
A weak correlation between R Ci/CO and U is seen from the SF to SB disks, with a large scatter.For reference, we plot a R Ci/CO = 0.07 U 0.25 line to indicate this weak trend.The lowest-U data point in NGC 3627 (mainly Region 22) and the centre of NGC 7469 do not follow this trend, possibly indicating other C i enhancement mechanisms than U .
We use the popular PDR modelling software, PDR Toolbox (version 2.2.9;Kaufman et al. 2006;Pound & Wolfire 2008,  2011) 6 to compute the theoretical R Ci/CO in uniform grids of ISRF intensity in the Habing (1968) G 0 unit and H nucleus density n H (both in log space).We adopt the latest wk2020 model set and a solar metallicity ([O/H 2 ] = 3.2 × 10 −4 , [C/H 2 ] = 1.6 × 10 −4 ).The ISRF intensity grid can be directly scaled to U by multiplying a factor 0.88 (Draine et al. 2007), corresponding to the x-axis of Fig. 6.The n H grid is not directly linked to the left y-axis of Fig. 6, but because the model contours show fairly flat R Ci/CO versus U trends, we can match the model contours and indicate the n H grid in the right y-axis of Fig. 6.
The PDR modelling predicts generally flat R Ci/CO versus U trends similar to the distribution of the observed data points, although our U are diluted values at a much larger physical scale than in the PDR models.However, in PDR models, R Ci/CO changes moderately with n H , and a relatively high n H ∼ 10 3.5-3.7 cm −3 is needed to reproduce an R Ci/CO ∼ 0.1-0.2, or equivalently ∼ 1.0-1.9 for the ratio of [C i] (1-0) and CO (2-1) fluxes in erg s −1 cm −2 sr −1 units used by the PDR Toolbox 7 .Furthermore, the PDR model needs a ∼ 0.5 dex lower n H to explain the 2× elevated R Ci/CO in SB disks, which seems somewhat counterintuitive.
This apparent tension is likely caused by the observations probing a much coarser resolution that consists of a large number of PDRs with different n H and U (e.g., increasing U from 10 to 400 can increase R Ci/CO from 0.2 to 0.5 based on our PDR Toolbox modelling).It may also relate to an abundance variation driven by non-PDR mechanisms -X-ray dominated region (XDR; Maloney et al. 1996;Meijerink & Spaans 2005;Meijerink et al. 2007; see also review by Wolfire et al. 2022) and Cosmic Ray Dominated Region (CRDR; Papadopoulos 2010; Papadopoulos et al. 2011Papadopoulos et al. , 2018;;Bisbas et al. 2015Bisbas et al. , 2017) ) are well-known to play an important role in starbursting and AGN environments, respectively.
As pointed out by Papadopoulos et al. (2011) and later works, ultra-luminous infrared galaxies (ULIRGs) with IR luminosities L IR, 8−1000µm ≥ 10 12 L and corresponding SFRs > 1000 M yr −1 have cosmic-ray (CR) energy densities ∼ 10 3 -10 4 × larger than the Milky Way.In our sample, the two SB galaxies are not as intensively star-forming as ULIRGs, but their SFRs are still a factor of 2 to 20 higher than our SF galaxies.More importantly, the SB disks have considerably higher (e.g., one to two orders of magnitude) SFR surface densities than the SF disks, plausibly leading to a factor of a few tens to hundreds that the multiplication of CI_609/CO_43 and CO_43/CO_21 does not produce the same result as the direct CI_609/CO_21 ratio in the new PDRToolbox v2.2.9.We use the new version for this work.
10 100 1000 100 1000 10000 higher CR intensities.The CRs emitted from OB star clusters and supernova remnants can penetrate deeply into the ISM, and thus nearly uniformly illuminate the ISM and effectively dissociate CO into C i.
A detailed simulation by Bisbas et al. (2017) showed that 10× and 100× higher CR ionization rates (ξ CR ) than the Milky Way level results into a factor of about 2 and 3 mildly-increased C i abundance and about 1/2 and 1/30 strongly-decreased CO abundance, respectively (see their Fig.3), in a 10 pc GMC with a mean gas density n H = 760 cm −3 , resembling Milky Way and typical SF galaxies' ISM conditions.Our SB galaxies probably have a ξ CR around or slightly above 10, but not reaching 100 given their mild U ∼ 20 − 30 compared to the U 50 in extreme SB galaxies, e.g., local ULIRGs and high-redshift hyperluminous IR-bright (HyLIRG) and submm galaxies (Daddi et al. 2015;Silverman et al. 2018;Liu et al. 2021).In Sect.4, we show that the mildly-increased [C i/CO] abundance ratio from a CRDR and temperature/density-driven excitation can lead to the factor of 2 increase in R Ci/CO when going from SF to SB disks.

R Ci/CO versus Metallicity
Our current sample covers a relatively small range in metallicity for a metallicity-line ratio study.Given their stellar masses and metallicity maps (e.g., Kreckel et al. 2019Kreckel et al. , 2020;;Williams et al. 2022), our targets probe a near-solar metallicity.NGC 3627 has a well-measured metallicity gradient of Z = 8.53+0.06×[r/R 25 ] from Kreckel et al. (2019) 8 , but the metallicity of individual H ii regions scatter between ∼ 8.45-8.65. Williams et al. (2022) produced the 2D spatial distributions of gas-phase metallicity in PHANGS-MUSE galaxies, including NGC 3627 and NGC 4321, by applying the Gaussian Process Regression to map the smooth, higher-order metallicity variation from the individual, sparse H ii regions detected in PHANGS-MUSE (Emsellem et al. 2022;Santoro et al. 2022).This method is significantly superior to a naïve nearest neighbor interpolation method.We refer the reader to Williams et al. (2022) for the technical details and validation demonstration.
In Fig. 7 we present the trend between R Ci/CO and metallicity in NGC 3627 and NGC 4321.Due to our limited metallicity range, only a very weak decreasing trend can be seen.A similar trend has been reported by Bolatto et al. (2000b), who observed CO(1-0) and [C i] (1-0) in the Region N27 of the Small Magellanic Cloud and combined with literature data to present a correlation between R Ci/CO and metallicity over ∼ 0.2-2 Z .We overlay the Bolatto et al. (2000b) empirical fitting 9 as well as a theoretical PDR model prediction line from Bolatto et al. (1999) in Fig. 7.Although our data covers only a small metallicity range, it is encouraging to see the overall good agreement with the observationally-driven trend in Bolatto et al. (2000b).More detailed understanding of the deviation from PDR models will likely need more C i mapping in lower-metallicity galaxies.

Estimation of abundance ratios under representative ISM conditions
The observed R Ci/CO depends on not only the C i and CO abundances ([CO/H 2 ] and [Ci/H 2 ]), but also the ISM conditions that determine the excitation (level population) and radiative transfer of the C atoms and CO molecules -i.e., the gas kinetic temperature T kin , volume density of H 2 (n H 2 ; as the primary collision partners of C and CO), turbulent line width ∆v, and the species' column densities N CO and N Ci .The column density divided by the line width (N/∆v) further determines the optical depth τ of each line.2021).Right panel shows the scatter plot of R Ci/CO versus U at a spatial resolution of ∼ 18 in the manually-selected apertures shown in the left panels.We computed a beam-weighted mean R Ci/CO and a centre pixel R Ci/CO for each aperture, displayed as the solid and open symbols, respectively.Error bars are the corresponding beam-weighted mean or central pixel R Ci/CO uncertainties in each region.The region diameter (∼ 18 ) corresponds to the coarsest resolution of the Herschel data used in the SED fitting.The gray dotted line shows R Ci/CO = 0.07 U 0.25 as a guiding line.The colored dotted contours (labeled 0.1-0.5)are the model contours of R Ci/CO computed from the PDR Toolbox (using the latest wk2020 model and solar metallicity; Kaufman et al. 2006).We convert G 0 from the model contour to the x-axis U by multiplying with a factor of 0.88 (Draine et al. 2007).The right y-axis indicates the PDR model's n H 2 , whose range is manually adjusted so that the model contours roughly match the left x-axis.The PDR model contours show little dependence of R Ci/CO on U in the horizontal direction but a strong correlation with n H 2 as seen by the vertical gradient.
In this work, having only a CO(2-1) and a [C i] (1-0) line is insufficient to numerically solve all the T kin , n H 2 , N CO /∆v, [CO/H 2 ] and [Ci/H 2 ].Therefore, we assume several representative ISM conditions where T kin , ∆v and [CO/H 2 ] are fixed to reasonable values, then solve the N CO and N Ci .We assume that C i and CO are spatially mixed at our resolution, i.e., they have the same filling factor.This is a relatively strong assumption and is likely not the real case, because although we see spatial colocalization at our ∼ 200 pc resolutions, the spectral profiles of CO and C i show discrepancies (Appendix D), indicating different underlying distributions at smaller scales.With the assumption of the same filling factor (and 3D distribution), the column density ratio N Ci /N CO equals the [C i/CO] abundance ratio.
We perform both LTE and non-LTE calculations.LTE assumes that all transitions are thermalized locally, and the level population follows Boltzmann statistics, so that each transition's excitation temperature T ex as defined by the level population equals T kin (e.g., Spitzer 1998;Tielens 2010;Draine 2011).The LTE assumption facilitates the computation as we do not need to solve the level population and thus n H 2 is not used.However, LTE is not easily reached in typical ISM conditions.Subthermalized lines need non-LTE calculations.The commonly adopted non-LTE approach is the Large Velocity Gradient (LVG; Scoville & Solomon 1974;Goldreich & Kwan 1974), where an escape fraction is used to facilitate the local radiation intensity calculation hence the detailed balance of radiative and collisional (de-)excitations.We use RADEX (van der Tak et al. 2010) to compute the non-LTE C i and CO line fluxes for each fixed pa-rameter set (T kin , n H 2 , N CO , N Ci , ∆v), then study how [C i/CO] determines the line flux ratio R Ci/CO .
Our choice of T kin considered multi-line measurements in the literature, e.g., Heyer & Dame (2015) for Galactic GMC ISM condition, and Teng et al. (2022) for NGC 3627 nucleus, as well as typical values in simulations (e.g., Offner et al. 2014;Glover et al. 2015;Glover & Clark 2016;Clark et al. 2019;Hu et al. 2021).Our ∆v is based on the observed line width in the highest-resolution ALMA data, i.e., from the original PHANGS-ALMA 1 data (Leroy et al. 2021b) or from Salak et al. (2019) and Izumi et al. (2020) for NGC 1808 and NGC 7469, respectively.The observed line width should at least provide an upper limit, if not directly probing the turbulent ∆v.The [CO/H 2 ] is more like a fiducial value and does not obviously affect later [C i/CO] results.We also take a fiducial Σ mol which matches the observed CO(1-0) line brightness temperature and a CO-to-H 2 conversion factor α CO ∼ 4-0.8 M (K km s −1 ) −1 depending on the assumed conditions (e.g., Bolatto et al. 2013;and our paper II).This Σ mol directly converts to an N H 2 if assuming that 36% of the mass is helium.Then N H 2 and [CO/H 2 ] together indicate a N CO .Only with these assumptions can we unambiguously link an abundance ratio [C i/CO] (i.e., N Ci /N CO ) to a line flux ratio R Ci/CO .Additionally, with the same filling factor assumption, R Ci/CO is also equivalent to the ratio of the CO-to-H 2 and C i-to-H 2 conversion factors, i.e., R Ci/CO = α CO /α Ci , because they trace the same molecular gas at our resolution.In a companion paper (paper II; Liu et al. 2022 in prep.), we study the behaviors of α CO and α Ci in more details, also under both LTE and non-LTE conditions.

LTE results
We determine the LTE solution of [C i/CO] in Fig. 8.To illustrate the complexity of this problem and the degeneracy of the T kin , N CO /∆v and [C i/CO] parameters, we show how R Ci/CO (y-axis) varies with N Ci /N CO (x-axis) and N CO (color bar) in the upper panels, and how R Ci/CO (contour) can be pinpointed by determining the Σ mol (right y-axis).At the top of each column the key parameters of the aforementioned representative ISM conditions are shown as text.The assumed T kin and ∆v generally increase from panels (1) to (7) mimicking real conditions.Observationaldriven I CO (1-0) ranges and the assumed Σ mol for these ISM environments are also listed, and can be converted one into the other via α CO .We also list the chosen [CO/H 2 ] for clarity.
As shown in the upper panels, R Ci/CO increases with both N Ci /N CO and N CO under low-T kin conditions.However, a maximum R Ci/CO is reached in the high C i column density regime, where an increasing abundance ratio will no longer increase R Ci/CO .The line ratio saturates at about 0.5 when T kin ∼ 15 K, or at about 1 for T kin ∼ 100 K.The observed line ratios (horizontal dashed lines) are all lower than the saturation limit of R Ci/COmax ∼ 0.7-0.9, and corresponds to a range of N CO or Σ mol .This range is very wide in the low-T kin conditions, with log N CO /cm −2 ∼ 17 − 19, but is very narrow at T kin ∼ 100 K, that is, R Ci/CO can be determined from [C i/CO] without much dependence on N CO .
In the lower panels, the contours of R Ci/CO = 0.05, 0.1, 0.2 and 0.5 are shown in the plane of CO column density N CO versus C i/CO abundance ratio.For a fixed R Ci/CO , a higher N Ci /N CO corresponds to a lower N CO in the low-T kin regime, but R Ci/CO becomes insensitive to N CO at a high T kin .In the latter case, N Ci /N CO alone determines R Ci/CO .By matching to the expected Σ mol for each ISM condition, we can obtain an unambiguous N Ci /N CO from the observed R Ci/CO as shown in Fig. 8.
The strongly enhanced C i/CO abundance ratios in SB disk and AGN environments are clear signatures of CRDR and XDR, respectively.The NGC 4321 centre is a low-ionization nuclear emission-line region (LINER, e.g., García-Burillo et al. 2005), and given our assumptions it requires N Ci /N CO ∼ 0.4 in order to explain the observed low line ratio.The gradually increasing C i/CO abundance ratios in the nuclei of NGC 4321, NGC 3627, NGC 1808 and NGC 7469 agree well with the increasing power of AGN activities, e.g., their X-ray luminosities (Table 1).
Our assumptions for the Galactic cloud ISM condition are mainly based on Heyer et al. (2009), Heyer & Dame (2015) and Roman-Duval et al. (2010).The obtained [C i/CO] ∼ 0.1 also agrees well with previous studies using C i and CO isotopologues (see Sect. 1.For example, Ikeda et al. (1999) reported N Ci /N CO ∼ 0.05-0.2for the Orion cloud; Oka et al. (2001) found N Ci /N CO = 0.064 ± 0.035 but a peak of N Ci /N CO ∼ 0.15 in the DR 15 H ii region and two IR dark clouds; Kramer et al. (2004) found N Ci /N CO ∼ 0.11-0.12 in the massive SF region W3; Genzel et al. (1988) found N Ci /N CO ∼ 0.2 for the massive SF region W51; Zmuidzinas et al. (1988) reported N Ci /N CO ∼ 0.05-0.15 in several dense clouds.Recently, Izumi et al. (2021) derived an overall N Ci /N CO ∼ 0.1 in the massive SF region RCW38 at A V ∼ 10-100, consistent with the Orion cloud, but also saw high [C i/CO] ∼ 0.2-0.6 in some even lower column density regions.
Accurate determination of all the ISM properties instead of simple assumptions as in this work requires comprehensive, multi-J CO isotopologue observations (e.g., Teng et al. 2022).Our representative ISM conditions only provide a qualitative diagnostic to understand how to infer the C i/CO abundance ratio from the observed R Ci/CO quantity.Their properties may not be very accurate, and are subject to change with future multitracer/CO isotopologue studies.

Non-LTE results
We perform the non-LTE calculations using the RADEX code (van der Tak et al. 2007) with the Large Velocity Gradient (LVG; i.e., expanding sphere) escape probability (Scoville & Solomon 1974;Goldreich & Kwan 1974).Comparing to LTE, the volume density of H 2 is now an additional property that plays a role.The excitation temperature of each line no longer equals T kin .We compute the non-LTE line ratios in a grid of (T kin , ∆v, [CO/H 2 ], N Ci /N CO , N CO , n H 2 ) parameter sets, then study the line and abundance ratios similarly as in the previous section.
In Fig. 9 we show the non-LTE version of Fig. 8.The non-LTE result does not qualitatively differ from that of LTE.The most obvious difference is that at high temperature (T kin > 50 K, given our other assumed conditions), the required [C i/CO] is lower by ×2 for the same observed R Ci/CO , whereas at low temperature (T kin ∼ 15 K) a 2× [C i/CO] is expected for the observed line ratio.
As listed in the title of each panel, we assume log n H 2 /cm −3 = 2.5 to 3.0 for our conditions (1) to (7).This is based on the fact that low-J CO and C i trace the bulk of molecular gas with log n H 2 /cm −3 ∼ 2-3 (e.g., Leroy et al. 2017;Liu et al. 2021), with a typical value n H 2 ∼ 300 cm −3 widely adopted/found for Galactic and nearby galaxy GMCs (e.g., Koda et al. 2012;Pety et al. 2017, see also simulations of Glover & Clark 2012), whereas a smaller portion of their emission is from the dense gas with log n H 2 /cm −3 ∼ 4-6 (e.g., Gao & Solomon 2004a).The fraction of the dense gas is expected to increase in more starbursting environment (Gao & Solomon 2004b), so does the mean gas density (Liu et al. 2021).Therefore, a log n H 2 /cm −3 of 3.0 is perhaps a reasonable assumption, or at least a lower limit for our galaxy centres.A 10× higher density does not qualitatively change Fig. 9 except for shifting the expected N Ci /N CO to be 50% higher in the T kin > 50 K conditions (making Fig. 9 more like Fig. 8).Therefore, our non-LTE result strongly confirms the LTE result in the previous section, that the intrinsic C i/CO abundance has to vary significantly from about 0.1-0.2 to about 1-3, in order to explain the observed line ratios from Galactic GMCs (hundred pc scales) and NGC 3627/4321 SF disk GMCs to NGC 1808/7469 SB disk/AGN (r 250 pc) environments.

Discussion
From this study we see that the line ratio R Ci/CO increases from 0.1 to 0.2 and > 0.5 from the SF disk to the SB disk and the vicinity (few hundred pc) of X-ray luminous AGN.Although we may be limited by our sample, we suggest that this trend can be extended to general cases at both local and high redshift, because the C i/CO abundance ratio enhancement is well-predicted by the CRDR and XDR theories.Papadopoulos (2010) and Papadopoulos et al. (2011) pointed out the key role of CRDR in local (U)LIRGs.With UV photons alone, cold gas in starbursting environments can still be well shielded, hence as cold as ∼ 10 K.However, with CRs produced from O, B star clusters and supernova remnants, the gas can be uniformly illuminated and heated up.The extremely starbursty local ULIRGs have a CR density (ξ CR ) 10 3 -10 4 times that of the Milky Way (Papadopoulos 2010), and thus will naturally dissociate more CO into C i. Bisbas et al. (2017) and Papadopoulos et al. (2018) conducted simulations of molecular gas at tens of pc scale with varying CR strength, and found that when ξ CR is ∼ 1-100 × Galactic value, the gas temperature in the dense gas remains largely unaffected (∼ 10 K).However, when ξ CR reaches ∼ 1000 × Galactic value as in (U)LIRGs, the gas temperature increases to 30-50 K.The Bisbas et al. (2017) simulation also shows that the C i and C ii fractional abundances at log A V ∼ 0.5 increase by one and three orders of magnitude, respectively, and the CO fractional abundance decreases correspondingly by over two orders of magnitude, with an increasing ξ CR from the Galactic value to 100 × of that.These results are in line with our calculations in Sect.4.1.
Meanwhile, the AGN-driven XDR chemistry and molecule/atom excitation have been firstly modeled in great detail by Maloney et al. (1996), Meijerink & Spaans (2005) and Meijerink et al. (2007).In XDRs, photo-ionization becomes the dominant heating source instead of photo-electric emission from dust grains in PDRs, thus more molecules are dissociated and/or ionized.Meijerink & Spaans (2005)'s Model 1 shows that CO molecules are fully dissociated when N H < 10 23 cm −2 , then the fractional abundance is gradually increased by six orders of magnitudes from N H = 10 23 cm −2 to 10 24 cm −2 .By contrast, the C i and C ii abundances are about constant until reaching the densest part (N H > 10 24.5 cm −2 ).Unlike in PDRs, there is no clear C ii/C i/CO stratification in XDRs (Meijerink et al. 2007;Wolfire et al. 2022).Their Model 2 with a 100× higher energy deposition rate shows that the CO distribution contracts to the highest gas density spots, and the total cumulative [C i/CO] is increased by about 0.2 dex, reaching N Ci /N CO ∼ 1.6.This value agrees well with our calculations shown in Figs. 8 and 9 (the ISM condition 7).
Therefore, both in theoretical and observational views, it is clear that the R Ci/CO line ratio is a good indicator of the starburst's CRDR and AGN's XDR.If a R Ci/CO line ratio (to be specific, I [C i] (1-0) /I CO (2-1) ) reaches 0.2-0.5, then the CRDR plays a crucial role and the abundance ratio [C i/CO] could be ∼ 1.If it exceeds 0.5-1.0,then it may only be enhanced by an XDR, e.g., caused by an X-ray luminous AGN, which boosts the abundance ratio [C i/CO] to 1-2.

Conclusion
We presented new beam-and uv-coverage-matched [C i] (1-0) and CO(2-1) maps in the galactocentric radius r 7 kpc SF disk of NGC 3627 at ∼ 190 pc resolution, and in the r 3 kpc disk of NGC 4321 at ∼ 270 pc resolution.They are among the largest C i mosaic mappings available in nearby galaxies.
We combined the imaging of [C i] (1-0) and CO(2-1) in two nearby, more starbursty and AGN-hosting galaxies, i.e., NGC 1808 central ∼ 1 kpc, and NGC 7469 central ∼ 1.5 kpc at 140-160 pc resolutions.Together, we studied the spatial distributions of C i and CO, and their line ratio R Ci/CO as functions of various galaxy resolved properties, including CO brightness (Fig. 4), ISRF U (Fig. 6), and metallicity (Fig. 7).
In order to obtain the underlying intrinsic [C i/CO] abundance ratio from the observed line ratio R Ci/CO , we assumed 7 representative ISM conditions and did both LTE and non-LTE calculations.
We summarize our findings as below: -The majority of the molecular gas in the SF disks of NGC 3627 and NGC 4321 has a uniform R Ci10/CO21 = 0.10 ± 0.05.[C i] (1-0) and CO (2-1) exhibit very similar spatial distributions at our 140-270 pc resolution.
-The majority of the molecular gas in the SB disk of NGC 1808 has R Ci10/CO21 = 0.20 ± 0.05.This elevated line ratio compared to the SF disk is likely caused by a ∼ 3-5× increased [C i/CO] abundance ratio, based on our LTE and non-LTE calculations and parameter assumptions (Figs. 8  and 9).-The centres (r 250 pc) of NGC 4321, NGC 3627, NGC 1808 and NGC 7469 show a strongly increasing trend in R Ci10/CO21 , from 0.05 to 0.5, as a function of I CO (2-1)i.e., R Ci10/CO21 ≈ 0.035 (I CO (2-1) /[100 K km s −1 ]) 0.8 .The inferred [C i/CO] abundance ratio from our non-LTE calculation changes from ∼ 0.2 to ∼ 3, in line with the increasing X-ray luminosity of the centres/AGNs.
-The observed R Ci/CO versus ISRF trend is fairly flat, i.e., R Ci/CO is insensitive to ISRF intensity changes, similar to the PDR model prediction calculated with PDRToolbox (although our resolution is much coarser than the physical scale of PDR models; Sect.3.3).-We find a R Ci/CO versus metallicity trend consistent with previous observations that covered lower-metallicity environments (∼ 0.2 Z ; Bolatto et al. 2000b), within the metallicity range (∼ 0.7-1 Z ) of our sample (Sect.3.4).-Our inferred [C i/CO] abundance ratios agree well with the CRDR and XDR's theoretical studies.Given the drastically changed R Ci/CO , we propose that the C i to CO line ratio can be a promising indicator of CRDR (starbursting) and XDR (AGN).For instance, an observed line flux ratio of R Ci10/CO21 ∼ 0.2-0.5 and 0.5 effectively indicates CRDR and XDR, respectively.-We do not find ubiquitous "CO-dark", C i-bright gas at the outer r ∼ 1-3 kpc disks of NGC 3627 and NGC4321, nor systematic enhancement in R Ci/CO at CO-faintest sightlines via stacking, mainly due to the incompleteness of our data.Nevertheless, a few sightlines, e.g., Region 22 in NGC 3627 show high R Ci/CO ∼ 0.5 which do not fully comply with the noise behavior seen in our CASA simulations.
factor of 0.2 in Kelvin units for [C i] (1-0) (representing a constant R Ci/CO = 0.2).Note that we do not use the short-spacingcorrected PHANGS-ALMA CO(2-1) cube because it has a finite beam size of ∼ 1 , and simulation based on that may result in more extended emissions.On the contrary, the clean model cube has an infinite resolution, and because both "multiscale" and "singlescale" cleaning are used in the PHANGS-ALMA imaging, the clean model contains point as well as extended sources.
We also set the actual observing dates and integration times as the real observations.Then we run CASA simobserve to simulate the visibilities of each observing block, concatenate them together, and process them into image cubes and moment maps in the same way as described in Appendix A. Meanwhile, by convolving the input skymodel with a 2D Gaussian beam corresponding to our working resolution, and integrating over the spectral axis, we obtain a "true" skymodel moment-0 map which can be directly compared to the imaged moment-0 maps from the simulation.Their ratio map shows a missing flux of ∼ 50% to ∼ 7% in the strict mask, and more specifically, 24%, 58%, 85% and 96% pixels selected in our method as for real data have a missing flux < 20%, < 30%, < 40% and < 50%, respectively.Fig. B.1 illustrates the results from methods (b) (upper panels) and (c) (lower panels).Left panels show the CO (2-1) moment-0 maps, analyzed at the same working resolution as our main analysis in Appendix A. The simulated data show a fainter emission because we used the clean mode cube as the skymodel, which was made by cleaning the PHANGS-ALMA 12m+7m data down to 1 σ within the clean mask voxels.This will indeed lose some fluxes.However, the right panels show that the two methods lead to consistent missing flux results.
Note that in the last panel of Fig. B.1, some regions appear black in the simulation-based analysis because we require line ratios to be analyzed in the strict mask, while these regions are not in it anymore.This is because the skymodel loses some fluxes as mentioned above, thus the strict mask is smaller, and also because in the simulation we assumed a constant R Ci/CO which may lead to fainter C i than real data along faint sightlines, e.g., Region 22, thus further shrinking the strict mask where both lines need to have sufficient S/N.
Bearing these in mind, we conclude that methods (b) and (c) have consistent results on the missing flux, and method (a) is less accurate due to the unconsidered uncertainties.
These checks show that the overall missing flux in NGC 3627 is 31% based on the full-and trimmed-uv CO (2-1) analysis, or ∼ 36 − 88% based on the more-uncertain Herschel data.At the pixel level, it varies from 7% to ∼50% from bright to faint sightlines at our working resolution of 3.46 .These values are as expected given our ACA-only observation setup.We emphasize again that by applying the uv-matching in our data reduction, our C i/CO line ratio study is not obviously affected by the missing flux issue, except perhaps in the faintest pixels or "cloud edges".

Fig. 1 .
Fig. 1.Top panels: NGC 3627 [C i] (1-0) (left) and CO (2-1) (right) integrated line intensity (moment-0) maps within our broad mask (representing high-completeness signals; see Appendix A).Bottom left panel: [C i] (1-0)/CO (2-1) line ratio map within the combined signal mask.Bottom right panel: HST F814W-F555W-F438W RGB composite image (PHANGS-HST; Lee et al. 2022).Regions marked with the white boxes and labeled with numbers are mainly for illustration purpose and for discussion in the text.Their zoom-in images and extracted spectra are provided in our online data release.Examples for Regions 1, 4 and 22 are presented in Appendix D. All images have the same field of view of 179 × 179 , or 9.83 kpc × 9.83 kpc.

Fig. 3 .
Fig. 3. [C i] (1-0) (left) and R Ci/CO (right) maps of NGC 1808 in the upper panels and NGC 7469 in the lower panels.Similar to Figs. 1 and 2, the line intensity maps are defined in the broad mask and the R Ci/CO maps are restricted to the combined signal mask.Boxes in the left panels indicate the manually selected central regions where the pixels are highlighted in Fig. 4. Contours in the right panels indicate the 5, 10, 20 and 50-σ levels of the left-panel C i intensity maps.

Fig. 4 .
Fig. 4. Left: I [C i](1-0) versus I CO in the four galaxies we studied: NGC 3627, NGC 4321, NGC 1808 and NGC 7469.I [C i] (1-0) and I CO (2-1) integrated fluxes are moment-0 intensities extracted from their beam and pixel size matched data cubes.Only data points with propagated moment-0 S/N ≥ 5 are shown.The blue (with crossed hatch), green (with horizontal line hatch), gray (with dotted hatch) and light gray shadings indicate the 5 σ thresholds for the four galaxies, respectively.Typical line widths of 10 km s −1 (NGC 3627 and NGC 4321) and 30 km s −1 (NGC 1808 and NGC 7469) in sigma are used to calculate the thresholds together with the channel rms in Table1.For NGC 3627 and NGC 4321, their (uvmatched) CO(2-1) 5 σ thresholds are out of the plotting range.Top axis shows the simply-converted H 2 surface density using a constant α CO = 4.3 and CO excitation (i.e., CO(2-1) to CO(1-0) line ratio) R 21 = 0.65.Right: Same data but shown as the line ratio R Ci/CO = I [C i] (1-0) /I CO (2-1) versus I CO (2-1) .The blue horizontal line indicates a star-forming (SF) disk-like R Ci/CO ∼ 0.1, i.e., the spiral arms of NGC 3627 and NGC 4321.The gray horizontal line indicates a starburst (SB) disk-like R Ci/CO ∼ 0.2, i.e., the ∼ 1 kpc SB disk/ring in NGC 1808 and NGC 7469.The red line with a slope of 0.8 in logarithm (0.035 (I CO (2-1) /[100 K km s −1 ]) 0.8 ) visually guides the trend of increasing R Ci/CO in the central regions of these galaxies (see boxes in Figs. 1 to 3, sizes ∼ 300-500 pc).

Fig. 5 .
Fig. 5. Left panel shows the histograms of the sum of CO(2-1) surface brightness over pixels in each C i/CO line ratio interval in the four galaxies.The bold curves are the cumulative distributions of the corresponding histograms.Right panel shows the mean C i/CO line ratio for pixels in each bin of CO(2-1) surface brightness (or molecular gas surface density as indicated by the top-axis assuming a constant α CO = 4.3 and R 21 = 0.65).Shaded area enveloping each mean trend indicates the 16-to 84-th percentiles in each bin, and dashed line indicates the 50-th percentile.This panel thus more clearly illustrates the ridge and scatter of the trends in the right panel of Fig. 4. Detection limits are shown consistently with Fig. 4.

Fig. 6 .
Fig.6.Line ratio versus mean ISRF intensity U measured from SED fitting to the near-to far-IR(Spitzer, Herschel)  data.Left four panels are the U maps of NGC 3627, NGC 4321, NGC 1808 and NGC 7469, and as described in Sect.3.3 are obtained from the works done by J.Chastenet et  al. (in prep.)  as a continued effort ofChastenet et al. (2021).Right panel shows the scatter plot of R Ci/CO versus U at a spatial resolution of ∼ 18 in the manually-selected apertures shown in the left panels.We computed a beam-weighted mean R Ci/CO and a centre pixel R Ci/CO for each aperture, displayed as the solid and open symbols, respectively.Error bars are the corresponding beam-weighted mean or central pixel R Ci/CO uncertainties in each region.The region diameter (∼ 18 ) corresponds to the coarsest resolution of the Herschel data used in the SED fitting.The gray dotted line shows R Ci/CO = 0.07 U 0.25 as a guiding line.The colored dotted contours (labeled 0.1-0.5)are the model contours of R Ci/CO computed from the PDR Toolbox (using the latest wk2020 model and solar metallicity;Kaufman et al. 2006).We convert G 0 from the model contour to the x-axis U by multiplying with a factor of 0.88(Draine et al. 2007).The right y-axis indicates the PDR model's n H 2 , whose range is manually adjusted so that the model contours roughly match the left x-axis.The PDR model contours show little dependence of R Ci/CO on U in the horizontal direction but a strong correlation with n H 2 as seen by the vertical gradient.

Fig. 7 .
Fig. 7. Left two panels show the metallicity maps of NGC 3627 and NGC 4321 from Williams et al. (2022) based on PHANGS-MUSE (Emsellem et al. 2022) data.Boxes are regions as defined in Figs. 1 and 2. Right panel shows the R Ci/CO versus metallicity scatter plot for regions in NGC 3627 and NGC 4321.The bottom-axis is metallicity in solar units, top-axis is metallicity in 12 + log(O/H), left-axis is the line ratio of [C i] (1-0) and CO (2-1) in K km s −1 units, and right-axis is the ratio of fluxes in erg s −1 cm −2 sr −1 units.The left-and right-axes are matched by a factor of (ν Ci1-0 /ν CO1-0 ) 3 × R 21 , where we adopt a R 21 = 0.65 (Leroy et al. 2022).Data points are the mean metallicity and R Ci/CO in each region shown in left panels, with error bars indicating the typical uncertainty in R Ci/CO , and vertical color shading indicating the 16-to 84-th percentiles of pixels' R Ci/CO distribution inside each region.The Bolatto et al. (2000b) empirical fitting and Bolatto et al. (1999) model prediction lines are overlaid.The three outliers in NGC 3627 with R Ci/CO ∼ 0.19, 0.28 and 0.35 are Regions 20, 21 and 22, respectively.The outlier point in NGC 4321 with R Ci/CO ∼ 0.06 is its centre Region 1.

Fig. 8 .
Fig. 8. Observed/calculated R Ci/CO as a function of CO column density N CO and [C i/CO] abundance (column density) ratio N Ci /N CO for 7 representative ISM conditions, based on LTE calculations.Upper panels: Curved lines are R Ci/CO versus N Ci /N CO for different N CO as indicated by the color.Dashed means CO(2-1) optically thick and [C i] (1-0) optically thin, dotted means both are both optically thin, and solid line means both optically thick.Horizontal thick line is the observed R Ci/CO value.Lower panels: Black contours are R Ci/CO = 0.05, 0.10, 0.20 and 0.50 in the CO column density versus C i/CO abundance ratio plane.Vertical color shading indicates the final N Ci /N CO that can match both the observed R Ci/CO and the assumed Σ mol value noted at the top of each panel.